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

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

2、技術(shù)對(duì)信號(hào)截取，然后對(duì)截取的局部信號(hào)作Fourier變換。然后不斷地移動(dòng)窗口函數(shù)中心的位置，就可以得到信號(hào)的局部區(qū)域的瞬時(shí)頻率,因此，對(duì)于連續(xù)的信號(hào)，它的窗口傅里葉變換為:Sf(u,)窗口傅里葉逆變換為:f(t)2Sf (u, )g(t u)e td du而對(duì)于離散的信號(hào)，它的窗口傅里葉變換為: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ù)要進(jìn)行窗口傅里葉變換，首先要要選擇窗口函數(shù)，窗口函數(shù)有很多，例如高斯窗、hamming窗和Hanning窗等等。其中高斯窗函數(shù)被

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

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

5、對(duì)每一小段的做快速傅里葉變換就可以了。方法源程序: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' )得到結(jié)果:-10.5a0.5I'D43JTi-ne結(jié)果分析:我們用的是高斯窗口，得到了一個(gè)很好的結(jié)果，無(wú)論是2D圖還是3D圖，都與書上的圖十分相似。但有一個(gè)十分大的缺陷，就是無(wú)法重構(gòu)原來(lái)的信號(hào)，因?yàn)槲覀兪侵苯诱{(diào)用的spectrogram函數(shù)，并不太知道里面的具體程序是什么樣的，所以無(wú)法還原原信號(hào)，也無(wú)法計(jì)算誤差。因此我們就需要第二種方法。方法二:源程序: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ù)長(zhǎng) window lengthL=19;%

9、窗函數(shù)每次移動(dòng)的樣點(diǎn)數(shù)，重疊寬度Ts=rou nd(N-Nw)/L)+1;%計(jì)算把數(shù)據(jù)x共分成多少段n fft=2ceil(log2(Nw);%FFT的長(zhǎng)度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('時(shí)間');ylabel('頻率')sub plot(3,2,4);mesh(abs(TF);%E維繪圖axis tighttitle( 'STFT');xlabel('時(shí)間');ylabel('頻率');%短時(shí)傅里葉變換X=ff

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

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

13、lot(y);title（'信號(hào)波形圖）； 略名xlabel('時(shí)間');%x軸ylabel('振幅);%y軸hold onN=le ngth(y);Nw=600;%窗函數(shù)長(zhǎng) window lengthL=Nw/2;%窗函數(shù)每次移動(dòng)的樣點(diǎn)數(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('時(shí)間&