实验二-三

1、用三种不同的DFT程序计算x(n)=R8(n)的傅里叶变换X(ejw)，并比较三种程序计算机运行时间。(1)用forloop语句的M函数文件dft1.m，用循环变量逐点计算X(k)；(2)编写用MATLAB矩阵运算的M函数文件dft2.m，完成上述运算；(3)编写函数dft3.m，调用FFT库函数，直接计算X(k)；(4)分别利用上述三种不同方式编写的DFT程序计算序列x(n)的傅立叶变换X(ejw)，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。M函数文件如下：dft1.m：function[Am,pha]=dft1(x)N=length(x);w=exp(-j*2*pi/N);fork=1:Nsum=0;forn=1:Nsum=sum+x(n)*w^((k-1)*(n-1));endAm(k)=abs(sum);pha(k)=angle(sum);enddft2.m：function[Am,pha]=dft2(x)N=length(x);n=[0:N-1];k=[0:N-1];w=exp(-j*2*pi/N);nk=n'*k;wnk=w.^(nk);Xk=x*wnk;Am=abs(Xk);pha=angle(Xk);dft3.m：function[Am,pha]=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);源程序及运行结果：(1)x=[ones(1,8),zeros(1,248)];t=cputime;[Am1,pha1]=dft1(x);t1=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(1)subplot(2,1,1),plot(w,Am1,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha1,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians');(2)x=[ones(1,8),zeros(1,248)];t=cputime;[Am2,pha2]=dft2(x);t2=cputime-tn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)subplot(2,1,1),plot(w,Am2,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha2,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians');(3)x=[ones(1,8),zeros(1,248)];t=cputime;[Am3,pha3]=dft3(x);t3=cputime-t;n=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(3)subplot(2,1,1),plot(w,Am3,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');subplot(2,1,2),plot(w,pha3,'r');grid;title('PhasePart');xlabel('frequencyinradians');ylabel('argX[exp(jw)]/radians')从以上运行结果可以看出，调用FFT库函数直接计算X(k)速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。因此FFT算法大大提高了DFT的实用性。2、利用DFT实现两序列的卷积运算，并研究DFT点数与混叠的关系。给定x(n)=nR16(n)，h(n)=R8(n)，用FFT和IFFT分别求线性卷积和混叠结果输出，并用函数stem(n,y)画出相应图形。源程序如下：N=16;x=[0:N-1];h=ones(1,8);线性卷积：Xk1=fft(x,23);%做23点fftHk1=fft(h,23);Yk1=Xk1.*Hk1;y1=ifft(Yk1);n=0:22;figure(1)stem(n,y1)混叠：Xk2=fft(x);Hk2=fft(h,16);%做16点fftYk2=Xk2.*Hk2;y2=ifft(Yk2);n=0:15;figure(2)stem(n,y2)将两个信号的DFT相乘，再做反变换，相当于时域做循环卷积。而循环卷积是线性卷积以N为周期周期性延拓的结果，其中N为离散傅立叶变换的点数。若x1(n)是N1点序列，x2(n)是N2点序列，则其线性卷积为N1+N2-1点，做N=max(N1,N2)点的傅立叶变换，则结果中将有N2-1点是混叠的。题目中N1=16，N2=8，由程序运行结果可以看出，第二张图的前7（N2-1）个点为混叠结果，后9个点是线性卷积的结果。所以，要用DFT来做线性卷积，DFT的点数必须大于或等于线性卷积的结果。本题中取N1+N2-1=23为DFT的点数，即可得到正确的线性卷积结果。3、研究高密度频谱与高分辨率频谱。设有连续信号Xa(t)=cos(2π×6.5×10^3t)+cos(2π×7×10^3t)+cos(2π×9×10^3t)以采样频率fs=32kHz对信号Xa(t)采样，分析下列三种情况的幅频特性。(1)采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。(2)采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。(3)采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。源程序如下：fs=32000;%采样频率N=16;%采集16点n=0:N-1;%做16点DFTxa=cos(2*pi*6.5*10^3*n/fs)+cos(2*pi*7*10^3*n/fs)+cos(2*pi*9*10^3*n/fs);[Am1,pha1]=dft3(xa);n=[0:(length(xa)-1)];w=(2*pi/length(xa))*n;figure(1)plot(w,Am1,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');x=[xa(1:16),zeros(1,240)];%补零到256[Am2,pha2]=dft2(x);%做256点DFTn=[0:(length(x)-1)];w=(2*pi/length(x))*n;figure(2)plot(w,Am2,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');N0=256;%采集256点n0=0:N0-1;%做256点DFTxa0=cos(2*pi*6.5*10^3*n0/fs)+cos(2*pi*7*10^3*n0/fs)+cos(2*pi*9*10^3*n0/fs);[Am3,pha3]=dft3(xa0);n=[0:(length(xa0)-1)];w=(2*pi/length(xa0))*n;figure(3)plot(w,Am3,'b');grid;title('Magnitudepart');xlabel('frequencyinradians');ylabel('|X(exp(jw))|');运行结果：N=16点DFT幅频特性：高密度频谱：高分辨率频谱：观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。将16点DFT补零到256点，做256点DFT，由结果可以看出，频率间隔缩小（由减小为），得到了比较平滑的连续曲线，此为高密度频谱。它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。第三次做的是256点的DFT，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。