信息工程概論-窗口傅里葉變換

1、信息工程概論-窗口傅里葉變換作者: 日期:信息工程概論作業(yè)窗口傅里葉變換姓名：白子軒學號：2130602008班級：信計31一、傳統(tǒng)的傅里葉變換而我們一般所得到的信號表我們都知道，信號分析中最重要的兩個參數(shù)是時間和頻率,示形式都是f t的形式，而我們可以通過傳統(tǒng)的傅里葉變換?()f (t)e j 'dt，可以把信號變?yōu)轭l域表示f ()。但是，傳統(tǒng)的傅里葉變換只對平穩(wěn)的信號有用。對于非平穩(wěn)的信號需要用時間和品率的 聯(lián)合函數(shù)來表示信號。因此，我們需要短時傅里葉變換，也就是窗口傅里葉變換。二、窗口傅里葉變換對于信號的頻率是隨時間變化的信號。為了獲得它的隨時間變化的頻譜，最采用的處理辦法是加窗

2、技術對信號截取，然后對截取的局部信號作Fourier變換。然后不斷地移動窗口函數(shù)中心的位置，就可以得到信號的局部區(qū)域的瞬時頻率,因此，對于連續(xù)的信號，它的窗口傅里葉變換為:Sf(u,)窗口傅里葉逆變換為:f(t)2Sf (u, )g(t u)e td du而對于離散的信號，它的窗口傅里葉變換為:Sf m,lN1f ng n m exp n 0N窗口傅里葉逆變換為:.N 1 N 1-Cl1 Ii2 Inn Sf m,l g n m exp N m 0I 0N三、窗口函數(shù)要進行窗口傅里葉變換，首先要要選擇窗口函數(shù)，窗口函數(shù)有很多，例如高斯窗、hamming窗和Hanning窗等等。其中高斯窗函數(shù)被

