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

实验四 利用DFT 对信号进行频谱分析一、目的要求（1） 进一步加深对线性卷积的理解和分析能力；（2） 通过编程，上机调试程序，进一步增强使用计算机解决问题的能力； （3） 掌握线性卷积与循环卷积软件实现的方法，并验证二者之间的关系； （4） 掌握利用快速卷积的计算机实现方法；（5） 熟练掌握用MATLAB 实现FFT 及信号的频谱分析。

二、实验环境PC 机；MATLAB 语言环境三、实验原理及实例分析FFT 是DFT 的快速算法，凡是可以利用离散傅里叶变换DFT 来进行计算的场合，都可以用利用FFT 算法及数字信号处理技术加以实验。

FFT 在数字通信、语言信号处理、数字图像处理、匹配滤波器以及功率谱计算、系统分析与仿真等诸多领域都得到了广泛的应用。

（一）利用MATLAB 验证两个序列的线性卷积和圆周卷积的关系（1） 线性卷积线性时不变系统（Linear Time-Invariant System, or L. T. I 系统）输入、输出间的关系为：当系统输入序列为)(n x ，系统的单位脉冲响应为)(n h ，输出序列为)(n y ，则系统输出为：∑∞-∞==-=m n h n x m n h m x n y )(*)()()()(或∑+∞-∞==-=m n x n h m n x m h n y )(*)()()()(上式称为离散卷积或线性卷积。

图4.1示出线性时不变系统的输入、输出关系。

)(n δ→ L. T. I —→)(n h —→ —→图4.1 线性时不变系统的输入、输出关系（2） 循环卷积设两个有限长序列)(1n x 和)(2n x ，均为N 点长)(1n x )(1k X )(2n x )(2k X如果)()()(213k X k X k X ⋅=则 )()(~)(~)(10213n R m n x m x n x N N m ⎥⎦⎤⎢⎣⎡-=∑-=[]∑---=1021)()(N m N m n x m x)(1n x =N 10)(2-≤≤N n n x上式称为循环卷积或圆周卷积注：)(~1n x 为)(1n x 序列的周期化序列；)()(~1n R n x N 为)(~1n x 的主值序列。

上机编程计算时，)(3n x 可表示如下：∑∑-+==-++-=1121213)()()()()(N n m nm m n N xm x m n x m x n x（3） 两个有限长序列的线性卷积序列)(1n x 为L 点长，序列)(2n x 为P 点长，)(3n x 为这两个序列的线性卷积，则)(3n x 为∑+∞-∞=-=m m n xm x n x )()()(213且线性卷积)(3n x 的最大长1-+P L ，也就是说当1-≤n 和1-+≥P L n 时0)(3=n x 。

（4） 循环卷积与线性卷积的关系序列)(1n x 为L 点长，序列)(2n x 为P 点长，若序列)(1n x 和)(2n x 进行N 点的循环卷)(n xL. T. Ih(n)∑+∞-∞=-=m m n h m x n y )()()(D F TD F T积，其结果是否等于该两序列的线性卷积，完全取决于循环卷积的长度：当1-+≥P L N 时循环（圆周）卷积等于线性卷积，即)(1n x N )(*)()(212n x n x n x =当1-+<P L N 时，循环卷积等于两个序列的线性卷积加上相当于下式的时间混叠，即⎪⎩⎪⎨⎧-≤≤+=∑+∞-∞=nN n rN n x n x r N 其它010)()(33【实例4-1】已知两序列：⎪⎩⎪⎨⎧-≤≤+=∑+∞-∞=nN n rN n x n x r N 其它010)()(33求它们的线性卷积yl(n)=h(n)*x(n)和N 点的循环卷积yc=[h(n)* x(n)]n ，并研究两者之间的关系。

【MATLAB 实现的参考程序】 %循环卷积定义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 length if length(x1)>Nerror('N must not be less than length of x1'); endif length(x2)>Nerror('N must not be less than length of x2'); end%以上语句判断两个序列的长度是否小于N x1=[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=cirshiftd(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)>Nerror('length of x must be less than N');endx=[x,zeros(1,N-length(x))];n=[0:1:N-1];y=x(mod(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,16,0,4]);运行结果：01234线性卷积01234循环卷积图4.1 线性卷积和循环卷积的比较【任务一】请同学修改程序，使得循环卷积可以代替线性卷积。

（二）利用快速卷积法计算两个序列的卷积快速卷积法的计算框图如下所示：【实例】已知序列)(9.0)(),()4.0sin()(2015n R n h n R n n x n==，试利用快速卷积法计算这两个序列的卷积)(*)()(n h n x n y =。

L程序清单如下：Nx=15;Nh=20;h (n )x (n )L 点FFTL 点IFFTL 点FFTy (n )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); % 产生序列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(n)');【任务二】修改程序，比较用直接调用conv函数计算结果与FFT函数计算的结果是否一致。

（三）利用FFT 对信号进行频谱分析【实例4-3】已知连续信号，)20cos()16cos()8cos()(t t t t x πππ++=选择采样频率64=s f 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); %将零频率移到频谱中心 Tp=N*T;F=1/Tp;%频率分辨率F k=-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点|DFT[x(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;%FFT 的变换区间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点DFT X6k32=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点|DFT[x(nT)]|');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点DFT X6k64=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点|DFT[x(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])【任务二】若分析修改128,70,40=s f ，分析运行结果。