数字信号处理 实验六 离散傅立叶变换及其快速算法

实验六 离散傅立叶变换及其快速算法
一、 实验目的：掌握快速傅立叶变换的应用方法；掌握离散余弦变换的应用方法；
掌握Z 变换的应用方法；了解Chip z 变换的基本概念；掌握Hilbeit 变换的初步应用；了解倒谱变换的基本概念。

二、 实验仪器：电脑一台，MA TLAB6.5或更高级版本软件一套。

三、 实验内容：
（一） 快速傅立叶变换（FFT ）
在信号处理中，DFT 的计算具有举足轻重的地位，，信号的相关、滤波、谱估计等都要
通过DFT 来实现。

然而，当N 很大的时候，求一个N 点的DFT 要完成N N ⨯次复数乘法和)1(-N N 次复数加法，其计算量相当大。

1965年J.W .Cooley 和J.W .Tukey 巧妙地利用N W 因子的周期性和对称性，构造了一个DFT 快速算法，即快速傅立叶变换(FFT)。

概念
通过前面的知识，已经知道有限列长为N 的序列)(n x 的DFT 变换为
nk
n
N n W n x k X ∑
-==
1
)()( 12,1,0-=N k
其逆变换为
∑
-=-=
1
)(1)(N k nk
N
W k X N
n x 1,1,0-=N n
由于MA TLAB 软件本身的特点，序列或向量元素下标从1开始记录，而不是从0开始。

因此，上述两式在MA TLAB 中相应的表达式为
nk
n
N n W n x k X ∑
-==
1
)()( 12,1,0-=N k
∑
-==
1
1
)(1)(N k nk
N
W k X N
n x 1,2,1-=N n
而下面所讨论使用的快速傅立叶变换)(FFT 并不是与DFT 不同的另外一种变换，而是为减少DFT 计算次数的一种快速有效的算法。

这种快速算法，主要是利用了nk
N W 下面两个
特性使长序列的DFT 分解为更小点数的DFT 所实现的。

（1） 利用nk
N
W 的对称性使DFT 运算中有些项合并
*
--==)()
(kn
N
kn
N
n N k N
W W W
（2） 利用nk
N
W 的周期性和对称性使长序列的DFT 分解为更小点数的DFT
)()
()
(n
N k N
N n k N
nk
N
W W W ++==
快速傅立叶变换)(FFT 算法正是基于这一基本思想而发展起来的。

快速傅立叶变换算法形式很多，但是基本上可以分为两大类，即按时间抽取（Decimation －In －Time ，简称
DIT ）法和按频率抽取（Decimation －In －Frequency ）法。

在这里，以时间抽取)(DIT 的FFT 算法（库利－图算法）为例，简单说明一下FFT 算法的算法原理。

为了讨论方便，设γ2=N ，其中γ为整数。

如果不满足这个条件，可以认为得加上若干零点来达到。

由DFT 的定义知
nk
n
N n W n x k X ∑
-==
1
)()( 12,1,0-=N k
其中)(n x 是列长为)11,0(-=N n N 的输入序列，把它按n 的奇偶分成两个序列
⎩⎨
⎧=+=)
()12()()2(21r x r x r x r x 12
,1,0-=N r
又由于2
2
222N
N
j
N
j
N
W e
e W
===--ππ，
则∑
∑
-=-=+=+
=
1
211
)()()()()(N n n k
N nk
N
nk
n
N n n k X W k X W n x W n x k X 为奇数
为偶数
上式表明了一个N 点的DFT 可以被分解为两个N/2点的DFT 。

同时，这两个N/2点的DFT 按照上式又可以合成为一个N 点的DFT 。

为了要用点数为N/2点的X 1(k)、X 2(k)来表达N 点的X(k)值还必须要用W 系数的周期性，即)
2(2
2
N k r N
rk N W W
+
=
这样可得
∑∑==+
=
+2
2
2
1
)
2(2
11)()()2
(
N
r N
r rk
N N k r N
W r x
W r x k N X
即 X 1（
)()2
1k X k N =+
同理可得 X 2（)()2
2k X k N =+
另外再加上W k
N 的对称性
k
N k N N
N k N
N W W W W -=∙=+2)
2
(
就可以将X(k)的表达式分为前后两个部分：
前半部分 )()()(21k X W k X k X k
N += 12
,1,0-=N r
后半部分 )2
(
)2
(
)2
(
2)
2
(
1k N X W k N X k N X k N N
+++=++
)()(21k X W k X k
N -= 12
,1,0-=N r
由以上分析可见，只要求出区间⎥⎦
⎤
⎢⎣
⎡-12,
0N
内各个整数k 值所对应的X 1(k)、
X 2(k)的值，即可求出[]1,0-N 区间内的全部X （k ）值，这一点恰恰是FFT 能大量节省计算的关键所在。

（二） 函数应用
MA TLAB 为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有Fft 、Ifft 、Fft2 、Ifft2, Fftn 、ifftn 和Fftshift 、Ifftshift 等。

当所处理的数据的长度为2的幂次
时，采用基-2算法进行计算，计算速度会显著增加。

所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。

1.Fft和Ifft函数
调用方式
（1）Y=fft(X)
参数说明
·如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换；
·如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；
维数组，则是对第一个非单元素的维进行离散傅立叶变换；
·如果X是（N)
D
（2）Y＝fft(X,N)
参数说明
N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。

