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

实验报告

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

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

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

一、实验目的

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

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

二、实验原理

2.1关于DFT 的相关知识

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

n

j n j e

n x e X Ω-∞

-∞

=Ω

∑=

)()(，

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

n j N n j e n x e X Ω--=Ω

∑=1

)()(，

x (n )的离散傅里叶变换（DFT ）表达式为

)1,...,1,0()()(21

-==--=∑N k e

n x k X nk N

j

N n π，

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

)()()(22k X k X k X I R +=

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

)

()

(arctan

)(k X k X k R I =?。

离散傅里叶反变换（IDFT ）定义为

)1,...,1,0()(1)(21

-==∑-=N n e k X N n x nk N

j N n π

。

2.2关于FFT 的相关知识

快速傅里叶变换（FFT ）是DFT 的快速算法，并不是一个新的映射。FFT 利用了n N

j e

π2-函数的周期性

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

装

订

线

理的速度大大提高。

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

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

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

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

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

3.1 求有限长离散时间信号x (n )的离散时间．．傅里叶变换（DTFT ）X (e j Ω

)并绘图。 （1）已知??

?≤≤-=其他

0221)(n n x ；（2）已知1002

)(≤≤=n n x n

。

【解答】

思路：这是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;%Define the bracket of omega for plotting. X=zeros(size(w));%Define 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);%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 (radian)');hold off ; %Plot phase angle on the right. end

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

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

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

-505100

0.511.5

22.533.5

4

4.55ΩA m p l i f i c a t i o

n

-50510

-4

-3-2

-1012

34Ω

P h a s e A n g l e (r a d i a n )

第（2）小题

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

>> DTFT(0,10,x);

-5

5

10

600

8001000120014001600

180020002200ΩA m p l i f i c a t i o

n

-50510

-4

-3-2-1012

34Ω

P h a s e A n g l e (r a d i a n )

【分析】

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

()

Ω

-Ω-Ω-Ω

-Ω-=Ω-∞

+-∞

=Ω---=Ω?--==

=

Ω∑∑j j j j j n n

j n n j X n x X e

1e 1)(e 1e 1e e e )()(5522

2

，可见幅度频谱拥有主极大

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

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

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

对于第（2）小题，由于是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。而它的表达式：(

)

Ω

-Ω-Ω-=Ω-+∞

-∞

=Ω--≈Ω?--===

Ω∑∑j j j n n

j n n

j X n x X e

212)(e 21e 212e

e

)()(11111110

，因此主极大之间只有|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);%Define variable k for DFT.

