數(shù)字信號處理實驗DFT分析連續(xù)信號頻譜

1、數(shù)字信號matlab上機仿真報告題目：利用DFT分析x(t)=Acos(2pf1t)+Bcos(2pf2t)的頻譜，其中f1=100Hz，f2=120Hz。 (1)A=B=1; (2)A=1,B=0.2 要求選擇不同的DFT參數(shù)及窗函數(shù)(2-3類)，并對實驗結果進行比較,總結出選擇合適DFT參數(shù)的原則.1、a) 矩形窗截斷N=30; %數(shù)據(jù)的長度L=512; %DFT的點數(shù)f1=100; f2=120;fs=600; %抽樣頻率T=1/fs; %抽樣間隔ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);X=fftshift(fft(x

2、,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X);ylabel('幅度譜');title('矩形窗截斷'); b) 使用hamming窗截斷N=30;%數(shù)據(jù)的長度L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);wh=(hamming(N)'x=x.*wh;X=fftshift(fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X);

3、ylabel('幅度');xlabel('頻率');title('hamming窗口截斷')c) 使用blackman截斷N=30;%數(shù)據(jù)的長度L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);wh=(blackman(N)'x=x.*wh;X=fftshift(fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X);ylabel('幅度');xl

4、abel('頻率');title('blackman窗口截斷')2、a) 矩形窗截斷：N=30; %數(shù)據(jù)的長度L=512; %DFT的點數(shù)f1=100; f2=120;fs=600; %抽樣頻率T=1/fs; %抽樣間隔ws=2*pi*fs;t=(0:N-1)*T;f=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);F=fftshift(fft(f,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);hd=plot(w,abs(F);ylabel('幅度譜');title('使用矩形窗截斷');當

5、采樣點增加到300時對應的頻譜圖：旁瓣高頻十分多無法找的0.2*cos(2*pi*f2*t)的幅度低的無法分辨；b) Hamming窗截斷N=30;%數(shù)據(jù)的長度L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);wh=(hamming(N)'x=x.*wh;X=fftshift(fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X);ylabel('幅度');xlabel('頻率&

6、#39;);title('使用hamming截斷')c) 使用blackman截斷N=30;%數(shù)據(jù)的長度L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);wh=(blackman(N)'x=x.*wh;X=fftshift(fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X);ylabel('幅度');xlabel('頻率');title('使用

7、blackman截斷')使用hamming和blackman截斷可以清楚的分辨120hz低幅度的分量；總結：由于矩形窗在兩端變化太陡所以高頻分量多，使幅度低的頻率部分無法再頻譜圖分辨出來，所以在時域用該選取變化相對平緩的窗口函數(shù)，來避免；選擇合適DFT參數(shù)的原則：1、 抽樣頻率/時間間隔2、 時域抽樣點數(shù)或抽樣時間矩形窗時取c=1，哈明窗時取c=2 3、DFT點數(shù)思考題(1) 既然可以直接計算FT，為什么利用DFT分析連續(xù)信號譜?答：根據(jù)定義是可以根據(jù)傅里葉變換的定義直接計算連續(xù)信號的傅里葉變換，但是定義區(qū)間是無限長，這在計算上是不可實施的，無論是人工計算還是通過計算機進行計算。而DF

8、T是有限長的序列的傅里葉變換，在計算機上容易實現(xiàn)。再者，在數(shù)字信號處理中，希望能夠利用數(shù)字方法直接計算常見的四種信號的頻譜函數(shù)，這是需要的時域信號為有限長，其頻譜也為有限項。因此常常利用DFT對序列進行頻譜分析。(2) 在利用DFT分析連續(xù)信號頻譜時，會出現(xiàn)哪些誤差？如何克服或減弱？A.混疊現(xiàn)象：對于帶限連續(xù)信號，只要提高抽樣頻率使之滿足時域抽樣定理；對于非帶限信號，可以根據(jù)實際情況對其進行低通濾波，使之成為帶限信號。工程中的信號一般都不是帶限信號，連續(xù)信號在抽樣前通常都進過一個低通濾波器（即抗混疊濾波器）進行低通濾波，以減少混疊誤差，提高頻譜分析精度。B.泄漏現(xiàn)象：在選擇矩形窗口的長度時，適

9、當增加窗的長度，可以提高頻譜分辨率，但是不能減小旁瓣引起的頻譜泄露，因此可以選擇旁瓣幅度很小甚至為零的非矩形窗對信號進行加窗處理，就可以降低頻譜泄露。C.柵欄現(xiàn)象：改善柵欄現(xiàn)象最常用的方法是在離散序列之后補零，得到一個比原有序列更長的序列，這樣就可以增加頻譜圖中的很多細節(jié)，降低柵欄現(xiàn)象。(3) 在利用DFT分析連續(xù)信號頻譜時，如何選擇窗函數(shù)？答：在用DFT分析連續(xù)信號頻譜時，選擇窗函數(shù)一般首選矩形窗，因為對信號進行加窗處理的目的是去截斷信號，故一般情況下選擇矩形窗就可以了。但是在對頻譜分析精度要求高和重要的頻率點上的幅度較小的情況下，就要合理選擇非矩形窗，選擇旁瓣幅值小甚至為零的非矩形窗以滿足

10、要求，提高頻譜分析精度。(4) 若截短信號造成泄漏而導致頻譜分辨率下降，可否通過在截短序列后補零得到改善？答：可以對于序列補0后可以提高頻譜的分辨率。(5) 窗函數(shù)對頻譜分辨率有何影響嗎？怎樣提高頻譜分辨率？答：窗函數(shù)對頻譜分析精度的影響主要是有窗函數(shù)的主瓣寬度和旁瓣幅度影響的。具體而言是：主瓣寬度越窄，精確度越高；旁瓣幅度越小，精確度越高。因此要想提高頻譜分析精度，就必須選擇合適的窗函數(shù)，應該使所選窗函數(shù)的長度盡量長，這樣就可以減小窗函數(shù)主瓣的寬度；此外盡量選則旁瓣幅度較小甚至為零的窗函數(shù)，這樣也可以減小頻譜泄露，提高頻譜分子精度。(6) 如何選擇取樣頻率？答：抽樣頻率應滿足抽樣定理即：抽樣

11、的最小頻率要大于被測信號最高頻率的2倍才能保證不發(fā)生頻譜的混疊。(7) 既然頻譜分辨率與信號采集時間成反比，是否意味著在實際中頻譜分辨率可以很容易實現(xiàn)？答：采集的時間t=nT又根據(jù)當采樣的時間變長時采樣的點數(shù)N會增多使計算量變大，可見頻率分辨率的增加時運算量會增多。(8) 討論補零對計算結果的影響。答：在有限長的序列Xk后面補零對DTFT沒有影響，但使DFT的抽樣點增多得到可以得到更多的頻譜細節(jié)提高頻譜分辨率。課后習題：M 2-1x = 1,-3,4,2,0,-2;h = 3,0,1,-1,2,1;的循環(huán)卷積程序：x = 1,-3,4,2,0,-2;h = 3,0,1,-1,2,1;%6點循環(huán)

12、卷積Yc1N1=6;Xm=fft(x,N1);Hm=fft(h,N1);Y=ifft(Xm.*Hm)運算結果：Y =6.0000 -3.0000 17.0000 -2.0000 7.0000 -13.0000與手算結果相同。x = 1,-3,4,2,0,-2;h = 3,0,1,-1,2,1;N1=6;Xm=fft(x,N1);Hm=fft(h,N1);k=1:6;Y=ifft(Xm.*Hm)stem(k,Y)axis(0 10,-20 20)xlabel('k')ylabel('Y')title('循環(huán)卷積')M2-2：n=1018;k1=0:

13、1023;k2=0:5;x = cos(pi*k1/2)h =power(3,k2)%3的k2次冪y=conv(x,h)l=length(y)stem(0:l-1,y)title('直接卷積')對于長序列使用分段卷積的兩種方法：1序列輸入部分%k1=0:29;k2=0:5;x = cos(pi*k1/2)h =power(3,k2)%3的k2次冪a=length(x)b=length(h)y=overlap_add(x,h,6)figure stem(1:a+b-1,y)title('使用重疊相加法30點')figurey=overlap_save(x,h,6)

14、stem(1:a+b-1,y)title('使用重疊保留法30點')%2利用FFT循環(huán)卷積計算線性卷積的子函數(shù)%function y=lin_conv(x1,x2,L)% 利用FFT循環(huán)卷積計算線性卷積y1=fft(x1,L);y2=fft(x2,L);Yk=y1.*y2;y=ifft(Yk);end%3 重疊相加法實現(xiàn)長短序列的線性卷積%function y=overlap_add(x,h,M)% x為較長的輸入序列，h為較短的系統(tǒng)函數(shù)，M為分段大小 %-前期處理-N=length(h); % 序列h(n)的長度Lx=length(x); % 序列x(n)的長度if N>

15、;M % 算法要求 N<=M M=N+1;endL=M+N-1; % 用每段循環(huán)卷積計算線性卷積所需點數(shù)T=ceil(Lx/M); % 分段數(shù)，ceil向上取整t=zeros(1,N-1); % 緩存序列初始化x=x,zeros(1,(T+1)*M-Lx); % 最后一個不足M的分段補零y=zeros(1,(T+1)*M); % 生成輸出序列y(n) %-核心算法- for i=0:T xi=i*M+1; x_seg=x(xi:xi+M-1); % 低點數(shù)計算時的分段x(n) y_seg=lin_conv(x_seg,h,L); % 計算 i 分段和 h 的循環(huán)線性卷積 y_seg(1:

16、N-1)=y_seg(1:N-1)+t(1:N-1); % 與前一段卷積的后N-1位重疊相加 t(1:N-1)=y_seg(M+1:L); % 緩存序列更新 y(xi:xi+M-1)=y_seg(1:M); % 每卷積一段輸出M個點 endy=y(1:Lx+N-1) ; % 最終輸出序列end%4重疊保留法 實現(xiàn)長短序列的線性卷積%function y=overlap_save(x,h,M)% x為較長的輸入序列，h為較短的系統(tǒng)函數(shù)，M為分段大小 %-前期處理-N=length(h); % 序列h(n)的長度Lx=length(x); % 序列x(n)的長度if N>M % 算法要求 N<=M M=N+1;endL=M+N-1; % 用每段循環(huán)卷積計算線性卷積所需點數(shù)T=ceil(Lx/M); % 分段數(shù)，ceil向上取整t=zeros(1,N-1); % 緩存序列初始化x=x,zeros(1,(T+1)*M-Lx); % 最后一個不足M的分段補零y=zeros(1,(T+1)*M); % 生成輸出序列y(n) %-核心算法- for i=0:T; xi=i*M+1; x_seg=t,x(xi:xi+M-1); %