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

1、第 1 页 共 27 页实 验 报 告课程名称： 信号分析与处理 指导老师： 成绩：_实验名称：离散傅里叶变换和快速傅里叶变换 实验类型： 基础实验 同组学生姓名： 第二次实验 离散傅里叶变换和快速傅里叶变换一、实验目的1.1 掌握离散傅里叶变换（DFT）的原理和实现；1.2 掌握快速傅里叶变换（FFT）的原理和实现，掌握用 FFT 对连续信号和离散信号进行谱分析的方法。1.3 会用 Matlab 软件进行以上练习。二、实验原理2.1 关于 DFT 的相关知识序列 x(n)的离散事件傅里叶变换（ DTFT）表示为，njnjexeX)()(如果 x(n)为因果有限长序列， n=0,1,.,N-1

2、，则 x(n)的 DTFT 表示为，njNnjee10x(n)的离散傅里叶变换（DFT）表达式为，)1,.()()(210kexkXnNjNn序列的 N 点 DFT 是序列 DTFT 在频率区间0,2 上的 N 点灯间隔采样，采样间隔为 2/N。通过DFT，可以完成由一组有限个信号采样值 x(n)直接计算得到一组有限个频谱采样值 X(k)。X (k)的幅度谱为，其中下标 R 和 I 分别表示取实部、虚部的运算。X(k)的相位谱为)()(22kXkIR。)(arctnRI离散傅里叶反变换（IDFT）定义为。)1,.0()(1)(20 NnekXNnxNj2.2 关于 FFT 的相关知识快速傅里叶

3、变换（FFT）是 DFT 的快速算法，并不是一个新的映射。 FFT 利用了 函数的周期nNje2性和对称性以及一些特殊值来减少 DFT 的运算量，可使 DFT 的运算量下降几个数量级，从而使数字信号装 订 线第 2 页 共 27 页处理的速度大大提高。若信号是连续信号，用 FFT 进行谱分析时，首先必须对信号进行采样，使之变成离散信号，然后就可以用 FFT 来对连续信号进行谱分析。为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器，且抗混叠滤波器的截止频率不得高于与采样频率的一半。比较 DFT 和 IDFT 的定义，两者的区别仅在于指数因子的指数部分的符号差异和幅度尺度变换，因此可用

4、FFT 算法来计算 IDFT。3、实验内容与相关分析（共 6 道）说明：为了便于老师查看，现将各题的内容写在这里题目按照 3.1、3.2、.、3.6 排列。每道题包含如下内容：题干、解答（思路、M 文件源代码、命令窗口中的运行及其结果） 、分析。其中“命令窗口中的运行及其结果”按照小题顺序排列，各小题包含命令与结果（图形或者序列） 。3.1 求有限长离散时间信号 x(n)的离散时间傅里叶变换（DTFT）X(e j )并绘图。（1）已知 ；（2）已知 。其 他01)(nx 102)(nx【解答】思路：这是 DTFT 的变换，按照定义编写 DTFT 的 M 文件即可。考虑到自变量 是连续的，为了方

5、便计算机计算，计算时只取三个周期-2,4 中均匀的 1000 个点用于绘图。理论计算的各序列 DTFT 表达式，请见本题的分析。M 文件源代码（我的 Matlab 源文件不支持中文注释，抱歉）：function DTFT(n1,n2,x)%This is a DTFT function for my experiment of Signal Processing %Define the bracket of omega for plotting.X=zeros(size(w);%Define the initial values of X.for i=n1:n2X=X+x(i-n1+1)*ex

6、p(-1)*j*w*i);%It is the definition of DTFT.endAmp=abs(X);%Acquire the amplification.Phs=angle(X);%Acquire the phase angle (radian).subplot(1,2,1);plot(w,Amp,r); xlabel(Omega);ylabel(Amplification);hold on;%Plot amplification on the left.subplot(1,2,2);plot(w,Phs,b);xlabel(Omega);ylabel(Phase Angle (

7、radian);hold off;%Plot phase angle on the right.end命令窗口中的运行及其结果（理论计算的各序列 DTFT 表达式，请见本题的分析）：第（1）小题 n=(-2:2); x=1.n; DTFT(-2,2,x);第 3 页 共 27 页-5 0 5 1000.511.522.533.544.55Amplification-5 0 5 10-4-3-2-101234Phase Angle(radian)第（2）小题 n=(0:10); x=2.n; DTFT(0,10,x);-5 0 5 10600800100012001400160018002000

8、2200Amplification-5 0 5 10-4-3-2-101234Phase Angle(radian)【分析】对于第（1）小题，由于序列 x(n)只在有限区间（-2 ，-1，-，1，2）上为 1，所以是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。事实上，我们可计算出它的表达式：，可见幅度频谱拥有主极 jjjnjnj XxX e)(e1e)()( 5522图 3.1.1 在-2,4 范围内 3 个周期的幅度谱和相位谱（弧度制）图 3.1.2 在-2,4 范围内 3 个周期的幅度谱和相位谱（弧度制）第 4 页 共 27 页大和次极大，两个主极大间有|5-1 |=4 个极小

9、，即有 3 个次级大。而对于它的相位频谱，则是周期性地在-、 0、 之间震荡。对于第（2）小题，由于是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。而它的表达式： ，因此主极大之间只有 jjnnjnj XxX e21)(e212ee)()(10|0-1|=1 个极小，不存在次级大。而对于它的相位频谱，则是在一个长为 2 的周期内有|11-1| =10 次振荡。而由 DTFT 的定义可知，频谱都是以 2 为周期向两边无限延伸的。由于 DTFT 是连续谱，对于计算机处理来说特别困难，因此我们才需要离散信号的频谱也离散，由此构造出 DFT（以及为加速计算 DFT的 FFT） 。3.2 已知

10、有限长序列 x(n)=8,7,9,5,1,7,9,5,试分别采用 DFT 和 FFT 求其离散傅里叶变换 X(k)的幅度、相位图。【解答】思路：按照定义编写 M 文件即可。M 文件源代码：i) DFT 函数：function DFT(N,x)%This is a DFT function for my experiment of Signal Processing %Define variable k for DFT.X=zeros(size(k);%Define the initial valves of X.for i=0:N-1X=X+x(i+1)*exp(-1)*j*2*k*pi/N*

11、i);%It is the definition of DFT.endAmp=abs(X);%Acquire the amplification.Phs=angle(X);%Acquire the phase angle (radian).subplot(1,2,1);stem(k,Amp,.,MarkerSize,18); xlabel(k);ylabel(Amplification);hold on;%Plot amplification on the left.subplot(1,2,2);stem(k,Phs,*);xlabel(k);ylabel(Phase Angle (radia

12、n);hold off;%Plot phase angle on the right.endii) 基 2-FFT 函数function myFFT(N,x)%This is a base-2 FFT function.lov=(0:N-1);j1=0;for i=1:N %indexed addressingif i1digit=digit+1;k=k/2;endn=N/2;% Now we start the “butterfly-shaped“ process.for mu=1:digitdif=2(mu-1);%Differnce between the indexes of the

13、target variables.idx=1;for i=1:nidx1=idx;idx2=1;for j1=1:N/(2*n)r=(idx2-1)*2(digit-mu);wn=exp(j*(-2)*pi*r/N);%It is the “circulating coefficients“.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;endn=n/2;endAmp=abs(x);%Acquire the amplifica

14、tion.Phs=angle(x);%Acquire the phase angle (radian).subplot(1,2,1);stem(lov,Amp,.,MarkerSize,18);xlabel(FFT k);ylabel(FFT Amplification);hold on;%Plot the amplification.subplot(1,2,2);stem(lov,Phs,*);xlabel(FFT k);ylabel(FFT Phase Angle (radian);hold off;end第 6 页 共 27 页命令窗口中的运行及其结果：DFT： x=8,7,9,5,1,

15、7,9,5; DFT(8,x);0 2 4 6 80102030405060kAmplification0 2 4 6 8-3-2-10123kPhase Angle(radian)FFT： x=8,7,9,5,1,7,9,5; myFFT(8,x);0 2 4 6 80102030405060FFT kFFTAmplification0 2 4 6 8-3-2-10123FFT kFFTPhase Angle(radian)【分析】DFT 是离散信号、离散频谱之间的映射。在这里我们可以看到序列的频谱也被离散化。事实上，我们可以循着 DFT 构造的方法验证这个频谱：首先，将序列做 N=8 周期

16、延拓，成为离散周期信号。然后利用 DFS 计算得到延拓后的频谱：图 3.2.1 DFT 的幅度谱和相位谱（弧度制）图 3.2.1 DFT 的幅度谱和相位谱（相位是弧度制的）图 3.2.2 FFT 算法的幅度谱和相位谱（弧度制）第 7 页 共 27 页，从而取 DFS 延 拓以 上 作 87,.10,496,37,51e)(e)()(7041082 NniinxnxXjkNnjk的主值区间得到 DFT，与图一致。因此计算正确。而对于 FFT，我们可以看到它给出和 DFT 一样的结果，说明了 FFT 算法就是 DFT 的一个等价形式。不过，由于序列不够长，FFT 在计算速度上的优越性尚未凸显。3.

17、3 已知连续时间信号 x(t)=3cos8t, X()= ，该信号从 t=0 开始以采样周期)8()(3Ts=0.1 s 进行采样得到序列 x(n)，试选择合适的采样点数，分别采用 DFT 和 FFT 求其离散傅里叶变换X(k)的幅度、相位图，并将结果与 X(k)的幅度、相位图，并将结果与 X()相比较。【解答】思路：此题与下一题都是一样的操作，可以在编程时统一用变量 g（0 或 1）来控制是否有白噪声。这里取 g=0（无白噪声） 。另外，分别取 12 点、20 点、28 点采样，以考察采样长度的选择与频谱是否泄漏的关系。M 文件源代码：i)采样函数：function xs=sampJune3

18、(N,Ts,g)%This is a function applied to Problem 3 Ts: sampling period; g=0: No gaussian noise; g=1: gussian noise exists.n=1:N;for i=1:N%Note that i must start at 0 in analysis. Thus I made a adaptation.sample(i)=3*cos(8*pi*Ts*(i-1)+g*randn;%In Matlab, index starts at 1, which is not consistent with

19、our habit. Thus I made a adaptation in codes.%It is a sampling process with(g=1)/without(g=0) noise.endxs=sample;plot(n-1,sample,.,MarkerSize,18);xlabel(n);ylabel(value);hold off;% Must use (n-1), because in Matlab, index starts at 1.endii)DFT和基2-FFT函数的代码，请见第3.2节。不需再新编一个。命令窗口中的运行及其结果：12 点采样： xs=samp

20、June3(12,0.1,0);%末尾的0表示无噪声。 DFT(12,xs); myFFT(12,xs);第 8 页 共 27 页0 2 4 6 8 10 12-3-2-10123nvalue0 2 4 6 8 10 1202468101214161820kAmplification0 2 4 6 8 10 12-2.5-2-1.5-1-0.500.511.522.5kPhase Angle(radian)图 3.3.1 进行 12 点采样得到的序列图 3.3.2 DFT 的幅度谱和相位谱（弧度制） ，出现了泄漏第 9 页 共 27 页0 2 4 6 8 10 1202468101214161

21、820FFT kFFTAmplification0 2 4 6 8 10 12-2.5-2-1.5-1-0.500.511.522.5FFT kFFTPhaseAngle (radian)20 点采样： xs=sampJune3(20,0.1,0);%末尾的0表示无噪声。 DFT(20,xs); myFFT(20,xs);0 2 4 6 8 10 12 14 16 18 20-3-2-10123nvalue图 3.3.3 FFT 的幅度谱和相位谱（弧度制） 。出现了频谱泄漏。图 3.3.4 进行 20 点采样得到的序列第 10 页 共 27 页0 5 10 15 200510152025303

22、5kAmplification0 5 10 15 20-3-2-101234kPhase Angle(radian)0 5 10 15 2005101520253035FFT kFFTAmplification0 5 10 15 20-4-3-2-101234FFT kFFTPhaseAngle (radian)28 点采样： xs=sampJune3(28,0.1,0);%末尾的 0 表示无噪声。 DFT(28,xs); myFFT(28,xs);图 3.3.5 DFT 的幅度谱和相位谱（弧度制） 。频谱无泄漏。图 3.3.6 FFT 的幅度谱和相位谱（弧度制） 。频谱无泄漏。第 11 页

23、共 27 页0 5 10 15 20 25 30-3-2-10123nvalue0 10 20 300510152025303540kAmplification0 10 20 30-3-2-101234kPhase Angle (radian)图 3.3.7 进行 28 点采样得到的序列图 3.3.8 DFT 的幅度谱和相位谱（弧度制） 。再次出现频谱泄漏。第 12 页 共 27 页0 10 20 300510152025303540FFT kFFTAmplification0 10 20 30-3-2-101234FFT kFFTPhaseAngle(radian)【分析】分别取 12 点、

24、20 点、28 点采样，以考察采样长度的选择与频谱是否泄漏之间的关系。现在与原信号频谱 比较后可以得出如下结论：)(X原信号的频谱是 ，在8 上各有一强度为 3 的谱线，在其余)8()(3)( X频率上为 0。可见原信号被 0.1 s 采样周期的采样信号离散化之后，谱线以 20 为周期重复，并且只在（20k8） (k 为整数）处非 0。那么，在 20 点 DFT（采样时间原信号周期的整数倍）中，只有第 8 根、第 12 根谱线非 0。而在 12 点、28 点 DFT 中，由于采样时间不是原信号周期的整数倍，谱线将向两边泄漏。不过，对比 12 点采样和 28 点采样，我们还可以发现，28 点采样

25、频谱的主谱线高度是次谱线高度的4 倍，儿 12 点采样频谱的主谱线高度是次谱线高度的 3 倍。可见，在无法保证采样时间是信号周期整数倍的情况下，增加采样时间有助于减轻频谱泄漏的程度。图 3.3.9 FFT 的幅度谱和相位谱（弧度制） 。再次出现频谱泄漏。图 3.3.10 原信号的频谱（由两个冲激函数组成）第 13 页 共 27 页3.4 对第 3 步中所述连续时间信号叠加高斯白噪声信号，重复第 3 步过程。【解答】思路：此题与上一题都是一样的操作，可以在编程时统一用变量 g（0 或 1）来控制是否有白噪声。这里取 g=1（有白噪声） 。另外，仍然分别取 12 点、20 点、28 点采样，以考察

26、采样长度的选择与频谱是否泄漏的关系。M 文件源代码：不需要再新编程序。可以直接引用上面的函数：sampJune3(N,Ts,g)，取 g=1，以体现存在白噪声DFT(N,x)myFFT(N,x)命令窗口中的运行及其结果：12 点采样： xs=sampJune3(12,0.1,1);%末尾的1表示有噪声。 DFT(12,xs); myFFT(12,xs);0 2 4 6 8 10 12-5-4-3-2-1012345nvalue图 3.4.1 进行 12 点采样得到的含噪声的序列第 14 页 共 27 页0 2 4 6 8 10 12024681012141618kAmplification0

27、2 4 6 8 10 12-3-2-10123kPhase Angle(radian)0 2 4 6 8 10 12024681012141618FFT kFFTAmplification0 2 4 6 8 10 12-3-2-10123FFT kFFTPhaseAngle (radian)20 点采样： xs=sampJune3(20,0.1,1);%末尾的1表示有噪声。 DFT(20,xs); myFFT(20,xs);图 3.4.2 含噪声序列 DFT 的幅度谱和相位谱（弧度制） 。图 3.4.3 含噪声 FFT 的幅度谱和相位谱（弧度制） 。第 15 页 共 27 页0 2 4 6 8

28、 10 12 14 16 18 20-4-3-2-101234nvalue0 5 10 15 20051015202530kAmplification0 5 10 15 20-3-2-10123kPhaseAngle(radian)图 3.4.4 进行 20 点采样得到的含噪声序列图 3.4.5 含噪声 DFT 的幅度谱和相位谱（弧度制） 。第 16 页 共 27 页0 5 10 15 20051015202530FFT kFFT Amplification0 5 10 15 20-3-2-10123FFT kFFT Phase Angle (radian)28 点采样： xs=sampJun

29、e3(28,0.1,0);%末尾的1表示有噪声。 DFT(28,xs); myFFT(28,xs);0 5 10 15 20 25 30-4-3-2-101234nvalue图 3.4.6 含噪声 FFT 的幅度谱和相位谱（弧度制） 。图 3.4.7 进行 28 点采样得到的含噪声序列第 17 页 共 27 页0 10 20 300510152025303540kAmplification0 10 20 30-4-3-2-10123kPhaseAngle(radian)0 10 20 300510152025303540FFT kFFT Amplification0 10 20 30-3-2-

30、101234FFT kFFT Phase Angle (radian)【分析】依然分别取 12 点、20 点、28 点采样。仍然与原信号的频谱图 3.4.8 含噪声 DFT 的幅度谱和相位谱（弧度制） 。图 3.4.9 含噪声 FFT 的幅度谱和相位谱（弧度制） 。第 18 页 共 27 页（图 3.3.10）比较，可以得到结论：)8()(3)( X由于叠加了噪声，所以频谱都受到了一定的干扰。由于白噪声在各个频率的功率相等，因此频谱上各处的干扰也是均匀随机的。不过，通过对比我们可以发现，20 点采样（无噪声时不发生泄漏的采样方法）在存在噪声时，仍然可以明显区分出原信号的谱线。第二好的是 28

31、点采样，因为采样时间较长，即使存在频谱泄漏也能较好地区分原信号的谱线。而最差的是 12 点采样，由于噪声的存在和严重的频谱泄漏，它的次谱线与主谱线的高度相差不大，使原信号不明显。3.5 已知序列 ，X(k) 是 x(n)的 6 点 DFT，设 。3)2()1(3)4( nnnx kNjkeW2（1）若有限长序列 y(n)的 6 点 DFT 是 ，求 y(n)。()(46WkYk（2）若有限长序列 w(n)的 6 点 DFT W(k)是 的实部，求 w(n)。（3）若有限长序列 q(n)的 3 点 DFT 是 ，k=0,1,2，求 q(n)。)2(XQ【解答】思路：这是对 DFT 进行变换后求

32、IDFT。考虑到 IDFT 和 DFT 定义的对称性，可以在 DFT 的基础上略加调整既可用于计算。首先， ，)3()2()1(3)4( nnnx 它的 6 点采样是序列是 。值得指出的是，在 Matlab 中，数组的序号是从 1 开始的0,4x（而在信号分析中习惯从 0 开始） ，不过我在上面编程时已考虑到这一情况，具体可见实验报告最后的“附录” 。首先生成 x(n)的 6 点 DFT，再按照各小题分别转换，最后求相应的 IDFT。M 文件源代码：i) 输出x( n)的6点DFT的函数：function X = ExportDFT(N,x)%This is a DFT function th

33、at exports the solution to vector X.k=(0:N-1); %Define variable k for DFT.X=zeros(size(k); %Define the initial valves of X.for i=0:N-1X=X+x(i+1)*exp(-1)*j*2*k*pi/N*i); %It is the definition of DFT.endendii)算第（1）小题的 Y(k)的函数：function Y = Convertor1(X)%This is a mathematical convertor for the subproble

34、m (1).for k=1:6Y(k)=exp(-j)*2*pi*(-4*(k-1)/6)*X(k);%Here we must use (k-1),because in Matlab index starts at 1.end第 19 页 共 27 页endiii)算第（2）小题的 W(k)的函数：function W = Convertor2(X)%This is a mathematical convertor for the subproblem (2).W=real(X);%Acquire the real part of X.endiv)算第（3）小题的 Q(k)的函数：funct