X=zeros(size(k));%Define 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);%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 (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

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 "butterfly-shaped" process.

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 "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;

end

idx=idx1+2*dif;

end

n=n/2;

end

Amp=abs(x);%Acquire the amplification.

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

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

DFT ：

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

2

46

8

0102030

405060

k

A m p l i f i c a t i o n

02

468

-3

-2

-1

1

2

3

k

P h a s e A n g l e (r a d i a n )

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

2

46

8

010203040

5060

FFT k

F F T A m p l i f i c a t i o n

02

468

-3

-2

-1

1

2

3

FFT k

F F T P h a s e A n g l e (r a d i a n )

【分析】

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

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

??

?==++-===Ω∑∑=--=-延拓

以上作87

,...,1,07,49,6,3,7,49,7,51e

)(e

)()(7

4

1

8

2N n i i n x n x X n n

jk N n n jk

π

π，从而取DFS 图3.2.1 DFT 的幅度谱和相位谱（弧度制）

图3.2.2 FFT 算法的幅度谱和相位谱（弧度制）

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

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

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

024

681012

-3

-2

-1

1

2

3

n

v a l u e

2

4

68

10

12

024********

161820k

A m p l i f i c a t i o n

2

4

68

10

12

-2.5

-2-1.5-1-0.500.51

1.52

2.5k

P h a s e A n g l e (r a d i a n )

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

图3.3.2 DFT 的幅度谱和相位谱（弧度制），出现了泄漏

2

4

68

10

12

0246810121416

1820FFT k

F F T A m p l i f i c a t i o n

2

4

68

10

12

-2.5

-2-1.5-1-0.500.511.5

22.5FFT k

F F T P h a s e A n g l e (r a d i a n )

20点采样：

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

02468

101214161820

-3

-2

-1

1

2

3

n

v a l u e

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

图3.3.4 进行20点采样得到的序列

5

1015

20

05101520

253035k

A m p l i f i c a t i o n

05

101520

-3

-2

-1

1

2

3

4

k

P h a s e A n g l e (r a d i a n )

5

1015

20

5

10

15

20

25

30

35

FFT k

F F T A m p l i f i c a t i o n

05

101520

-4

-3-2-1012

34

FFT k

F F T P h a s e A n g l e (r a d i a n )

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

图3.3.5 DFT 的幅度谱和相位谱（弧度制）。频谱无泄漏。 图3.3.6 FFT 的幅度谱和相位谱（弧度制）。频谱无泄漏。

0510

15202530

-3

-2

-1

1

2

3

n

v a l u e

10

20

30

510152025303540k

A m p l i f i c a t i o n

010

2030

-3

-2

-1

1

2

3

4k

P h a s e A n g l e (r a d i a n )

图3.3.7 进行28点采样得到的序列

图3.3.8 DFT 的幅度谱和相位谱（弧度制）。再次出现频谱泄漏。

10

20

30

510152025303540FFT k

F F T A m p l i f i c a t i o n

010

2030

-3

-2

-1

1

2

3

4

FFT k

F F T P h a s e A n g l e (r a d i a n )

【分析】

分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏之间的关系。现在与原信号频谱)(ωX 比较后可以得出如下结论：

原信号的频谱是)]8()8([3)(πωδπωδπω++-=X ，在±8π上各有一强度为3π的谱线，在其余频率上为0。

可见原信号被0.1 s 采样周期的采样信号离散化之后，谱线以20π为周期重复，并且只在（20k ±8）π (k 为整数）处非0。那么，在20点DFT （采样时间原信号周期的整数倍）中，只有第8根、第12根谱线非0。而在12点、28点DFT 中，由于采样时间不是原信号周期的整数倍，谱线将向两边泄漏。

不过，对比12点采样和28点采样，我们还可以发现，28点采样频谱的主谱线高度是次谱线高度的4倍，儿12点采样频谱的主谱线高度是次谱线高度的3倍。可见，在无法保证采样时间是信号周期整数倍的情况下，增加采样时间有助于减轻频谱泄漏的程度。

图3.3.9 FFT 的幅度谱和相位谱（弧度制）。再次出现频谱泄漏。

图3.3.10 原信号的频谱（由两个冲激函数组成）

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)

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

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

2

4

68

10

12

-5-4-3-2-101

2345n

v a l u e

图3.4.1 进行12点采样得到的含噪声的序列

024

6810120

24

68101214

16

18k

A m p l i f i c a t i o n

024

681012

-3

-2

-1

1

2

3

k

P h a s e A n g l e (r a d i a n )

2

4

68

10

12

24

6810121416

18FFT k

F F T A m p l i f i c a t i o n

024

681012

-3

-2

-1

1

2

3

FFT k

F F T P h a s e A n g l e (r a d i a n )

20点采样：

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

图3.4.2 含噪声序列DFT 的幅度谱和相位谱（弧度制）。

图3.4.3 含噪声FFT 的幅度谱和相位谱（弧度制）。

02468

101214161820

-4

-3-2-101

234n

v a l u e

5

1015

20

051015

202530

k

A m p l i f i c a t i o n

05

101520

-3

-2

-1

1

2

3

k

P h a s e A n g l e (r a d i a n )

图3.4.4 进行20点采样得到的含噪声序列

图3.4.5 含噪声DFT 的幅度谱和相位谱（弧度制）。

5

1015

20

05101520

2530FFT k

F F T A m p l i f i c a t i o n

05

101520

-3

-2

-1

1

2

3

FFT k

F F T P h a s e A n g l e (r a d i a n )

28点采样：

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

0510

15202530

-4

-3-2-101

234n

v a l u e

图3.4.6 含噪声FFT 的幅度谱和相位谱（弧度制）。

图3.4.7 进行28点采样得到的含噪声序列

10

20

30

510152025303540k

A m p l i f i c a t i o n

010

2030

-4

-3

-2

-1

1

2

3

k

P h a s e A n g l e (r a d i a n )

10

20

30

510152025303540FFT k

F F T A m p l i f i c a t i o n

010

2030

-3

-2

-1

1

2

3

4

FFT k

F F T P h a s e A n g l e (r a d i a n )

【分析】

依然分别取12点、20点、28点采样。仍然与原信号的频谱)]8()8([3)(πωδπωδπω++-=X （图

图3.4.8 含噪声DFT 的幅度谱和相位谱（弧度制）。

图3.4.9 含噪声FFT 的幅度谱和相位谱（弧度制）。

3.3.10）比较，可以得到结论：

由于叠加了噪声，所以频谱都受到了一定的干扰。由于白噪声在各个频率的功率相等，因此频谱上各处的干扰也是均匀随机的。

不过，通过对比我们可以发现，20点采样（无噪声时不发生泄漏的采样方法）在存在噪声时，仍然可以明显区分出原信号的谱线。第二好的是28点采样，因为采样时间较长，即使存在频谱泄漏也能较好地区分原信号的谱线。而最差的是12点采样，由于噪声的存在和严重的频谱泄漏，它的次谱线与主谱线的高度相差不大，使原信号不明显。

3.5已知序列)3()2(2)1(3)(4)(-+-+-+=n n n n n x δδδδ，X (k )是x (n )的6点DFT ，设k N

j

k N

e

W π2=-。

（1）若有限长序列y (n )的6点DFT 是)()(46k X W k Y k -=，求y (n )。 （2）若有限长序列w (n )的6点DFT W (k )是)(k X 的实部，求w (n )。 （3）若有限长序列q (n )的3点DFT 是)2()(k X k Q =，k=0,1,2，求q (n )。

【解答】

思路：这是对DFT 进行变换后求IDFT 。考虑到IDFT 和DFT 定义的对称性，可以在DFT 的基础上略加调整既可用于计算。

首先，∵)3()2(2)1(3)(4)(-+-+-+=n n n n n x δδδδ，

∴它的6点采样是序列是}0,0,1,2,3,4{)(=n x 。值得指出的是，在Matlab 中，数组的序号是从1开始的（而在信号分析中习惯从0开始），不过我在上面编程时已考虑到这一情况，具体可见实验报告最后的“附录”。 首先生成x (n )的6点DFT ，再按照各小题分别转换，最后求相应的IDFT 。 M 文件源代码：

i) 输出x (n )的6点DFT 的函数：

function X = ExportDFT(N,x)

%This is a DFT function that 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-1

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

ii)算第（1）小题的Y (k )的函数： function Y = Convertor1(X)

