离散傅里叶变换和快速傅里叶变换

.专业整理.实验报告课程名称： 信号分析与处理 指导老师： 成绩：__________________实验名称：离散傅里叶变换和快速傅里叶变换 实验类型： 基础实验 同组学生姓名 ：第二次实验 离散傅里叶变换和快速傅里叶变换一、实验目的1.1掌握离散傅里叶变换 （DFT）的原理和实现 ；1.2掌握快速傅里叶变换 （FFT）的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法 。1.3会用Matlab 软件进行以上练习 。装二、实验原理订2.1关于DFT的相关知识线序列x(n)的离散事件傅里叶变换 （DTFT）表示为X(ej)x(n)ejn，n如果()为因果有限长序列，n=0,1,...,N-1，则()的DTFT表示为xnxnN1X(ej)x(n)ejn，n0x(n)的离散傅里叶变换（DFT）表达式为N1j2nkX(k)x(n)eN(k0,1,...,N1)，n0.学习帮手..专业整理.序列的N点DFT是序列DTFT在频率区间[0,2π]上的N点灯间隔采样，采样间隔为2π/N。通过DFT，可以完成由一组有限个信号采样值 x(n)直接计算得到一组有限个频谱采样值 X(k)。X(k)的幅度谱为X(k)22(k)，其中下标R和I分别表示取实部、虚部的运算。Xk的相位谱为XR(k)XI()(k)arctanXI(k)。XR(k)离散傅里叶反变换 （IDFT）定义为1N1j2nkx(n)NNn0X(k)e(n0,1,...,N1)。2.2关于FFT的相关知识j2 n快速傅里叶变换 （FFT）是DFT的快速算法，并不是一个新的映射 。FFT利用了e N 函数的周期性和对称性以及一些特殊值来减少 DFT的运算量，可使DFT的运算量下降几个数量级 ，从而使数字信号处理的速度大大提高 。若信号是连续信号 ，用FFT进行谱分析时 ，首先必须对信号进行采样 ，使之变成离散信号 ，然后就可以用 FFT来对连续信号进行谱分析 。为了满足采样定理 ，一般在采样之前要设置一个抗混叠低通滤波器，且抗混叠滤波器的截止频率不得高于与采样频率的一半 。比较DFT和IDFT的定义，两者的区别仅在于指数因子的指数部分的符号差异和幅度尺度变换 ，因此可用FFT算法来计算 IDFT。三、实验内容与相关分析 （共6道）说明：为了便于老师查看 ，现将各题的内容写在这里 ——题目按照 3.1、3.2、...、3.6排列。每道题包含如下内容 ：题干、解答（思路、M文件源代码、命令.学习帮手..专业整理.窗口中的运行及其结果 ）、分析。其中“命令窗口中的运行及其结果 ”按照小题顺序排列 ，各小题包含命令与结果（图形或者序列）。3.1 求有限长离散时间信号 x(n)的离散时间傅里叶变换（DTFT）X(ejΩ)并绘图。．．12n20n10。（1）已知x(n)；（2）已知x(n)2n0其他【解答】思路：这是DTFT的变换，按照定义编写 DTFT的M文件即可。考虑到自变量Ω是连续的，为了方便计算机计算，计算时只取三个周期 [-2π,4π]中均匀的 1000个点用于绘图 。理论计算的各序列 DTFT表达式，请见本题的分析 。M文件源代码（我的Matlab 源文件不支持中文注释 ，抱歉）：function DTFT(n1,n2,x)%ThisisaDTFTfunctionformyexperimentofSignalProcessing&Analysis.w=0:2*pi/1000:2*pi; %Definethebracketofomegaforplotting.X=zeros(size(w));%DefinetheinitialvaluesofX.fori=n1:n2X=X+x(i-n1+1)*exp((-1)*j*w*i); %ItisthedefinitionofDTFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);plot(w,Amp,'r');xlabel('\Omega' );ylabel('Amplification' );hold on;%Plotamplificationontheleft..学习帮手..专业整理.subplot(1,2,2);plot(w,Phs,'b');xlabel('\Omega' );ylabel('PhaseAngle(radian)' );holdoff;%Plotphaseangleontheright.end命令窗口中的运行及其结果 （理论计算的各序列 DTFT表达式，请见本题的分析 ）：第（1）小题n=(-2:2);x=1.^n;DTFT(-2,2,x);544.53423.5)na1n3doat(aeci2.5l0fgnlpAmeA2sa-1hP1.5-210.5-300510-40510-5-5图在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）第（2）小题n=(0:10);x=2.^n;DTFT(0,10,x);.学习帮手..专业整理.220042000318002)1600n1andoat(aegf14000nlpAmeAs1200a-1hP1000-2800-36000510-40510-5-5图在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）【分析】对于第（1）小题，由于序列 x(n)只在有限区间（-2，-1，-，1，2）上为1，所以是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号 。事实上，我们可计算出它的表达式 ：22j1e5j1ejnjne

5jX()x(n)eejX()nn21e1e

j

，可见幅度频谱拥有主极大和次极大，两个主极大间有 |5-1|=4个极小，即有3个次级大。而对于它的相位频谱 ，则是周期性地在π、0、π之间震荡。对于第（2）小题，由于是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。而它的表10n1111j11达式：X()x(n)ejn2ej12eX()2，因此主极大之间只有jnn012e12ej|0-1|=1个极小，不存在次级大。而对于它的相位频谱，则是在一个长为2π的周期内有|11-1|=10次振荡。而由DTFT的定义可知，频谱都是以 2π为周期向两边无限延伸的 。由于DTFT是连续谱，对于计算机处理来说特别困难 ，因此我们才需要离散信号的频谱也离散 ，由此构造出 DFT（以及为加速计算 DFT的.学习帮手..专业整理.FFT）。3.2已知有限长序列 x(n)={8,7,9,5,1,7,9,5},试分别采用 DFT和FFT求其离散傅里叶变换 X(k)的幅度、相位图。【解答】思路：按照定义编写 M文件即可。M文件源代码：i)DFT函数：function DFT(N,x)%ThisisaDFTfunctionformyexperimentofSignalProcessing&Analysis.k=(0:N-1);%DefinevariablekforDFT.X=zeros(size(k));%DefinetheinitialvalvesofX.fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i); %ItisthedefinitionofDFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);stem(k,Amp,'.',’MarkerSize’,18);xlabel('k');ylabel('Amplification' );hold on;%Plotamplificationontheleft.subplot(1,2,2);stem(k,Phs,'*');xlabel('k');ylabel('PhaseAngle(radian)' );holdoff;.学习帮手..专业整理.%Plotphaseangleontheright.endii)基2-FFT函数function myFFT(N,x)%Thisisabase-2FFTfunction.lov=(0:N-1);j1=0;fori=1:N %indexedaddressingifi<j1+1temp=x( j1+1);x(j1+1)=x(i);x(i)=temp;endk=N/2;whilek<=j1j1=j1-k;k=k/2;endj1=j1+k;enddigit=0;k=N;.学习帮手..专业整理.whilek>1digit=digit+1;k=k/2;endn=N/2;%Nowwestartthe"butterfly-shaped"process.formu=1:digitdif=2^(mu-1); %Differncebetweentheindexesofthetargetvariables.idx=1;fori=1:nidx1=idx;idx2=1;forj1=1:N/(2*n)r=(idx2-1)*2^(digit-mu);wn=exp(j*(-2)*pi*r/N); %Itisthe"circulatingcoefficients".temp=x(idx);x(idx)=temp+x(idx+dif)*wn;x(idx+dif)=temp-x(idx+dif)*wn;idx=idx+1;idx2=idx2+1;endidx=idx1+2*dif;end.学习帮手..专业整理.n=n/2;endAmp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.',’MarkerSize’,18);xlabel('FFTk');ylabel('FFTAmplification' );hold on;%Plottheamplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)' );holdoff;end命令窗口中的运行及其结果 ：DFT：x=[8,7,9,5,1,7,9,5];DFT(8,x);.学习帮手..专业整理.60 350 240)1nandoat(ae30g0fnlpAmeAsah20P-110 -202468-3246800kk图的幅度谱和相位谱（弧度制）FFT：x=[8,7,9,5,1,7,9,5];myFFT(8,x);603502)n40a1ndoat(aecginl300pAmeAsTahFPF20T-1FF10-20-30246802468FFTkFFTk图3.2.2FFT算法的幅度谱和相位谱（弧度制）【分析】DFT是离散信号、离散频谱之间的映射 。在这里我们可以看到序列的频谱也被离散化 。事实上，我们可以循着 DFT构造的方法验证这个频谱 ：.学习帮手..专业整理.首先，将序列做 N=8周期延拓，成为离散周期信号 。然后利用 DFS计算得到延拓后的频谱 ：N127jknjknX()x(n)e8x(n)e4n0n0

51,7,9 4i,7,3,6,9 4i,7 n 0,1,...,7，从而取 DFS以上作N 8延拓的主值区间得到 DFT，与图一致。因此计算正确。而对于FFT，我们可以看到它给出和 DFT一样的结果，说明了FFT算法就是 DFT的一个等价形式 。不过，由于序列不够长 ，FFT在计算速度上的优越性尚未凸显 。3.3已知连续时间信号x(t)=3cos8πt,X(ω)=3[(8)(8)]，该信号从t=0开始以采样周期Ts=0.1s进行采样得到序列x(n)，试选择合适的采样点数，分别采用DFT和FFT求其离散傅里叶变换X(k)的幅度、相位图，并将结果与X(k)的幅度、相位图，并将结果与X(ω)相比较。【解答】思路：此题与下一题都是一样的操作，可以在编程时统一用变量g（0或1）来控制是否有白噪声。这里取g=0（无白噪声）。另外，分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏的关系。M文件源代码：i)采样函数：function xs=sampJune3(N,Ts,g)%ThisisafunctionappliedtoProblem3&4.%N:numberofsamplingpoints;Ts:samplingperiod;g=0:Nogaussiannoise;g=1:gussiannoiseexists.n=1:N;fori=1:N%Notethatimuststartat0inanalysis.ThusImadeaadaptation.sample(i)=3*cos(8*pi*Ts*(i-1))+g*randn;.学习帮手..专业整理.%InMatlab,indexstartsat1,whichisnotconsistentwithourhabit.ThusImadeaadaptationincodes.%Itisasamplingprocesswith(g=1)/without(g=0)noise.endxs=sample;plot(n-1,sample, '.',’MarkerSize’,18);xlabel('n');ylabel('value');hold off;Mustuse(n-1),becauseinMatlab,indexstartsat1.endii)DFT和基2-FFT函数的代码，请见第3.2节。不需再新编一个。命令窗口中的运行及其结果 ：点采样：>>xs=sampJune3(12,0.1,0);% 末尾的0表示无噪声。DFT(12,xs);myFFT(12,xs);.学习帮手..专业整理.321eula 0v-1-2-30 2 4 6 8 10 12n图进行12点采样得到的序列202.5182161.514)1nan12d0.5aort(aecf10g0nlpAmeA8s-0.5ah6P-14-1.52-20-2.5024681012024681012kk图3.3.2DFT的幅度谱和相位谱（弧度制），出现了泄漏.学习帮手..专业整理.202.5182161.514)1nanidoai12r0.5t(aecgiin100pAmeAsTa8h-0.5FFPT6F-1F4-1.52-20-2.5024681012024681012FFTkFFTk图3.3.3FFT的幅度谱和相位谱（弧度制）。出现了频谱泄漏。点采样：>>xs=sampJune3(20,0.1,0);% 末尾的0表示无噪声。DFT(20,xs);myFFT(20,xs);321eula 0v-1-2-30 2 4 6 8 10 12 14 16 18 20n图进行20点采样得到的序列.学习帮手..专业整理.354303252)nandoari201aeclgfnlpA0m15eAsahP10-15-205101520-3510152000kk图3.3.5DFT的幅度谱和相位谱（弧度制）。频谱无泄漏。35430325)2nandoa1irt20(aeignl0pAmeA15sTahFP-1FTF10F-25-305101520-4510152000FFTkFFTk图3.3.6FFT的幅度谱和相位谱（弧度制）。频谱无泄漏。28点采样：>>xs=sampJune3(28,0.1,0);% 末尾的0表示无噪声。DFT(28,xs);myFFT(28,xs);.学习帮手..专业整理.321eul0av-1-2-3051015202530n图3.3.7进行28点采样得到的序列404353302)n25anidoair1t(aecl20gfnpA0meAs15ahP10-15-20102030-310203000kk图的幅度谱和相位谱（弧度制） 。再次出现频谱泄漏。.学习帮手..专业整理.40435330)2nanido25airt(ae1cligin20pAmeAsTa0hF15PFTFF-110-250102030-310203000FFTkFFTk图的幅度谱和相位谱（弧度制） 。再次出现频谱泄漏。【分析】分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏之间的关系 。现在与原信号频谱X( )比较后可以得出如下结论 ：图原信号的频谱（由两个冲激函数组成）原信号的频谱是X()3[(8)(8)]，在±8π上各有一强度为3π的谱线，在其余频率上为0。可见原信号被0.1s采样周期的采样信号离散化之后，谱线以20π为周期重复，并且只在（20k±8）π(k为整数）处非0。那么，在20点DFT（采样时间原信号周期的整数倍）中，只有第8根、第12根谱.学习帮手..专业整理.线非0。而在12点、28点DFT中，由于采样时间不是原信号周期的整数倍 ，谱线将向两边泄漏 。不过，对比12点采样和 28点采样，我们还可以发现 ，28点采样频谱的主谱线高度是次谱线高度的4倍，儿12点采样频谱的主谱线高度是次谱线高度的 3倍。可见，在无法保证采样时间是信号周期整数倍的情况下，增加采样时间有助于减轻频谱泄漏的程度 。3.4对第3步中所述连续时间信号叠加高斯白噪声信号 ，重复第3步过程。【解答】思路：此题与上一题都是一样的操作 ，可以在编程时统一用变量 g（0或1）来控制是否有白噪声 。这里取g=1（有白噪声）。另外，仍然分别取 12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏的关系 。M文件源代码：不需要再新编程序 。可以直接引用上面的函数 ：sampJune3(N,Ts,g)，取g=1，以体现存在白噪声DFT(N,x)myFFT(N,x)命令窗口中的运行及其结果 ：点采样：>>xs=sampJune3(12,0.1,1);% 末尾的1表示有噪声。DFT(12,xs);myFFT(12,xs);.学习帮手..专业整理.notiacfiilpmA

eulav1816141210864200 2

543210-1-2-3-4-50 2 4 6 8 10 12n图进行12点采样得到的含噪声的序列32)1nadar(eg0nAesah-1-24681012-3246810120kk图含噪声序列 DFT的幅度谱和相位谱（弧度制） 。.学习帮手..专业整理.18316214)nan12i1doat(aec10giin0pAmeA8sTahFPF6T-1FF4-220-3024681012024681012FFTkFFTk图3.4.3含噪声FFT的幅度谱和相位谱（弧度制）。点采样：>>xs=sampJune3(20,0.1,1);% 末尾的1表示有噪声。DFT(20,xs);myFFT(20,xs);4321eula 0v-1-2-3-40 2 4 6 8 10 12 14 16 18 20n图进行20点采样得到的含噪声序列.学习帮手..专业整理.30325220)1nanidoar(tae15g0fnlpAmeAsah10P-15-20-30510152005101520kk图3.4.5含噪声DFT的幅度谱和相位谱（弧度制）。303252)nan201doairt(aecligfn150pAmeAsTahFPF10T-1FF5-20-30510152005101520FFTkFFTk图3.4.6含噪声FFT的幅度谱和相位谱（弧度制）。.学习帮手..专业整理.点采样：>>xs=sampJune3(28,0.1,0);% 末尾的1表示有噪声。DFT(28,xs);myFFT(28,xs);4321eula 0v-1-2-3-40 5 10 15 20 25 30n图进行28点采样得到的含噪声序列.学习帮手..专业整理.403530n 25oitacifil 20pmA1510500403530noi 25tacfiilp 20mATFF 1510500【分析】

321)naida 0(elgnA-1sahP-2-3102030-41020300kk图含噪声DFT的幅度谱和相位谱（弧度制）。43)2naidar(e1lgnAesa0hPTFF-1-2102030-31020300FFTkFFTk图3.4.9含噪声FFT的幅度谱和相位谱（弧度制）。.学习帮手..专业整理.依然分别取 12点、20点、28点采样。仍然与原信号的频谱 X( ) 3[ ( 8) ( 8)]（图3.3.10）比较，可以得到结论：由于叠加了噪声，所以频谱都受到了一定的干扰。由于白噪声在各个频率的功率相等，因此频谱上各处的干扰也是均匀随机的。不过，通过对比我们可以发现，20点采样（无噪声时不发生泄漏的采样方法）在存在噪声时，仍然可以明显区分出原信号的谱线。第二好的是28点采样，因为采样时间较长，即使存在频谱泄漏也能较好地区分原信号的谱线。而最差的是12点采样，由于噪声的存在和严重的频谱泄漏，它的次谱线与主谱线的高度相差不大，使原信号不明显。j2k3.5已知序列x(n)4(n)3(n1)2(n2)(n3)，X(k)是x(n)的6点DFT，设WNkeN。（1）若有限长序列y(n)的6点DFT是Y(k)4kX(k)，求yn。W6()（2）若有限长序列w(n)的6点DFTW(k)是X(k)的实部，求w(n)。（3）若有限长序列q(n)的3点DFT是Q(k)X(2k)，k=0,1,2，求q(n)。【解答】思路：这是对DFT进行变换后求IDFT。考虑到IDFT和DFT定义的对称性，可以在DFT的基础上略加调整既可用于计算。首先，∵x(n)4(n)3(n1)2(n2)(n3)，∴它的6点采样是序列是x(n){4,3,2,1,0,0}。值得指出的是，在Matlab中，数组的序号是从1开始的（而在信号分析中习惯从0开始），不过我在上面编程时已考虑到这一情况，具体可见实验报告最后的“附录”。首先生成x(n)的6点DFT，再按照各小题分别转换 ，最后求相应的 IDFT。M文件源代码：输出x(n)的6点DFT的函数：.学习帮手..专业整理.function X=ExportDFT(N,x)%ThisisaDFTfunctionthatexportsthesolutiontovectorX.k=(0:N-1); %DefinevariablekforDFT.X=zeros(size(k)); %DefinetheinitialvalvesofX.fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i); %ItisthedefinitionofDFT.endendii)算第（1）小题的Y(k)的函数：function Y=Convertor1(X)%Thisisamathematicalconvertorforthesubproblem(1).fork=1:6Y(k)=exp((-j)*2*pi*(-4*(k-1))/6)*X(k);%Herewemustuse(k-1),becauseinMatlabindexstartsat1.endend算第（2）小题的W(k)的函数：function W=Convertor2(X)%Thisisamathematicalconvertorforthesubproblem(2).W=real(X);%AcquiretherealpartofX.endiv)算第（3）小题的Q(k)的函数：.学习帮手..专业整理.function Q=Convertor3(X)%Thisisamathematicalconvertorforthesubproblem(3).Q=zeros(3);% Detailedexplanationgoesherefortmp=1:3Q(tmp)=X(2*tmp);endendv)输出IDFT的函数：function x=ExportIDFT(N,X)%ThisistheIDFTfunctionformyexperiment.n=(0:N-1);%DefinevariablenforIDFT.x=zeros(size(n));%Definetheinitialvalvesofx.fork=0:N-1x=x+X(k+1)*exp(j*2*k*pi/N*n);endx=x/N;a=real(x);%WeMUSTusereal(x),thoughwemayALREADYknowxisreal.%It'scausedbynumericcalculation(notanalyticcalculation)inMatlab.stem(n,a,'.','MarkerSize',18);xlabel('n');ylabel('Solution');end命令窗口中的运行及其结果 ：>>x=[4,3,2,1,0,0];.学习帮手..专业整理.>>X=ExportDFT(6,x);第（1）小题Y=Convertor1(X);y=ExportIDFT(6,Y)y=Columns1through40.0000+0.0000i 0.0000+0.0000i 4.0000+0.0000i 3.0000+0.0000i%虚部都是 0，说明是实数Columns5through62.0000+0.0000i 1.0000-0.0000i %虚部都是 0，说明是实数%事实上，在Matlab 中，由于数值计算的 截断误差，对原复数做乘法后 ，答案的虚部可能有一极小的量。43.532.5n2ouo1.5S10.50-0.50 1 2 3 4 5n

