課程設(shè)計(jì)（論文）_基于MATLAB的快速傅立葉分析程序設(shè)計(jì)

1、基于MATLAB的快速傅立葉分析程序設(shè)計(jì)機(jī)械研07，徐微，07202110051.信號(hào)數(shù)據(jù)對(duì)一個(gè)人為產(chǎn)生的信號(hào)進(jìn)行采用FFT變換方法進(jìn)行功率譜分析。信號(hào)x(n)=120.0*COS(2*3.14*SF*n/FS)式中： n=0,1,2 N-1 SF-信號(hào)頻率 FS-采樣頻率這里，定義參數(shù)如下：fs=200；%設(shè)定采樣頻率N=512；sf=10；%設(shè)定余弦信號(hào)頻率采樣點(diǎn)=1024；2信號(hào)的時(shí)域波形和頻域波形圖1 余弦函數(shù)的時(shí)域波形圖圖2 余弦函數(shù)的幅頻譜波形圖3 余弦函數(shù)的功率譜密度波形上面三圖的程序如下所示：fs=200;%設(shè)定采樣頻率N=512;n=0:N-1;t=n/fs;sf=10;%設(shè)

2、定正弦信號(hào)頻率%生成信號(hào)x=120.0*cos(2*3.14*sf*t);figure;plot(t,x);%作余弦信號(hào)的時(shí)域波形xlabel(t);ylabel(y);title(x=120.0*cos(2*3.14*sf*t)時(shí)域波形);grid;%進(jìn)行FFT變換并做頻譜圖y=fft(x,N);%進(jìn)行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length(y);%進(jìn)行對(duì)應(yīng)的頻率轉(zhuǎn)換figure;plot(f,mag);%做頻譜圖xlabel(頻率(Hz);ylabel(幅值);title(x=120.0*cos(2*3.14*sf*t)幅頻譜圖N=5

3、12);grid;Py =2*(y.*conj(y)/N; %計(jì)算功率譜密度Pyfigure;plot(f,Py);xlabel(頻率(Hz);ylabel(功率譜密度);title(x=120.0*cos(2*3.14*sf*t)功率譜密度);grid;3信號(hào)經(jīng)過(guò)均值化處理或不經(jīng)過(guò)均值化處理的結(jié)果比擬圖4 余弦函數(shù)均值化處理后的時(shí)頻域波形結(jié)論：經(jīng)過(guò)均值化處理的波形發(fā)生上下移動(dòng)，但是頻譜不變。4采用不同窗函數(shù)時(shí)的譜結(jié)果矩形窗函數(shù), 漢寧窗函數(shù)，漢明窗等等1采用矩形窗函數(shù)時(shí)信號(hào)的頻譜分析結(jié)果圖圖5 矩形窗函數(shù)處理信號(hào)時(shí)域波形圖6 矩形窗函數(shù)處理信號(hào)頻域波形2采用漢明窗函數(shù)時(shí)信號(hào)的頻譜分析結(jié)果圖7

4、 漢明窗函數(shù)處理信號(hào)時(shí)域波形圖8 漢明窗函數(shù)處理信號(hào)頻域波形3采用漢寧窗函數(shù)時(shí)信號(hào)的頻譜分析結(jié)果圖9 漢寧窗函數(shù)處理信號(hào)時(shí)域波形圖10 漢寧窗函數(shù)處理信號(hào)頻域波形4采用三角窗函數(shù)時(shí)信號(hào)的頻譜分析結(jié)果圖11 三角窗函數(shù)處理信號(hào)時(shí)域波形圖12 三角窗函數(shù)處理信號(hào)頻域波形5采用Blackman窗函數(shù)時(shí)信號(hào)的頻譜分析結(jié)果圖13 Blackman窗函數(shù)處理信號(hào)時(shí)域波形圖14 Blackman窗函數(shù)處理信號(hào)頻域波形在前面程序的根底上，繼續(xù)輸入加窗處理程序，可以得到上面所示的結(jié)果，其加窗處理程序如下：w_han=(hanning(N);y1=x.*w_han;figure;plot(t,y1);y2=mag

5、.*w_han;figure;plot(f,y2);w_box=(boxcar(N);y3=x.*w_box;figure;plot(t,y3);y4=mag.*w_box;figure;plot(f,y4)w_ham=(hamming(N);y5=x.*w_ham;figure;plot(t,y5);y6=mag.*w_ham;figure;plot(f,y6)w_tri=(triang(N);y7=x.*w_tri;figure;plot(t,y7);y8=mag.*w_tri;figure;plot(f,y8)w_black=(blackman(N);y9=x.*w_black;figu

6、re;plot(t,y9);y10=mag.*w_black;figure;plot(f,y10)5典型函數(shù)的頻譜矩形窗函數(shù), 漢寧窗函數(shù)，三角窗函數(shù)，切比雪夫窗設(shè)計(jì)方法：主要應(yīng)用了MALTAB中的交互式圖形用戶界面以及直接編程來(lái)做信號(hào)的處理過(guò)程，其設(shè)計(jì)過(guò)程如下：1圖形用戶界面的啟動(dòng)：在MATLAB COMMAND窗口下，鍵入sptool，會(huì)彈出一個(gè)SPTool窗口。2在進(jìn)行處理之前，我們需要建立一個(gè)所要處理的信號(hào)公式，即信號(hào)x(n)=120.0*cos(2*3.14*SF*n/FS)，MATLAB提供了編程的代碼，其代碼如下：Fs=200;%設(shè)定采樣頻率 N=512; n=0:N-1; t=

7、n/Fs;sf=10;%設(shè)定余弦信號(hào)頻率x=120.0*cos(2*3.14*sf*t); %生成信號(hào)plot(t,x)gridsave hdata.mat x Fs %把信號(hào)保存到了工作空間中，以備調(diào)用這樣程序運(yùn)行結(jié)果會(huì)生成信號(hào)數(shù)據(jù)文件hdata.mat，存放信號(hào)x和采樣頻率的數(shù)據(jù)。3從SPTool窗口中的File菜單中選擇Import命令，彈出Import to SPTool窗。窗口如下所示：圖15 信號(hào)處理Import to SPTool窗口設(shè)置好后，點(diǎn)擊OK。4這樣就可以進(jìn)行信號(hào)的時(shí)域分析了，只需要點(diǎn)擊下列圖中Signals中所導(dǎo)入的信號(hào)就可以了。圖16 MATLAB信號(hào)處理SPToo

8、l窗口5如上圖所示的Spectra為信號(hào)的頻譜分析工具欄，選好所導(dǎo)入的信號(hào)，點(diǎn)擊Crate按鈕，進(jìn)入頻譜分析的窗口。窗口圖片如下所示：圖17 信號(hào)的頻譜分析工具欄Spectra窗口通過(guò)上面圖形所示，我們可以從左邊的工具欄上選擇各種窗函數(shù)，并且可以隨意定義Nfft和Nwind的大小，方法比擬簡(jiǎn)單直觀。同樣，我們還可以用編程的方法實(shí)現(xiàn)上面要求，其各種窗函數(shù)及程序代碼如下所示：1) 矩形窗函數(shù)圖18 矩形窗函數(shù)頻譜圖2漢寧窗函數(shù)圖19 漢寧窗函數(shù)頻譜圖3三角窗函數(shù)圖20 三角窗函數(shù)頻譜圖4切比雪夫窗函數(shù)圖21 切比雪夫窗函數(shù)頻譜圖典型窗函數(shù)的頻譜程序如下所示：生成一個(gè)長(zhǎng)度為50的矩形窗，并觀察其頻率

9、特性n=50;window=boxcar(n);h,w=freqz(window,1);subplot(2,1,1)stem(window);subplot(2,1,2);plot(w/pi,20*log(abs(h)/abs(h(1);生成一個(gè)長(zhǎng)度為60的漢寧窗，并觀察其頻率特性n=60;window=hanning(n);h,w=freqz(window,1);subplot(1,2,1)stem(window);subplot(1,2,2); plot(w/pi,20*log(abs(h)/abs(h(1);生成一個(gè)長(zhǎng)度為40的三角窗，并觀察其頻率特性n=40;window=triang

10、(n);h,w=freqz(window,1);subplot(1,2,1)stem(window);subplot(1,2,2);plot(w/pi,20*log(abs(h)/abs(h(1);生成一個(gè)長(zhǎng)度為40的切比雪夫窗，并觀察其頻率特性n=50;r=50;window=chebwin(n,r);h,w=freqz(window,1);subplot(1,2,1)stem(window);subplot(1,2,2) plot(w/pi,20*log(abs(h)/abs(h(1);6整周期和非整周期比擬以高度為1的矩形窗w(t)在一個(gè)完整周期內(nèi)截取采樣序列，在時(shí)域以w(t)截取采樣序列，在頻域內(nèi)等效于W(f)和X(f)的卷積。在fm/T0,處正好是X(f) *W(f)=W(f -f0)的零點(diǎn)處，故在頻譜上僅在f0處有單一頻線，無(wú)泄漏現(xiàn)象。因此，只要矩形窗寬包含序列的假設(shè)干完整周期，就都能在f=m/T0處使采樣序列的加窗譜不歪曲其真實(shí)譜。當(dāng)用矩形窗函數(shù)w(t)對(duì)周期信號(hào)進(jìn)行非整周期截取，如上圖所示產(chǎn)生明顯的“泄漏現(xiàn)象