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

实验报告

课程名称： 信号分析与处理 指导老师： 成绩：__________________

实验名称：离散傅里叶变换和快速傅里叶变换 实验类型： 基础实验 同组学生姓名：

第二次实验 离散傅里叶变换和快速傅里叶变换

一、实验目的

1.1掌握离散傅里叶变换（DFT）的原理和实现；

1.2掌握快速傅里叶变换（FFT）的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法。 1.3 会用Matlab软件进行以上练习。

二、实验原理

2.1关于DFT的相关知识

序列x(n)的离散事件傅里叶变换（DTFT）表示为

X(e)?装 j?n????x(n)e??j?n，

如果x(n)为因果有限长序列，n=0,1,...,N-1，则x(n)的DTFT表示为

订 j?X(e)??x(n)e?j?n，

n?0N?1线 x(n)的离散傅里叶变换（DFT）表达式为

X(k)??x(n)en?0N?1?j2?nkN(k?0,1,...,N?1)，

序列的N点DFT是序列DTFT在频率区间[0,2π]上的N点灯间隔采样，采样间隔为2π/N。通过DFT，可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。X(k)的幅度谱为

X(k)?2XR(k)?XI2(k)，其中下标R和I分别表示取实部、虚部的运算。X(k)的相位谱为

?(k)?arctanXI(k)。

XR(k)2?离散傅里叶反变换（IDFT）定义为

jnk1N?1x(n)??X(k)eNNn?0(n?0,1,...,N?1)。

2.2关于FFT的相关知识

快速傅里叶变换（FFT）是DFT的快速算法，并不是一个新的映射。FFT利用了e?j2?nN函数的周期性

和对称性以及一些特殊值来减少DFT的运算量，可使DFT的运算量下降几个数量级，从而使数字信号处

word文档 可自由复制编辑

理的速度大大提高。

若信号是连续信号，用FFT进行谱分析时，首先必须对信号进行采样，使之变成离散信号，然后就可以用FFT来对连续信号进行谱分析。为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器，且抗混叠滤波器的截止频率不得高于与采样频率的一半。

比较DFT和IDFT的定义，两者的区别仅在于指数因子的指数部分的符号差异和幅度尺度变换，因此可用FFT算法来计算IDFT。

三、实验内容与相关分析（共6道）

说明：为了便于老师查看，现将各题的内容写在这里——

题目按照3.1、3.2、...、3.6排列。每道题包含如下内容：题干、解答（思路、M文件源代码、命令窗口中的运行及其结果）、分析。其中“命令窗口中的运行及其结果”按照小题顺序排列，各小题包含命令与结果（图形或者序列）。

Ω

3.1 求有限长离散时间信号x(n)的离散时间傅里叶变换（DTFT）X(ej)并绘图。 ．．（1）已知x(n)???1?2?n?2n；（2）已知x(n)?2?0其他0?n?10。

【解答】

思路：这是DTFT的变换，按照定义编写DTFT的M文件即可。考虑到自变量Ω是连续的，为了方便计算机计算，计算时只取三个周期[-2π,4π]中均匀的1000个点用于绘图。

理论计算的各序列DTFT表达式，请见本题的分析。 M文件源代码（我的Matlab源文件不支持中文注释，抱歉）： function DTFT(n1,n2,x)

%This is a DTFT function for my experiment of Signal Processing & Analysis. w=0:2*pi/1000:2*pi;Tfine the bracket of omega for plotting. X=zeros(size(w));Tfine the initial values of X. for i=n1:n2

X=X+x(i-n1+1)*exp((-1)*j*w*i);%It is the definition of DTFT. end

Amp=abs(X);?quire the amplification.

Phs=angle(X);?quire 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 (radian)');hold off; %Plot phase angle on the right. end

命令窗口中的运行及其结果（理论计算的各序列DTFT表达式，请见本题的分析）： 第（1）小题

>> n=(-2:2); >> x=1.^n;

>> DTFT(-2,2,x);

word文档 可自由复制编辑

54.544323.5Phase Angle (radian)

Amplification32.521.510-1-210.50-3-4-50510-50510??第（2）小题

>> n=(0:10); >> x=2.^n;

>> DTFT(0,10,x);

2200200018001600图3.1.1在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）

432Phase Angle (radian)-5051010-1-2-3-4Amplification140012001000800600-50510??

图3.1.2在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）

【分析】

对于第（1）小题，由于序列x(n)只在有限区间（-2，-1，-，1，2）上为1，所以是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。事实上，我们可计算出它的表达式：

X(?)?n????x(n)e?j?n???n??2?e?j?n21?e?5j?e2j?1?e?5j???X(?)?，可见幅度频谱拥有主极大?j??j?1?e1?e?? word文档 可自由复制编辑

和次极大，两个主极大间有|5-1|=4个极小，即有3个次级大。而对于它的相位频谱，则是周期性地在-π、0、π之间震荡。

对于第（2）小题，由于是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。而它的表达式：X(?)?n????x(n)e???j?n??2en?010??j?n?1?211e?11j?211??X(?)?，因此主极大之间只有?j??j?1?2e1?2e|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)