%This is a mathematical convertor for the subproblem (1). for k=1:6

Y(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 end

iii)算第（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.

end

iv)算第（3）小题的Q(k)的函数：

function Q = Convertor3(X)

%This is a mathematical convertor for the subproblem (3).

Q=zeros(3);% Detailed explanation goes here

for tmp=1:3

Q(tmp)=X(2*tmp);

end

end

v)输出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));%Define the initial valves of x.

for k=0:N-1

x=x+X(k+1)*exp(j*2*k*pi/N*n);

end

x=x/N;

a=real(x);%We MUST use real(x), though we may ALREADY know x is real.

%It's 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,0];

>> X=ExportDFT(6,x);

第（1）小题

>> Y=Convertor1(X);

>> y=ExportIDFT(6,Y)

y =

Columns 1 through 4

0.0000 + 0.0000i 0.0000 + 0.0000i 4.0000 + 0.0000i 3.0000 + 0.0000i

%虚部都是0，说明是实数 Columns 5 through 6

2.0000 + 0.0000i 1.0000 - 0.0000i %虚部都是0，说明是实数

%事实上，在Matlab中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。

012

345

-0.5

00.511.522.533.54n

S o l u t i o n

第（2）小题

>> W=Convertor2(X); >> w=ExportIDFT(6,W) w =

Columns 1 through 4

4.0000 + 0.0000i 1.5000 + 0.0000i 1.0000 + 0.0000i 1.0000 + 0.0000i

%虚部都是0，说明是实数

Columns 5 through 6

1.0000 + 0.0000i 1.5000 + 0.0000i %虚部都是0，说明是实数；

%事实上，在Matlab 中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。

012

345

0.511.522.533.54n

S o l u t i o n

第（3）小题

>> Q=Convertor3(X); >> q=ExportIDFT(6,Q) q =

Columns 1 through 4

1.5000 - 0.0000i -0.1667 - 0.2887i 0.7500 - 1.2990i 0.8333 - 0.0000i Columns 5 through 6

-0.5000 - 0.8660i 1.0833 - 1.8764i

这里的答案都是幅值、相位均非0的复数，而教材（实验指导第109页）并未要求作图，这里略去。

答案：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 (n )。