頻譜分析實驗報告

1、基于matlab的頻譜分析一、實驗目的 1頻譜分析的基本概念及應用；2常用信號的頻譜：余弦信號、方波信號、三角波信號； 3應用dft進行頻譜分析，matlab編程仿真（1）整周期采樣、非整周期采樣頻譜泄漏；（2）填零運算高密度頻譜；（3）非整周期采樣時，加窗改進；4應用頻譜分析儀測量常用信號的頻譜。二、實驗原理1、利用dft計算連續(xù)信號的頻譜分析過程采樣截短dft從圖中看出，這是一次近似過程，首先，用離散采樣信號的 dtft 來近似連續(xù)信號 的傅立葉變換 ，其次，將截短，這一過程等效于用一矩形序列與相乘，其dtft為 最后，再對截短的信號作dft2、頻譜分析的近似過程出現(xiàn)問題及解決辦法（1）混

2、疊混疊現(xiàn)象對連續(xù)信號x(t)進行數字處理前，要進行采樣采樣序列的頻譜是連續(xù)信號頻譜的周期延拓，周期為fs，如采樣率過低，不滿足采樣定理，fs2fh，則導致頻譜混迭，使一個周期內的譜對原信號譜產生失真，無法恢復原信號，進一步的數字處理失去依據。另外，干擾信號的頻率高于折疊頻率的也會產生頻譜混疊。對于干擾信號的頻率高于產生頻譜混疊的解決辦法：在采樣前加一抗混疊濾波器，濾除高于折疊頻率的頻率成分。（2） 泄漏產生原因處理實際信號序列 x（n）時，一般總要將它截斷為一有限長序列，長為n點，相當于乘以一個矩形窗 w(n)=rn(n)。矩形窗函數，其頻譜有主瓣，也有許多副瓣，窗口越大，主瓣越窄，當窗口趨于

3、無窮大時，就是一個沖擊函數。我們知道，時域的乘積對應頻域的卷積，所以，加窗后的頻譜實際是原信號頻譜與矩形窗函數頻譜的卷積，卷積的結果使頻譜延伸到了主瓣以外，且一直延伸到無窮。當窗口無窮大時，與沖擊函數的卷積才是其本身，這時無畸變，否則就有畸變。例如，信號為 ，是一單線譜，但當加窗后，線譜與抽樣函數進行卷積，原來在0處的一根譜線變成了以0為中心的，形狀為抽樣函數的譜線序列，原來在一個周期（s）內只有一個頻率上有非零值，而現(xiàn)在一個周期內幾乎所有頻率上都有非零值，即 的頻率成份從0處“泄漏”到其它頻率處去了。 考慮各采樣頻率周期間頻譜“泄漏”后的互相串漏，卷積后還有頻譜混迭現(xiàn)象產生。解決辦法：增加截

4、短的長度n，以減小泄漏。（3）柵欄效應柵欄效應的產生n點dft是在頻率區(qū)間 0，2 上對信號頻譜進行n點等間隔采樣，得到的是若干個離散的頻譜點 x（k），且它們限制在基頻的整數倍上，這就好像在柵欄的一邊通過縫隙看另一邊的景象一樣，只能在離散點處看到真實的景象，其余部分頻譜成分被遮擋， 所以稱之為柵欄效應。減小柵欄效應方法：尾部補零，使譜線變密，增加頻域采樣點數，原來漏掉的某些頻譜分量就可能被檢測出來。（4） dft的分辨率 填補零值可以改變對dtft的采樣密度，人們常常有一種誤解，認為補零可以提高dft的頻率分辨率。事實上我們通常規(guī)定dft的頻率分辨率為，這里的n是指信號x(n)的有效長度，而

5、不是補零的長度。不同長度的x(n)其dtft的結果是不同的；而相同長度的x(n)盡管補零的長度不同其dtft的結果應是相同的，他們的dft只是反映了對相同的dtft采用了不同的采樣密度。 參數選擇的一般原則： （1）若已知信號的最高頻率，為防止混疊，選定采樣頻率 ；（2）根據頻率分辯率 ，確定所需dft的長度；（3）和n確定以后，即可確定相應模擬信號的時間長度（這里t是采樣周期）。三、實驗內容及結果1、余弦波頻譜分析對余弦函數x(t)=cos(2*25t)進行頻譜分析，采樣頻率為1000hz，對其進行整周期采樣，非整周期采樣，頻譜泄露及其改善方法。其程序代碼和運行結果如下。程序代碼：fs =