35、ion Q = Convertor3(X)%This is a mathematical convertor for the subproblem (3).Q=zeros(3);% Detailed explanation goes herefor tmp=1:3Q(tmp)=X(2*tmp);endendv)输出 IDFT 的函数：function x = ExportIDFT(N,X)%This is the IDFT function for my experiment.n=(0:N-1);%Define variable n for IDFT.x=zeros(size(n);%Defi

36、ne the initial valves of x.for k=0:N-1x=x+X(k+1)*exp(j*2*k*pi/N*n);endx=x/N;a=real(x);%We MUST use real(x), though we may ALREADY know x is real. %Its caused by numeric calculation (not analytic calculation) in Matlab.stem(n,a,.,MarkerSize,18);xlabel(n);ylabel(Solution);end命令窗口中的运行及其结果： x=4,3,2,1,0,

37、0; X=ExportDFT(6,x);第（1）小题 Y=Convertor1(X); y=ExportIDFT(6,Y)y =Columns 1 through 40.0000 + 0.0000i 0.0000 + 0.0000i 4.0000 + 0.0000i 3.0000 + 0.0000i%虚部都是 0，说明是实数Columns 5 through 62.0000 + 0.0000i 1.0000 - 0.0000i %虚部都是 0，说明是实数%事实上，在 Matlab 中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。第 20 页 共 27 页0 1 2 3

38、 4 5-0.500.511.522.533.54nSolution第（2）小题 W=Convertor2(X); w=ExportIDFT(6,W)w =Columns 1 through 44.0000 + 0.0000i 1.5000 + 0.0000i 1.0000 + 0.0000i 1.0000 + 0.0000i%虚部都是 0，说明是实数Columns 5 through 61.0000 + 0.0000i 1.5000 + 0.0000i %虚部都是 0，说明是实数；%事实上，在 Matlab 中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。0 1 2

39、 3 4 500.511.522.533.54nSolution第（3）小题 Q=Convertor3(X); q=ExportIDFT(6,Q)q =Columns 1 through 41.5000 - 0.0000i -0.1667 - 0.2887i 0.7500 - 1.2990i 0.8333 - 0.0000iColumns 5 through 6-0.5000 - 0.8660i 1.0833 - 1.8764i答案：y(n)=0 ,0,4,3,2,1图 3.5.1 输出的 y(n)，这是对 x(n)的圆周移位。答案：w(n)=0,0,4,3,2,1图 3.5.2 输出的 w(

40、n)。第 21 页 共 27 页这里的答案都是幅值、相位均非 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 的性质可以判断，对原频谱乘上旋转因子 之后进行 IDFT 得到的 y(n)，kW46就是对原序列做圆周移位： 。1,2340,

41、e)(1)(500knjkNny第（2）小题，由于对原频谱取了实部，那么根据 DFT 的奇偶虚实性知，得到的 w(n)也是实数的。第（3）小题，对原信号进行了尺度变换（“抽取” ） ，导致丢失了一些谱线，使得无法通过 IDFT 得到原来的序列 x(n)。说明频谱记录了原有信号的信息，若频谱发生变化，则对应的时域信号也随之改变。3.6 已知信号 ，其中 f1=4 Hz、f 2=4.02 Hz、f 3=5 Hz，采用采样频)2sin()si()2si() 321 tftftft 率为 20 Hz 进行采样，求（1）当采样长度 N 分别为 512 和 2048 情况下 x(t)的幅度频谱；（2）当采