答案：y(n)={0,0,4,3,2,1}图输出的y(n)，这是对x(n)的圆周移位。第（2）小题W=Convertor2(X);w=ExportIDFT(6,W)w=.学习帮手..专业整理.Columns1through44.0000+0.0000iColumns5through61.0000+0.0000i%事实上，在Matlab量。43.532.5notiul 2oS1.510.500 1 2

1.5000+0.0000i 1.0000+0.0000i 1.0000+0.0000i%虚部都是 0，说明是实数1.5000+0.0000i% 虚部都是0，说明是实数；中，由于数值计算的截断误差 ，对原复数做乘法后 ，答案的虚部可能有一极小的答案：w(n)={0,0,4,3,2,1}图输出的w(n)。3 4 5n第（3）小题Q=Convertor3(X);q=ExportIDFT(6,Q)q=Columns1through41.5000-0.0000i-0.1667-0.2887i0.7500-1.2990i0.8333-0.0000iColumns5through6-0.5000-0.8660i 1.0833-1.8764i这里的答案都是幅值 、相位均非 0的复数，而教材（实验指导第 109页）并未要求作图，这里略去。.学习帮手..专业整理.答案：q(n)={1.5，-0.1667-0.2887i，0.7500-1.2990i，0.8333，-0.5000-0.8660i，1.0833-1.8764i}【分析】对原序列进行DFT运算后，可以得到X(k)={10，3.5-4.33i，2.5-0.866i，2，2.5+0.866i，3.5+4.330i}。第（1）小题，根据DFT的性质可以判断，对原频谱乘上旋转因子W64k之后进行IDFT得到的y(n)，就是对原序列做圆周移位：y(n)15X(k)ejk0n0,0,4,3,2,1。Nk0第（2）小题，由于对原频谱取了实部，那么根据DFT的奇偶虚实性知，得到的w(n)也是实数的。第（3）小题，对原信号进行了尺度变换（“抽取”），导致丢失了一些谱线，使得无法通过IDFT得到原来的序列()。说明频谱记录了原有信号的信息，若频谱发生变化，则对应的时域信号也随之改变。xn3.6已知信号x(t) sin(2f1t) sin(2f2t) sin(2f3t)，其中f1=4Hz、f2=4.02Hz、f3=5Hz，采用采样频率为20Hz进行采样，求1）当采样长度N分别为512和2048情况下x(t)的幅度频谱；．．．．2）当采样长度N为32，且增补N个零点、4N个零点、8N个零点、16N个零点情况下x(t)的幅度频．．．谱。．【解答】思路：采样是有限且离散的 ，用DFT（FFT算法）计算频谱，以便得到离散的频谱 ，并且具有较高速度 。20Hz对应的采样周期 Ts=0.05s。M文件源代码：i)采样函数（其中Plus表示采样后补零的个数 ）function xs=sampJune6(N,Plus).学习帮手..专业整理.%ThisisafunctionappliedtoProblem6%N:samlingpoints;Plus:quantityofadditinalzeropoints.Ts=1/20;w1=2*pi*4;w2=2*pi*4.02;w3=2*pi*5;sample=zeros(1,N+Plus);n=(1:N);sample=sin(w1*Ts*(n-1))+sin(w2*Ts*(n-1))+sin(w3*Ts*(n-1));forp=(N+1):(N+Plus)sample(p)=0;%Addzeropoints.endxs=sample;%Returnendii)由于只要求显示幅度频谱 ，所以删去 myFFT函数中绘制相位频谱的命令 ，使它的最后部分修改如下 ：原来的：function myFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification' );hold on;%Plottheamplification..学习帮手..专业整理.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)' );holdoff;end修改后的：function myFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification' );%Plottheamplification.end命令窗口中的运行及其结果 ：第（1）小题x512=sampJune6(512,0);x2048=sampJune6(2048,0);myFFT(512,x512);myFFT(2048,x2048);.学习帮手..专业整理.512Points300250200ntacfiilmAFF1005000 50 100 150 200 250 300 350 400 450 500FFTk图进行512点采样得到的频谱2048Points12001000800noc600pATF40020000200400600800100012001400160018002000FFTk图进行2048点采样得到的频谱第（2）小题>>x32p1N=sampJune6(32,32*1);%32点采样，补零N=32个，共64个数据点>>x32p4N=sampJune6(32,32*4);%32点采样，补零4N=128个，共160个数据点>>x32p8N=sampJune6(32,32*8);%32点采样，补零8N=256个，共288个数据点>>x32p16N=sampJune6(32,32*16);%32点采样，补零16N=128个，共544个数据点>>myFFT(32+32*1,x32p1N);>>myFFT(32+32*4,x32p4N);.学习帮手..专业整理.myFFT(32+32*8,x32p8N);myFFT(32+32*16,x32p16N);353025n tafilpmAF151050010203040506070FFTk图采样N=32点，补零N=32点，共64点的频谱35302520afilpmATFF