%This is a DFT function for my experiment of Signal Processing & Analysis. k=(0:N-1);Tfine variable k for DFT.

X=zeros(size(k));Tfine the initial valves of X. for i=0:N-1

X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i);%It is the definition of DFT. end

Amp=abs(X);?quire the amplification.

Phs=angle(X);?quire 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 (radian)');hold off; %Plot phase angle on the right. end

ii) 基2-FFT函数

function myFFT(N,x)

%This is a base-2 FFT function. lov=(0:N-1); j1=0;

for i=1:N %indexed addressing if i

temp=x(j1+1); x(j1+1)=x(i); x(i)=temp; end

word文档 可自由复制编辑

k=N/2;

while k<=j1 j1=j1-k; k=k/2; end

j1=j1+k; end

digit=0; k=N;

while k>1

digit=digit+1; k=k/2; end

n=N/2;% Now we start the \ for mu=1:digit

dif=2^(mu-1);%Differnce between the indexes of the target variables. idx=1; for i=1:n idx1=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 \ temp=x(idx);

x(idx)=temp+x(idx+dif)*wn;

x(idx+dif)=temp-x(idx+dif)*wn; idx=idx+1; idx2=idx2+1; end

idx=idx1+2*dif; end n=n/2; end

Amp=abs(x);?quire the amplification.

Phs=angle(x);?quire 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

命令窗口中的运行及其结果：

word文档 可自由复制编辑

DFT：

>> x=[8,7,9,5,1,7,9,5]; >> DFT(8,x);

60350230Phase Angle (radian)024k68401Amplification020-110-20-3024k68图3.2.1 DFT的幅度谱和相位谱（弧度制）

FFT： >> x=[8,7,9,5,1,7,9,5]; >> myFFT(8,x);

603

50240FFT Phase Angle (radian)024FFT k681FFT Amplification30020-110-20-302 图3.2.2 FFT算法的幅度谱和相位谱（弧度制）

【分析】

DFT是离散信号、离散频谱之间的映射。在这里我们可以看到序列的频谱也被离散化。事实上，我们可以循着DFT构造的方法验证这个频谱：

首先，将序列做N=8周期延拓，成为离散周期信号。然后利用DFS计算得到延拓后的频谱：

4FFT k68X(?)??x(n)en?0N?1?jk2?n8??x(n)en?07?jkn4??51,7,?9?4i,7,3,6,9?4i,7n?0,1,...,7，从而取DFS???以上作N?8延拓 word文档 可自由复制编辑

的主值区间得到DFT，与图一致。因此计算正确。

而对于FFT，我们可以看到它给出和DFT一样的结果，说明了FFT算法就是DFT的一个等价形式。不过，由于序列不够长，FFT在计算速度上的优越性尚未凸显。

3.3已知连续时间信号x(t)=3cos8πt, X(ω)=3?[?(??8?)??(??8?)]，该信号从t=0开始以采样周期Ts=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(N,Ts,g)

%This is a function applied to Problem 3 & 4.

%N: number of sampling points; 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 our habit. Thus I made a adaptation in codes.

%It is a sampling process with(g=1)/without(g=0) noise. end

xs=sample;

plot(n-1,sample,'.',’MarkerSize’,18);xlabel('n');ylabel('value');hold off; % Must use (n-1), because in Matlab, index starts at 1. end

ii)DFT和基2-FFT函数的代码，请见第3.2节。不需再新编一个。

命令窗口中的运行及其结果： 12点采样：

>> xs=sampJune3(12,0.1,0);%末尾的0表示无噪声。 >> DFT(12,xs); >> myFFT(12,xs);

word文档 可自由复制编辑

word文档

321uela0v-1-2-3024681012n

图3.3.1 进行12点采样得到的序列

202.5182161.5141n)aindo12a0.5ita(r celifgi10lpn0mA A8eas-0.5hP6-14-1.52-20024681012-2.5024681012kk图3.3.2 DFT的幅度谱和相位谱（弧度制），出现了泄漏

可自由复制编辑

2018162.521.51210864200246FFT k81012FFT Phase Angle (radian)1410.50-0.5-1-1.5-2-2.50246FFT k81012FFT Amplification

20点采样：

>> xs=sampJune3(20,0.1,0);%末尾的0表示无噪声。 >> DFT(20,xs); >> myFFT(20,xs);

3图3.3.3 FFT的幅度谱和相位谱（弧度制）。出现了频谱泄漏。

21value0-1-2-30246810n1214161820图3.3.4 进行20点采样得到的序列

word文档 可自由复制编辑

35430325220Phase Angle (radian)0510k1520Amplification115010-15-20-30510k1520

35图3.3.5 DFT的幅度谱和相位谱（弧度制）。频谱无泄漏。 43230FFT Phase Angle (radian)0510FFT k152025FFT Amplification10-1-2-3-4201510500510FFT k1520

图3.3.6 FFT的幅度谱和相位谱（弧度制）。频谱无泄漏。 28点采样：

>> xs=sampJune3(28,0.1,0);%末尾的0表示无噪声。 >> DFT(28,xs); >> myFFT(28,xs);

word文档 可自由复制编辑

本文来源：https://www.bwwdw.com/article/4u15.html