信息工程概论-窗口傅里叶变换

第1页共9页信息工程概论作业——窗口傅里叶变换姓名：白子轩学号：2130602008班级：信计31一、传统的傅里叶变换我们都知道，信号分析中最重要的两个参数是时间和频率，而我们一般所得到的信号表示形式都是tf的形式，而我们可以通过传统的傅里叶变换dtetfftj)()(ˆ，可以把信号变为频域表示）（f。但是，传统的傅里叶变换只对平稳的信号有用。对于非平稳的信号需要用时间和品率的联合函数来表示信号。因此，我们需要短时傅里叶变换，也就是窗口傅里叶变换。二、窗口傅里叶变换对于信号的频率是随时间变化的信号。为了获得它的随时间变化的频谱，最采用的处理办法是加窗技术对信号截取，然后对截取的局部信号作Fourier变换。然后不断地移动窗口函数中心的位置，就可以得到信号的局部区域的瞬时频率，因此，对于连续的信号，它的窗口傅里叶变换为：(,)()()itSfuftgtuedt窗口傅里叶逆变换为：1()(,)()2itftSfugtueddu而对于离散的信号，它的窗口傅里叶变换为：102,expNnilnSfmlfngnmN窗口傅里叶逆变换为：110012,expNNmlilnfnSfmlgnmNN三、窗口函数要进行窗口傅里叶变换，首先要要选择窗口函数，窗口函数有很多，例如高斯窗、hamming窗和Hanning窗等等。其中高斯窗函数被设计为了分析瞬态信号，Hamming和Hann窗函数被设计为了分析窄带信号，Kaiser-Bessel窗函数可用于更好地分离两个频率成分非常第2页共9页接近但振幅完全不同的信号。在matlab中我们也可以直接调用一些窗口函数，调用的方法如下：四、实验实验1：题目：在这里我先做了一下书上的例子，对线性调频信号2sin22ftt进行频谱分析。题目分析：这个例子有两种做法，第一种方法是直接调用matlab中的spectrogram函数，第二种方法是按照定义选取窗口函数，然后对每一小段的做快速傅里叶变换就可以了。方法一：源程序：cleart=0:0.001:10;t1=t;f1=sin(2*pi*2*power(t,2));subplot(2,2,1);plot(t,f1);subplot(2,2,2);第3页共9页g=1/6*exp(-0.5*power(t,2)).*(t=-3&t=3)+0.*(t3|t-3);t=-4:0.01:4;g1=1/6*exp(-0.5*power(t,2)).*(t=-3&t=3)+0.*(t3|t-3);plot(t,g1);subplot(2,2,3);[S,F,T,P]=spectrogram(f1,gausswin(600),580,600,1E3);surf(T,F,10*log10(abs(P)),'edgecolor','none');axis([010050])view(0,90);xlabel('Time(Seconds)');ylabel('Hz');subplot(2,2,4);surf(T,F,10*log10(abs(P))+80,'edgecolor','none');axis([0100500200])%axistightxlabel('Time(Seconds)');ylabel('Hz');zlabel('enargy')得到结果：第4页共9页结果分析：我们用的是高斯窗口，得到了一个很好的结果，无论是2D图还是3D图，都与书上的图十分相似。但有一个十分大的缺陷，就是无法重构原来的信号，因为我们是直接调用的spectrogram函数，并不太知道里面的具体程序是什么样的，所以无法还原原信号，也无法计算误差。因此我们就需要第二种方法。方法二：源程序：cleart=0:0.001:10;t1=t;f1=sin(2*pi*2*power(t,2));subplot(3,2,1);plot(t,f1);subplot(3,2,2);g=1/6*exp(-0.5*power(t,2)).*(t=-3&t=3)+0.*(t3|t-3);t=-4:0.01:4;g1=1/6*exp(-0.5*power(t,2)).*(t=-3&t=3)+0.*(t3|t-3);plot(t,g1);N=length(f1);Nw=20;%窗函数长windowlengthL=19;%窗函数每次移动的样点数,重叠宽度Ts=round((N-Nw)/L)+1;%计算把数据x共分成多少段nfft=2^ceil(log2(Nw));%FFT的长度TF=zeros(Ts,nfft);%将存放三维谱图，先清零%fori=1:Tsi=0;flag=0;whileflag==0;第5页共9页i=i+1;if(i-1)*Nw+NwN%y((i-1)*L+1:i*L+L)%hamming(Nw)xw=f1((i-1)*L+1:i*L+L)*hamming(2*L);%取一段数据temp=fft(xw,nfft);%FFT变换temp=fftshift(abs(temp));%频谱以0频为中心%TF(i,:);TF(i,:)=temp;%把谱图存放在TF中elseflag=1;endendsubplot(3,2,3);%mesh(abs(TF));%view(0,90);%axistight%imagesc(TF);contour(abs(TF));xlabel('时间');ylabel('频率')subplot(3,2,4);mesh(abs(TF));%三维绘图axistighttitle('STFT');xlabel('时间');ylabel('频率');%%%%短时傅里叶变换X=fft(f1);X=fftshift(X);第6页共9页%%%%重构x1=ifftshift(X);x1=ifft(x1);subplot(3,2,5);plot(t1,x1);xlabel('时间');%x轴ylabel('振幅');%y轴%%%%误差e=x1-f1;subplot(3,2,6);plot(t1,abs(e));xlabel('时间');%x轴ylabel('振幅');%y轴得到结果：结果分析：我们发现得到的结果非常不理想，经过多次试验，我并没有找到是哪里出的问题，因此第7页共9页在接下来的实验中将要放弃这种方法。实验2：题目：录一段自己的声音，对这段声音进行频谱分析。实验分析：这个题目和实验1类似，只不过一个用的是连续的窗口傅里叶变换，一个用的是离散的窗口傅里叶变换，所以在这里我们仍然运用matlab中的spectrogram函数即可。在这里我们需要注意一下，mp3是双音轨的，所以数字化后是一个是一个二维数组，而我们只需要分析其中的一组即可。于是我就录了自己的一段笑声，并抽取其中的一组数组，进行分析。源程序：clear[x,fs]=audioread('C:\Users\lenovo\Desktop\小波分析\笑.mp3');y=x(:,1)';subplot(3,2,1);plot(y);title('信号波形图');%图名xlabel('时间');%x轴ylabel('振幅');%y轴holdonN=length(y);Nw=600;%窗函数长windowlengthL=Nw/2;%窗函数每次移动的样点数,重叠宽度%%%%%%窗口函数subplot(3,2,2);g=hamming(Nw);j=-Nw/2;fori=1:Nw第8页共9页t(i)=j;j=j+1;endplot(t,g,'.');subplot(3,2,3);[S,F,T,P]=spectrogram(y,hamming(600),580,600,1E3);surf(T,F,10*log10(P),'edgecolor','none');axistightview(0,90);xlabel('Time(Seconds)');ylabel('Hz');subplot(3,2,4);surf(T,F,10*log10(P),'edgecolor','none');axistightxlabel('Time(Seconds)');ylabel('Hz');zlabel('enargy')Y=fft(y);X=fftshift(Y);%%%%重构x1=ifftshift(X);x1=ifft(x1);subplot(3,2,5);plot(x1);xlabel('时间');%x轴ylabel('振幅');%y轴%%%%误差e=x1-y;第9页共9页subplot(3,2,6);plot(abs(e));xlabel('时间');%x轴ylabel('振幅');%y轴得到结果：结果分析：我们可以看到得到了一个比较比较清晰地频谱图，但是它的3D图并不是太清晰。结果不是太理想。