151050020406080100120140160FFTk图采样N=32点，补零 4N=128点，共160点的频谱.学习帮手..专业整理.353025n taficilpmATFF 151050050100150200250300FFTk图采样N=32点，补零 8N=32点，共288点的频谱353025ni 20fiilpmATFF 10500 100 200 300 400 500 600FFTk图采样N=32点，补零16N=32点，共544点的频谱【分析】请注意，题目只要求绘制幅度频谱 。第（1）小题：首先，由于采样时间都不是原有信号周期的整数倍 ，两个采样方式对应的频谱均发生了泄漏 。不过由于2048点采样对应的采样时间较长 ，它频谱泄漏的程度比 512点采样轻。其次，由于20Hz的2048点采样的频率分辨率为 20/2048=0.0098Hz<0.2Hz ，因此放大频谱图之.学习帮手..专业整理.后我们可以看到 4Hz、4.02Hz和5Hz对应的谱线，而512点采样的频率分辨率为 20/512=0.039Hz>0.2Hz，因此4Hz和4.02Hz对应的谱线无法区分 。第（2）小题：首先，由于采样时间都不是原有信号周期的整数倍 ，频谱均发生了泄漏 。而且由于采样时间较短 ，频谱泄漏比第（1）小题的两个频谱更加严重 。其次，由于都是 32点采样，因此实际的频率分辨率较低 ，无法区分 4Hz和4.02Hz对应的谱线。最后，虽然都是 32点采样，但由于补 0个数的不同，各频谱谱线间距各不相同 。例如，补零最多的序列一共有 544个数据点，因此谱线间距小 。由此还可以得出结论 ：数据点个数越多 ，则频谱越倾向于连续。可见，当采样时间不是原信号周期整数倍而且采样时间较短时 ，频谱泄漏相当严重的 ，所有的频率上都有了幅值即能量 ，可见当取样信号的样点数取得不够时 ，原信号所携带的信息就不能被完全取得 。而若将取样信号补零 ，由图可见信号的能量相应的泄漏到了几乎所有频率上了 ，这样所得的信号仍然严重失真，因此不能靠将信号补零这样的方法来取得更精确的采样信号 。要想获得不泄漏的频谱 ，在采样频率不变的前提下 ，必须使采样时间等于原信号周期的整数倍 ，或者尽量延长采样时间以减少泄漏 。四、实验体会4.1关于各个实验的分析，请见第3部分每道题的末尾。4.2在Matlab中，数组的序号是从1开始的，这与信号处理时通常的序号起点（0）不一致。我在编程充分注意到了这个问题。4.3由于Matlab进行数值计算的过程中存在截断误差，所以最后算得的值并不是准确值。例如，对一个复数z，即使f(z)的虚部为零，但由于截断误差的存在（特别是z的虚部为无穷小数时），最终f(z)值的虚部可能是一个极小的非零值，从而在显示时出现“零虚部”（例如，2+0.0000i)。.学习帮手..专业整理.4.4通过利用软件对离散信号的各种变换 DTFT、DFT以及其快速算法 FFT进行计算，使得在实验中比较难以实现的信号分析过程 （离散信号的采集和显示都是比较困难的 ）在计算机计算中实现 ，证明了理论的正确性，说明仿真计算是一种十分有效的辅助手段 。4.5通过这次实验和上次实验 《信号的采集与恢复 》我知道了，要想尽量不失真地取得一个信号的频谱 （低混叠、低泄漏），应该尽量满足以下条件 ：1）使用的开关函数要尽量接近理想冲激串；2）采样频率要高于原始信号的奈奎斯特频率。对于频谱不受限的信号，为了避免频谱混叠，应该使用低通滤波器进行滤波；3）对于频带不受限的信号，抗混叠滤波器要尽量接近理想滤波器。4）采样的持续时间最好能够是原信号周期的整数倍，一避免频谱泄漏。而当不知道原信号的周期（或者周期不稳定）时，就要通过延长采样时间来尽量减少泄漏，从而突出原信号的谱线。5）当信号混有白噪声时，就更应注意减少频谱的泄漏和混叠，否则信号分析更加困难，甚至可能会使原信号被误差“淹没”。6）若原信号有多个频率成分，应该尽量提高采样的频率分辨率，以区分出更细微的频率差异。4.6在实验中，在计算2048点采样时，初步体会到了 FFT算法的优越性 。在我的计算机上 ，FFT算法的确比原始的DFT更快。不过由于采样点数较少 ，这一差别仅限于几秒钟 。在采样点更多时 ，FFT在速度上的优越性应该能进一步突出 。4.7实验中遇到的问题及其解决 ：实验中有些 M 文件代码总是出错 。解决方