6、1000;t = 0:1/fs:1.5;%xx = cos(2*pi*25*t);subplot(321),plot(t,x), axis(0 0.2 -1.2 1.2); title(x(t)=cos(2*pi*25*t);grid on;% x(n)nx=0:99;x0=cos(2*pi*1/25*nx);subplot(322),stem(nx,x0);title(x(n);grid on;% 整周期采樣n1=100;x1=x0(1:n1);x1=fft(x1,n1);magx1=abs(x1(1:n1/2+1)*2/n1;w1=0:n1/2*2*pi/n1;subplot(323);s

7、tem(w1/pi,magx1);title(整周期采樣);grid on;% 265點dft，high resolution spectrumn2=256;x=fft(x0,n2);magx=abs(x(1:n2/2+1)*2/n2;w=0:n2/2*2*pi/n2;subplot(324);stem(w/pi,magx);title(256點dft，high resolution spectrum);grid on;%加長矩形窗長度mx=0:149;xm=cos(2*pi*1/25*mx);nm=256;y=fft(xm,nm);magy=abs(y(1:nm/2+1)*2/nm;wm=0

8、:nm/2*2*pi/nm;subplot(325);stem(wm/pi,magy);title(矩形窗加長到125);grid on;%改窗n=256;wn = hamming(length(x0);hn = x0.*wn;h = fft(hn,n);magh=abs(h(1:n/2+1)*2/n;wh=0:n/2*2*pi/n;subplot(326);stem(wh/pi,magh);title(加hamming窗);grid on; 運行結果：2、方波信號頻譜對方波函數x(t)= square (2*25t)進行頻譜分析，采樣頻率為1000hz，對其進行整周期采樣，非整周期采樣，頻譜

9、泄露及其改善方法。其程序代碼和運行結果如下。程序代碼：fs = 1000;t = 0:1/fs:1.5;%xx = square(2*pi*25*t);subplot(321),plot(t,x), axis(0 0.2 -1.2 1.2); title(x(t)=square(2*pi*25*t);grid on;% x(n)nx=0:99;x0=square(2*pi*1/25*nx);subplot(322),stem(nx,x0);title(x(n);grid on;% 整周期采樣n1=100;x1=x0(1:n1);x1=fft(x1,n1);magx1=abs(x1(1:n1/2

10、+1)*2/n1;w1=0:n1/2*2*pi/n1;subplot(323);stem(w1/pi,magx1);title(整周期采樣);grid on;% 265點dft，high resolution spectrumn2=256;x=fft(x0,n2);magx=abs(x(1:n2/2+1)*2/n2;w=0:n2/2*2*pi/n2;subplot(324);stem(w/pi,magx);title(256點dft，high resolution spectrum);grid on;%加長矩形窗長度mx=0:149;xm=square(2*pi*1/25*mx);nm=256

11、;y=fft(xm,nm);magy=abs(y(1:nm/2+1)*2/nm;wm=0:nm/2*2*pi/nm;subplot(325);stem(wm/pi,magy);title(矩形窗加長到125);grid on;%改窗n=256;wn = hamming(length(x0);hn = x0.*wn;h = fft(hn,n);magh=abs(h(1:n/2+1)*2/n;wh=0:n/2*2*pi/n;subplot(326);stem(wh/pi,magh);title(加hamming窗);grid on;運行結果：3、鋸齒波信號頻譜對方波函數x(t)= sawtooth

12、 (2*25t)進行頻譜分析，采樣頻率為1000hz，對其進行整周期采樣，非整周期采樣，頻譜泄露及其改善方法。其程序代碼和運行結果如下。程序代碼：fs = 1000;t = 0:1/fs:1.5;%xx = sawtooth(2*pi*25*t);subplot(321),plot(t,x), axis(0 0.2 -1.2 1.2); title(x(t)=sawtooth(2*pi*25*t);grid on;% x(n)nx=0:99;x0=sawtooth(2*pi*1/25*nx);subplot(322),stem(nx,x0);title(x(n);grid on;% 整周期采樣

13、n1=100;x1=x0(1:n1);x1=fft(x1,n1);magx1=abs(x1(1:n1/2+1)*2/n1;w1=0:n1/2*2*pi/n1;subplot(323);stem(w1/pi,magx1);title(整周期采樣);grid on;% 265點dft，high resolution spectrumn2=256;x=fft(x0,n2);magx=abs(x(1:n2/2+1)*2/n2;w=0:n2/2*2*pi/n2;subplot(324);stem(w/pi,magx);title(256點dft，high resolution spectrum);gri

14、d on;%加長矩形窗長度mx=0:149;xm=sawtooth(2*pi*1/25*mx);nm=256;y=fft(xm,nm);magy=abs(y(1:nm/2+1)*2/nm;wm=0:nm/2*2*pi/nm;subplot(325);stem(wm/pi,magy);title(矩形窗加長到125);grid on;n=256; %改窗wn = hamming(length(x0);hn = x0.*wn;h = fft(hn,n);magh=abs(h(1:n/2+1)*2/n;wh=0:n/2*2*pi/n;subplot(326);stem(wh/pi,magh);tit

15、le(加hamming窗);grid on;運行結果：4、三角波信號頻譜對方波函數x(t)= sawtooth (2*25t , 0.5)進行頻譜分析，采樣頻率為1000hz，對其進行整周期采樣，非整周期采樣，頻譜泄露及其改善方法。其程序代碼和運行結果如下。程序代碼：fs = 1000;t = 0:1/fs:1.5;%xx = sawtooth(2*pi*25*t,0.5);subplot(321),plot(t,x), axis(0 0.2 -1.2 1.2); title(x(t)=sawtooth(2*pi*25*t,0.5);grid on;% x(n)nx=0:99;x0=sawto

16、oth(2*pi*1/25*nx,0.5);subplot(322),stem(nx,x0);title(x(n);grid on;% 整周期采樣n1=100;x1=x0(1:n1);x1=fft(x1,n1);magx1=abs(x1(1:n1/2+1)*2/n1;w1=0:n1/2*2*pi/n1;subplot(323);stem(w1/pi,magx1);title(整周期采樣);grid on;% 265點dft，high resolution spectrumn2=256;x=fft(x0,n2);magx=abs(x(1:n2/2+1)*2/n2;w=0:n2/2*2*pi/n2

17、;subplot(324);stem(w/pi,magx);title(256點dft，high resolution spectrum );grid on;%加長矩形窗長度mx=0:149;xm=sawtooth(2*pi*1/25*mx,0.5);nm=256;y=fft(xm,nm);magy=abs(y(1:nm/2+1)*2/nm;wm=0:nm/2*2*pi/nm;subplot(325);stem(wm/pi,magy);title(矩形窗加長到125);grid on;%改窗n=256;wn = hamming(length(x0);hn = x0.*wn;h = fft(hn

18、,n);magh=abs(h(1:n/2+1)*2/n;wh=0:n/2*2*pi/n;subplot(326);stem(wh/pi,magh);title(加hamming窗);grid on; 運行結果：5、頻率分辨率頻率分辨率是指所用算法能把信號中兩個靠得很近的譜峰保持分開的能力。若兩個信號的譜峰靠的很近，很難分開。為了解決這個問題，可以進行填零運算，但是填零運算產生的是高密度頻譜，不能很好的解決。另一個方法是加長原始信號長度，這樣便是產生高分辨率頻譜，便可以分開很靠近的兩個譜峰。對函數nx=0:99;x=cos(0.48*pi*nx)+cos(0.52*pi*nx)進行頻譜分析，對其進行dtft，10點dft，100點dft，填零運算，高密度頻譜和高分辨率頻譜。其程序代碼和運行結果如下。程序代碼：clear;clc;% x(n)nx=0:99;x=cos(0.48*pi*nx)+cos(0.52*pi*nx);subplot(321);stem(nx,x);title(x(n)=cos(0.48*pi*n)+cos(0.52*pi*n);grid on;% x(ejw)w=0:200*pi/200;x=dtft(x,nx,w);magx