42、样长度 N 为 32，且增补 N 个零点、4N 个零点、8N 个零点、16N 个零点情况下 x(t)的幅度频谱。【解答】思路：采样是有限且离散的，用 DFT（FFT 算法）计算频谱，以便得到离散的频谱，并且具有较高速度。20 Hz 对应的采样周期 Ts=0.05s。M 文件源代码：i)采样函数（其中 Plus 表示采样后补零的个数）function xs=sampJune6(N,Plus)%This is a function applied to Problem 6%N:samling points;Plus:quantity of additinal zero points.Ts=1/20

43、;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);for p=(N+1):(N+Plus)sample(p)=0;%Add zero points.endxs=sample;%Returnendii)由于只要求显示幅度频谱，所以删去 myFFT 函数中绘制相位频谱的命令，使它的最后部分修改如下：原来的：function myFFT(N,x)第 22 页 共 27 页%This is a base-2 FFT

44、 function wrote by myself.Amp=abs(x);%Acquire the amplification.Phs=angle(x);%Acquire the phase angle (radian).subplot(1,2,1);stem(lov,Amp,.);xlabel(FFT k);ylabel(FFT Amplification);hold on;%Plot the amplification.subplot(1,2,2);stem(lov,Phs,*);xlabel(FFT k);ylabel(FFT Phase Angle (radian);hold off;

45、end修改后的：function myFFT(N,x)%This is a base-2 FFT function wrote by myself.Amp=abs(x);%Acquire the amplification.stem(lov,Amp,.);xlabel(FFT k);ylabel(FFT Amplification);%Plot the amplification.end命令窗口中的运行及其结果：第（1）小题 x512=sampJune6(512,0); x2048=sampJune6(2048,0); myFFT(512,x512); myFFT(2048,x2048);0

46、50 10 150 20 250 30 350 40 450 50050101502025030FT kFT Amplification512 Points图 3.6.1 进行 512 点采样得到的频谱第 23 页 共 27 页0 20 40 60 80 100 120 140 160 180 200020406080100120FT kFT Amplification2048 Points第（2）小题 x32p1N=sampJune6(32,32*1); %32点采样，补零N=32 个，共64 个数据点 x32p4N=sampJune6(32,32*4); %32点采样，补零4N=128个，

47、共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);0 10 20 30 40 50 60 7005101520253035FT kFT Amplification图 3.6.2 进行 2048 点采样得到的频谱图 3.6.3

48、采样 N=32 点，补零 N=32 点，共 64 点的频谱第 24 页 共 27 页0 20 40 60 80 10 120 140 16005101520253035FT kFT Amplification0 50 10 150 20 250 3005101520253035FT kFT Amplification图 3.6.4 采样 N=32 点，补零 4N=128 点，共 160 点的频谱图 3.6.5 采样 N=32 点，补零 8N=32 点，共 288 点的频谱第 25 页 共 27 页0 10 20 30 40 50 6005101520253035FT kFT Amplification【分析】请注意，题目只要求绘制幅度频谱。第（1）小题：首先，由于采样时间都不是原有信号周期的整数倍，两个采样方式对应的频谱均发生了泄漏。不过由于 2048 点采样对应的采样时间较长，它频谱泄漏的程度比 512 点采样轻。其次，由于 20 Hz 的 2048 点采样的频率分辨率为 20/2048=0.0098 Hz 0.2 Hz，因此 4 Hz 和 4.02 Hz 对应的谱线无法区分。第（2）小题：首先，由于采样时间都不是原有信号周期