数字信号处理实验DFT分析连续信号频谱

1、数字信号matlab上机仿真报告题目：利用DFT分析x(t)=Acos(2pf1t)+Bcos(2pf2t)的频谱，其中f1=100Hz，f2=120Hz。 (1)A=B=1; (2)A=1,B=0.2 要求选择不同的DFT参数及窗函数(2-3类)，并对实验结果进行比较,总结出选择合适DFT参数的原则.1、a) 矩形窗截断N=30; %数据的长度L=512; %DFT的点数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;%数据的长度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;%数据的长度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; %数据的长度L=512; %DFT的点数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;%数据的长度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;%数据的长度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低幅度的分量；总结：由于矩形窗在两端变化太陡所以高频分量多，使幅度低的频率部分无法再频谱图分辨出来，所以在时域用该选取变化相对平缓的窗口函数，来避免；选择合适DFT参数的原则：1、 抽样频率/时间间隔2、 时域抽样点数或抽样时间矩形窗时取c=1，哈明窗时取c=2 3、DFT点数思考题(1) 既然可以直接计算FT，为什么利用DFT分析连续信号谱?答：根据定义是可以根据傅里叶变换的定义直接计算连续信号的傅里叶变换，但是定义区间是无限长，这在计算上是不可实施的，无论是人工计算还是通过计算机进行计算。而DF

8、T是有限长的序列的傅里叶变换，在计算机上容易实现。再者，在数字信号处理中，希望能够利用数字方法直接计算常见的四种信号的频谱函数，这是需要的时域信号为有限长，其频谱也为有限项。因此常常利用DFT对序列进行频谱分析。(2) 在利用DFT分析连续信号频谱时，会出现哪些误差？如何克服或减弱？A.混叠现象：对于带限连续信号，只要提高抽样频率使之满足时域抽样定理；对于非带限信号，可以根据实际情况对其进行低通滤波，使之成为带限信号。工程中的信号一般都不是带限信号，连续信号在抽样前通常都进过一个低通滤波器（即抗混叠滤波器）进行低通滤波，以减少混叠误差，提高频谱分析精度。B.泄漏现象：在选择矩形窗口的长度时，适

9、当增加窗的长度，可以提高频谱分辨率，但是不能减小旁瓣引起的频谱泄露，因此可以选择旁瓣幅度很小甚至为零的非矩形窗对信号进行加窗处理，就可以降低频谱泄露。C.栅栏现象：改善栅栏现象最常用的方法是在离散序列之后补零，得到一个比原有序列更长的序列，这样就可以增加频谱图中的很多细节，降低栅栏现象。(3) 在利用DFT分析连续信号频谱时，如何选择窗函数？答：在用DFT分析连续信号频谱时，选择窗函数一般首选矩形窗，因为对信号进行加窗处理的目的是去截断信号，故一般情况下选择矩形窗就可以了。但是在对频谱分析精度要求高和重要的频率点上的幅度较小的情况下，就要合理选择非矩形窗，选择旁瓣幅值小甚至为零的非矩形窗以满足

10、要求，提高频谱分析精度。(4) 若截短信号造成泄漏而导致频谱分辨率下降，可否通过在截短序列后补零得到改善？答：可以对于序列补0后可以提高频谱的分辨率。(5) 窗函数对频谱分辨率有何影响吗？怎样提高频谱分辨率？答：窗函数对频谱分析精度的影响主要是有窗函数的主瓣宽度和旁瓣幅度影响的。具体而言是：主瓣宽度越窄，精确度越高；旁瓣幅度越小，精确度越高。因此要想提高频谱分析精度，就必须选择合适的窗函数，应该使所选窗函数的长度尽量长，这样就可以减小窗函数主瓣的宽度；此外尽量选则旁瓣幅度较小甚至为零的窗函数，这样也可以减小频谱泄露，提高频谱分子精度。(6) 如何选择取样频率？答：抽样频率应满足抽样定理即：抽样

11、的最小频率要大于被测信号最高频率的2倍才能保证不发生频谱的混叠。(7) 既然频谱分辨率与信号采集时间成反比，是否意味着在实际中频谱分辨率可以很容易实现？答：采集的时间t=nT又根据当采样的时间变长时采样的点数N会增多使计算量变大，可见频率分辨率的增加时运算量会增多。(8) 讨论补零对计算结果的影响。答：在有限长的序列Xk后面补零对DTFT没有影响，但使DFT的抽样点增多得到可以得到更多的频谱细节提高频谱分辨率。课后习题：M 2-1x = 1,-3,4,2,0,-2;h = 3,0,1,-1,2,1;的循环卷积程序：x = 1,-3,4,2,0,-2;h = 3,0,1,-1,2,1;%6点循环

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('循环卷积')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循环卷积计算线性卷积的子函数%function y=lin_conv(x1,x2,L)% 利用FFT循环卷积计算线性卷积y1=fft(x1,L);y2=fft(x2,L);Yk=y1.*y2;y=ifft(Yk);end%3 重叠相加法实现长短序列的线性卷积%function y=overlap_add(x,h,M)% x为较长的输入序列，h为较短的系统函数，M为分段大小 %-前期处理-N=length(h); % 序列h(n)的长度Lx=length(x); % 序列x(n)的长度if N>

15、;M % 算法要求 N<=M M=N+1;endL=M+N-1; % 用每段循环卷积计算线性卷积所需点数T=ceil(Lx/M); % 分段数，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); % 低点数计算时的分段x(n) y_seg=lin_conv(x_seg,h,L); % 计算 i 分段和 h 的循环线性卷积 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重叠保留法 实现长短序列的线性卷积%function y=overlap_save(x,h,M)% x为较长的输入序列，h为较短的系统函数，M为分段大小 %-前期处理-N=length(h); % 序列h(n)的长度Lx=length(x); % 序列x(n)的长度if N>M % 算法要求 N<=M M=N+1;endL=M+N-1; % 用每段循环卷积计算线性卷积所需点数T=ceil(Lx/M); % 分段数，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); %