（3）（3）Y＝fft(X,[],dim) 或Y＝fft(X,N,dim)
参数说明
·在参数dim指定的维上进行离散傅立叶变换；
·当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。

函数Ifft的参数应用与函数Fft完全相同。

应用说明
【实例6-1】fft的应用
X=[2 1 2 8];
Y=fft(X)
运行结果
Y＝
13.0000 0+7.0000i -5.0000 0-7.0000
【实例6-2】fft(X,N,dim)的应用
A=[2 5 7 8;
1 4 0 5;
3 8 5 1;
9 1 2 7];
Z=fft(A,[],1)
【实例6-3】fft在信号分析中的应用
使用频率分析方法从受噪声污染的信号x(t)中鉴别出有用的信号。

程序：
t=0:0.001:1; %采样周期为0.001s,即采样频率为1000Hz;
%产生受噪声污染的正县正弦波信号；
x=sin(2*pi*100*t)+sin(2*pi*200*t)+rand(size(t));
subplot(2,1,1)
plot(x(1:50)); %画出时域内的信号；
Y=fft(x,512); %对X进行512点的傅立叶变换；
f=1000*(0:256)/512; %设置频率轴（横轴）坐标，1000为采样频率；subplot(2,1,2)
plot(f,Y(1:257)); %画出频域内的信号
运行结果
图6-1 时域信号和频域信号的比较
由上面的图6-1可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。

但是通过对x(t)作傅立叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100Hz 和200Hz的两个频率分量。

2.Fft2和Ifft函数
调用方式
（1）Y=fft2(X)
参数说明
·如果X是向量，则此傅立叶变换即变成一维傅立叶变换fft；
·如果X是矩阵，则是计算该矩阵的二维快速傅立叶变换；
·数据二维傅立叶变换fft2（X）相当于fft(fft(X)'')，即先对X的列做一维傅立叶变换，然后再对变换结果的行做一维傅立叶变换。

（2）Y＝fft2(X,M,N)
通过对X进行补零或截断，使得成为（M*N）的矩阵。

函数Ifft2的参数应用与函数Fft2完全相同。

Fftn.、Ifftn是对数据进行多维快速傅立叶变换，其应用与Fft2、Ifft2类似；在此，不再一
一赘述。

【实例6-4】fft2、ifft2的应用
A=[2 5 7 8 9;
1 3 7 5 0;
2 6 1 4 9;
8 1 5 2 6];
Y=fft2(A);
B=ifft2(Y);
运行结果：
3.Fftshift和Ifftshift函数
调用方式
Z=fftshift(Y)
此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0处得值）移到频谱的中间位置。

参数说明
·如果X是向量，则变换Y的左右两边；
·如果X是矩阵，则交换Y的一、三象限和二、四象限；
·如果Y是多维数组，则在数组的每一维交换其“半空间”。

函数Iffshift的参数应用与函数Fftshift完全相同。

【实例6-5】fftshift的应用
X=rand(5,4);
y=fft(X);
z=fftshift(y);%只将傅立叶变换结果y 中的直流成分移到频谱的中间位置. 运行结果： y=
3.2250 2.5277 1.4820 1.6314 0.3294+0.2368i 0.0768+0.3092i 0.6453+0.4519i -0.7240-0.4116i -0.2867-0.6435i 0.5657+0.4661i -0.5515+0.2297i -0.0573-0.0881i -0.2867+0.6435i 0.5657-0.4661i -0.5515-0.2297i -0.0573+0.0881i 0.3294-0.2368i 0.0768-0.3092i 0.6453-0.4519i -0.7240+0.4116i Z=
-0.5515-0.2297i -0.0573+0.0881i -0.2867+0.6435i 0.5657-0.4661i 0.6453-0.4519i -0.7240+0.4116i 0.3294-0.2368i 0.0768-0.3092i 1.4820 1.6314 3.2250 2.5277 0.6453+0.4519i -0.7240-0.4116i 0.3294+0.2368i 0.0768+0.3092i -0.5515+0.2297i -0.0573-0.0881i -0.2867-0.6435i 0.5657+0.4661i （三） 离散傅立叶变换（IFFT ） 逆变换IFFT 为
nk
N
N k W n X N
n x --=∑
=
1
)(1)( 12
,1,0-=N r
通过下列修改，就可以用FFT 算法来实现逆FFT 运算：
·增加一个归一化因子
N
1；
·将用nk N W 其复数共轭nk
N W -代替
但是由因为，
*
**
1
*)]}([{1])([1)(k X FFT N
W n X N
n x nk
N
N k =
=
∑-=
所以，求X （k ）的逆FFT 可以分为以下3个步骤： （1）取X （k ）的共扼，得到)(*k X ；
（2）求)(*k X 的FFT ，得N )(*k x ；
（3）取)(*n x 的共扼，并除以N ，即得到x(n)。

采用这种方法可以直接用FFT 程序来计算逆FFT 。

有关IFFT 的具体应用，与FFT 一一对应，在此不再赘述。

（四） 编程练习 1． 试
用
Mablab
求其有限长序列
)
100()8.0()(1≤≤=n n x n
与
)180()6.0()(2≤≤=n n x n
的圆周卷积，
（N=20），并画出其结果图。