实验四-利用FFT对信号进行频谱分析1112

1、实验四 利用DFT对信号进行频谱分析一、目的要求（1） 进一步加深对线性卷积的理解和分析能力；（2） 通过编程，上机调试程序，进一步增强使用计算机解决问题的能力；（3） 掌握线性卷积与循环卷积软件实现的方法，并验证二者之间的关系；（4） 掌握利用快速卷积的计算机实现方法；（5） 熟练掌握用MATLAB实现FFT及信号的频谱分析。二、实验环境 PC机；MATLAB语言环境三、实验原理及实例分析 FFT是DFT的快速算法，凡是可以利用离散傅里叶变换DFT来进行计算的场合，都可以用利用FFT算法及数字信号处理技术加以实验。FFT在数字通信、语言信号处理、数字图像处理、匹配滤波器以及功率谱计算、系统分

2、析与仿真等诸多领域都得到了广泛的应用。（一）利用MATLAB 验证两个序列的线性卷积和圆周卷积的关系（1） 线性卷积线性时不变系统（Linear Time-Invariant System, or L. T. I系统）输入、输出间的关系为：当系统输入序列为，系统的单位脉冲响应为，输出序列为，则系统输出为：或 上式称为离散卷积或线性卷积。图4.1示出线性时不变系统的输入、输出关系。L. T. Ih(n)  L. T. I  图4.1 线性时不变系统的输入、输出关系（2） 循环卷积D F T设两个有限长序列和，均为点长 D F T 如果则 N上式称为循环卷积或圆周卷积注：为序列

3、的周期化序列；为的主值序列。上机编程计算时，可表示如下：（3） 两个有限长序列的线性卷积序列为点长，序列为点长，为这两个序列的线性卷积，则为且线性卷积的最大长，也就是说当和时。（4） 循环卷积与线性卷积的关系序列为点长，序列为点长，若序列和进行N点的循环卷积，其结果是否等于该两序列的线性卷积，完全取决于循环卷积的长度：当时循环（圆周）卷积等于线性卷积，即N当时，循环卷积等于两个序列的线性卷积加上相当于下式的时间混叠，即【实例4-1】已知两序列： 求它们的线性卷积yl(n)=h(n)*x(n)和N点的循环卷积yc=h(n)* x(n)n，并研究两者之间的关系。【MATLAB实现的参考程序】 %循

4、环卷积定义function yc=circonv(x1,x2,N)%realize circular convolution use direct method%y=circonv(x1,x2,N)%y:output sequences%x1,x2:input sequences%N:circulation lengthif length(x1)>Nerror('N must not be less than length of x1');endif length(x2)>Nerror('N must not be less than length of x

5、2');end%以上语句判断两个序列的长度是否小于Nx1=x1,zeros(1,N-length(x1);%填充序列x1(n)使其长度为N1+N2-1(序列%h(n)的长度为N1，序列x(n)的长度为N2)x2=x2,zeros(1,N-length(x2);%填充序列x2(n)使其长度为N1+N2-1n=0:1:N-1;x2=x2(mod(-n,N)+1);%生成序列x2(-n)NH=zeros(N,N);for n=1:1:NH(n,:)=cirshiftd(x2,n-1,N);%该矩阵的k行为x2(k-1-n)Nendyc=x1*H'%计算循环卷function y=ci

6、rshiftd(x,m,N)%directly realize circular shift for sequence x%y=cirshiftd(x,m,N);%x:input sequence whose length is less than N%m:how much to shift%N:circular length%y:output shifted sequenceif length(x)>N error('length of x must be less than N');endx=x,zeros(1,N-length(x);n=0:1:N-1;y=x(mo

7、d(n-m,N)+1); % 研究两者之间的关系clear all;n=0:1:11;m=0:1:5;N1=length(n);N2=length(m);xn=0.8.n;%生成x(n)hn=ones(1,N2);%生成h(n)yln=conv(xn,hn);%直接用函数conv计算线性卷积ycn=circonv(xn,hn,N1);%用函数circonv计算N1点循环卷积ny1=0:1:length(yln)-1;ny2=0:1:length(ycn)-1;subplot(2,1,1);%画图stem(ny1,yln);subplot(2,1,2);stem(ny2,ycn);axis(0,

8、16,0,4);  运行结果：图4.1 线性卷积和循环卷积的比较【任务一】请同学修改程序，使得循环卷积可以代替线性卷积。（二）利用快速卷积法计算两个序列的卷积快速卷积法的计算框图如下所示：h(n)x(n)L点FFTL点IFFTL点FFTy(n)【实例4-2】 已知序列，试利用快速卷积法计算这两个序列的卷积。MATLAB程序清单如下：Nx=15;Nh=20;n1=1:Nx-1;n2=0:Nh-1;xn=sin(0.4*n1).*(n1>=0 &n1<Nx);% 产生序列x(n)hn=0.9.n2.*(n2>=0 &n2<Nh); % 产

9、生序列h(n)L=pow2(nextpow2(Nx+Nh-1);% 序列x(n)和h(n) 卷积后得到序列y(n) 的长度Xk=fft(xn,L); % 对序列x(n)作L点DFTHk=fft(hn,L); % 对序列h(n)作L点DFTyn=ifft(Xk.*Hk) subplot(3,1,1);stem(xn,'.');title('x(n)');subplot(3,1,2);stem(hn,'.');title('h(n)');subplot(3,1,3);stem(yn,'.');title('y(

10、n)');【任务二】修改程序，比较用直接调用conv函数计算结果与FFT函数计算的结果是否一致。（三）利用FFT对信号进行频谱分析【实例4-3】已知连续信号，选择采样频率Hz，变换区间在N=16，32，64三种情况下进行谱分析，分别打印其幅频特性，并进行分析和讨论。%程序清单Fs=64;T=1/Fs;N=16;n=0:N-1;%FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样X6k16=fft(x6nT);%计算x6nT的16点DFTX6k16=fftshift(X6k16); %将零频

11、率移到频谱中心Tp=N*T;F=1/Tp;%频率分辨率Fk=-N/2:N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,1);stem(fk,abs(X6k16),'.'); %绘制16点DFT的幅频特性图box on; %坐标系右边和上边有边框，off了就没有了title('16点|DFTx(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1;%FF

12、T的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样X6k32=fft(x6nT);%计算x6nT的32点DFTX6k32=fftshift(X6k32);%将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生32点DFT对应的采样点频率（以零频率为中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %绘制32点DFT的幅频特性图title('32点|DFTx(nT)|

13、9;);xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32);N=64;n=0:N-1;%FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采X6k64=fft(x6nT); %计算x6nT的64点DFTX6k64=fftshift(X6k64); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F;%产生64点DFT对应的采样点频率（以零频率为中心）subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on %绘制64点DFT的幅频特性图title('( 64点|DFTx(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64)【任务二】若分析修改，分析运行结果。【任务三】若分析修改，分析运行结果。四、实验内容1有限长序列x(n)=2,1,0,