3、設計為了分析瞬態(tài)信號，Hamming和Hann窗函數(shù)被設計為了分析窄帶信號，Kaiser-Besse I窗函數(shù)可用于更好地分離兩個頻率成分非常接近但振幅完全不同的信號。在matlab中我們也可以直接調用一些窗口函數(shù)，調用的方法如下:(n),根據(jù)長疥r產生一個(1 矩形窗 C Rectangle %'i ndow)(2' 三角囲(.?rLangular Windov) 角S嘰w-tr iang (n),根據(jù)長度n嚴生一個(3?限寧窗(Hanning Tindow) Z寧窗W.w-liannitiE(n),根扌居長度n產生一調ffl稲式:whanmifis (n)，根據(jù)長度n產生一

4、個(4；每明窗(HfiEBiLni5 Window) 肚月茴叭調胳式!w=hl fn).桐據(jù)長度r產生一i5 * 韋竝克w窗(El ack.man Window) 卜布拉克曼窗w<(6.|愷撒'® (Kaiser Window)勾®®數(shù)旁瓣的呂參數(shù)產生一"卜愷撒餵0調用格w=kaiser(n,beta),根括長度n札影四、實驗實驗1：題目在這里我先做了一下書上的例子，對線性調頻信號2f t sin 2 2t進行頻譜分析。這個例子有兩種做法，第一種方法是直接調用matlab中的spectrogram函數(shù)，第二種方法是按照定義選取窗口函數(shù)，然后

5、對每一小段的做快速傅里葉變換就可以了。方法源程序:cleart=0:0.001:10;t1=t;f1=si n(2* pi*2* po wer(t,2);sub plot(2,2,1);plot(t,f1);sub plot(2,2,2);g=1/6*ex p(-0.5* po wer(t,2).*(t>=-3 & t<=3)+0.*(t>3 | t<-3);t=-4:0.01:4;g1=1/6*ex p(-0.5* po wer(t,2).*(t>=-3 & t<=3)+0.*(t>3 | t<-3);plot(t,g1);su

6、b plot(2,2,3);S,F,T ,P = sp ectrogram(f1,gausswi n(600),580,600,1E3);surf(T,F,10*log10(abs( P),'edgecolor'， none');axis(0 10 0 50)view(0,90););xlabel('Time (Seco nds)'ylabel( 'Hz');sub plot(2,2,4);surf(T,F,10*log10(abs( P)+80,'edgecolor'， none');axis(0 10 0 50

7、 0 200)%axis tightxlabel('Time (Sec on ds)');ylabel('Hz');zlabel('en argy' )得到結果:-10.5a0.5I'D43JTi-ne結果分析:我們用的是高斯窗口，得到了一個很好的結果，無論是2D圖還是3D圖，都與書上的圖十分相似。但有一個十分大的缺陷，就是無法重構原來的信號，因為我們是直接調用的spectrogram函數(shù)，并不太知道里面的具體程序是什么樣的，所以無法還原原信號，也無法計算誤差。因此我們就需要第二種方法。方法二:源程序:clear t=0:0.001:10

8、;t1=t;f1=si n(2*pi *2* po wer(t,2);sub plot(3,2,1);plot(t,f1);sub plot(3,2,2);g=1/6*ex p(-0.5* po wer(t,2).*(t>=-3 & t<=3)+0.*(t>3 | t<-3);t=-4:0.01:4;g1=1/6*ex p(-0.5* power(t,2).*(t>=-3 & t<=3)+0.*(t>3 | t<-3);plot(t,g1);N=le ngth(f1);Nw=20;%窗函數(shù)長 window lengthL=19;%

9、窗函數(shù)每次移動的樣點數(shù)，重疊寬度Ts=rou nd(N-Nw)/L)+1;%計算把數(shù)據(jù)x共分成多少段n fft=2ceil(log2(Nw);%FFT的長度TF=zeros(Ts ,n fft);%將存放三維譜圖，先清零%for i=1:Tsi=0;flag=O;while flag=O;i=i+1;if (i-1)*Nw+Nw<N%y(i-1)*L+1:i*L+L)%取一段數(shù)據(jù)%hammi ng(Nw)xw=f1(i-1)*L+1:i*L+L)*hammi ng(2*L);temp=fft(xw, nfft);%FFT變換temp=fftshift(abs(te mp);%頻譜以0頻為

10、中心%TF(i,:);TF(i,:)=tem p;%把譜圖存放在TF中elseflag=1;end end sub plot(3,2,3);%mesh(abs(TF);%view(0,90);%axis tight %imagesc(TF);con tour(abs(TF);xlabel('時間');ylabel('頻率')sub plot(3,2,4);mesh(abs(TF);%E維繪圖axis tighttitle( 'STFT');xlabel('時間');ylabel('頻率');%短時傅里葉變換X=ff

11、t(f1);X=fftshift(X);%重構x1=ifftshift(X);x1=ifft(x1);sub plot(3,2,5);plot(t1,x1);xlabel('時間');%x軸ylabel('振幅');%y軸%差e=x1-f1;sub plot(3,2,6);plot(t1,abs(e);xlabel('時間');%x軸ylabel('振幅);%y軸得到結果:500< 300蔗2001DOSTFT10 20時間3010”孑50(3oaor30IDS時間-15時間結果分析:我們發(fā)現(xiàn)得到的結果非常不理想,經過多次試驗,我并

12、沒有找到是哪里出的問題，因此在接下來的實驗中將要放棄這種方法。實驗2：題目 錄一段自己的聲音，對這段聲音進行頻譜分析。實驗分析:這個題目和實驗1類似，只不過一個用的是連續(xù)的窗口傅里葉變換，一個用的是離散的窗口傅里葉變換，所以在這里我們仍然運用matlab中的spectrogram函數(shù)即可。在這里我們需要注意一下，mp3是雙音軌的，所以數(shù)字化后是一個是一個二維數(shù)組，而我們只需要分析其中的一組即可。于是我就錄了自己的一段笑聲，并抽取其中的一組數(shù)組，進行分析。源程序:clearx,fs=audioread(小波分析笑.mp 3'）;y=x(:,1)'sub plot(3,2,1);p

13、lot(y);title（'信號波形圖）； 略名xlabel('時間');%x軸ylabel('振幅);%y軸hold onN=le ngth(y);Nw=600;%窗函數(shù)長 window lengthL=Nw/2;%窗函數(shù)每次移動的樣點數(shù)，重疊寬度%窗口函數(shù)sub plot(3,2,2);g=hammi ng(Nw);j=-Nw/2;for i=1:Nwt(i)=j;j=j+1;endplot(t,g,'.');sub plot(3,2,3);S,F,T, P = sp ectrogram(y,hammi ng(600),580,600,1E3

14、);surf(T,F,10*log10( P),'edgecolor' , 'none');axis tightview(0,90);xlabel( 'Time (Seconds)');ylabel( 'Hz');sub plot(3,2,4);surf(T,F,10*log10( P),'edgecolor' , 'none');axis tightxlabel( 'Time (Seconds)');ylabel( 'Hz');zlabel( 'enargy'Y=ft(y);X=ftshift(Y);%勾x1=ifftshift(X);x1=ifft(x1);sub plot(3,2,5);plot(x1);xlabel('時間&