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

1、數(shù)字信號matlab上機(jī)仿真報(bào)告題目：利用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類)，并對實(shí)驗(yàn)結(jié)果進(jìn)行比較,總結(jié)出選擇合適DFT參數(shù)的原則.1、a) 矩形窗截?cái)郚=30; %數(shù)據(jù)的長度L=512; %DFT的點(diǎn)數(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('矩形窗截?cái)?#39;); b) 使用hamming窗截?cái)郚=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ái)?#39;)c) 使用blackman截?cái)郚=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窗口截?cái)?#39;)2、a) 矩形窗截?cái)啵篘=30; %數(shù)據(jù)的長度L=512; %DFT的點(diǎn)數(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('使用矩形窗截?cái)?#39;);當(dāng)

5、采樣點(diǎn)增加到300時對應(yīng)的頻譜圖：旁瓣高頻十分多無法找的0.2*cos(2*pi*f2*t)的幅度低的無法分辨；b) Hamming窗截?cái)郚=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ái)?#39;)c) 使用blackman截?cái)郚=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截?cái)?#39;)使用hamming和blackman截?cái)嗫梢郧宄姆直?20hz低幅度的分量；總結(jié)：由于矩形窗在兩端變化太陡所以高頻分量多，使幅度低的頻率部分無法再頻譜圖分辨出來，所以在時域用該選取變化相對平緩的窗口函數(shù)，來避免；選擇合適DFT參數(shù)的原則：1、 抽樣頻率/時間間隔2、 時域抽樣點(diǎn)數(shù)或抽樣時間矩形窗時取c=1，哈明窗時取c=2 3、DFT點(diǎn)數(shù)思考題(1) 既然可以直接計(jì)算FT，為什么利用DFT分析連續(xù)信號譜?答：根據(jù)定義是可以根據(jù)傅里葉變換的定義直接計(jì)算連續(xù)信號的傅里葉變換，但是定義區(qū)間是無限長，這在計(jì)算上是不可實(shí)施的，無論是人工計(jì)算還是通過計(jì)算機(jī)進(jìn)行計(jì)算。而DF

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

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

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

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

12、卷積Yc1N1=6;Xm=fft(x,N1);Hm=fft(h,N1);Y=ifft(Xm.*Hm)運(yùn)算結(jié)果：Y =6.0000 -3.0000 17.0000 -2.0000 7.0000 -13.0000與手算結(jié)果相同。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點(diǎn)')figurey=overlap_save(x,h,6)

14、stem(1:a+b-1,y)title('使用重疊保留法30點(diǎn)')%2利用FFT循環(huán)卷積計(jì)算線性卷積的子函數(shù)%function y=lin_conv(x1,x2,L)% 利用FFT循環(huán)卷積計(jì)算線性卷積y1=fft(x1,L);y2=fft(x2,L);Yk=y1.*y2;y=ifft(Yk);end%3 重疊相加法實(shí)現(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)卷積計(jì)算線性卷積所需點(diǎn)數(shù)T=ceil(Lx/M); % 分段數(shù)，ceil向上取整t=zeros(1,N-1); % 緩存序列初始化x=x,zeros(1,(T+1)*M-Lx); % 最后一個不足M的分段補(bǔ)零y=zeros(1,(T+1)*M); % 生成輸出序列y(n) %-核心算法- for i=0:T xi=i*M+1; x_seg=x(xi:xi+M-1); % 低點(diǎn)數(shù)計(jì)算時的分段x(n) y_seg=lin_conv(x_seg,h,L); % 計(jì)算 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個點(diǎn) endy=y(1:Lx+N-1) ; % 最終輸出序列end%4重疊保留法 實(shí)現(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)卷積計(jì)算線性卷積所需點(diǎn)數(shù)T=ceil(Lx/M); % 分段數(shù)，ceil向上取整t=zeros(1,N-1); % 緩存序列初始化x=x,zeros(1,(T+1)*M-Lx); % 最后一個不足M的分段補(bǔ)零y=zeros(1,(T+1)*M); % 生成輸出序列y(n) %-核心算法- for i=0:T; xi=i*M+1; x_seg=t,x(xi:xi+M-1); %