详解FFT快速傅里叶变换FFT.docx

详解FFT快速傅里叶变换FFT.docx

《详解FFT快速傅里叶变换FFT.docx》由会员分享，可在线阅读，更多相关《详解FFT快速傅里叶变换FFT.docx（27页珍藏版）》请在冰点文库上搜索。

详解FFT快速傅里叶变换FFT

第四章快速傅里叶变换

有限长序列可以通过离散傅里叶变换（DFT）将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换（FFT）.1965年，Cooley和Tukey提出了计算离散傅里叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。

从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。

根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基２DIT和基２DIF。

FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。

快速傅里叶变换（FFT）是计算离散傅里叶变换（DFT）的快速算法。

DFT的定义式为

kn

N−1

X（k）=∑x（n）WN

RN（k）

n=0

N

在所有复指数值Wkn的值全部已算好的情况下，要计算一个X（k）需要N

次复数乘法和N－1次复数加法。

算出全部N点X（k）共需N2次复数乘法和N（N−1）次复数加法。

即计算量是与N2成正比的。

FFT的基本思想：

将大点数的DFT分解为若干个小点数DFT的组合，从而减少运算量。

WN因子具有以下两个特性，可使DFT运算量尽量分解为小点数的DFT

运算：

（1）周期性：

（k+N）n

N

=Wkn

=W（n+N）k

（2）对称性：

W（k+N/2）=−Wk

NN

利用这两个性质，可以使DFT运算中有些项合并，以减少乘法次数。

例子：

求当N＝4时，X

（2）的值

3

∑44444

X

（2）=

n=0

x（n）W2n=x（0）W0+x

（1）W2+x

（2）W4+x（3）W6

=[x（0）+x

（2）]W0+[x

（1）+x（3）]W2

（周期性）

4

4

＝[x（0）+x

（2）]－[x

（1）+x（3）]W0

4

（对称性）

通过合并，使乘法次数由4次减少到1次，运算量减少。

FFT的算法形式有很多种，但基本上可以分为两大类：

按时间抽取

（DIT）和按频率抽取（DIF）。

4.1按时间抽取（DIT）的FTT

为了将大点数的DFT分解为小点数的DFT运算，要求序列的长度N为复合数，最常用的是N=2M的情况（M为正整数）。

该情况下的变换称为基2FFT。

下面讨论基2情况的算法。

先将序列x（n）按奇偶项分解为两组

⎧x（2r）=x1（r）

⎨

⎩x（2r+1）=x2（r）

将DFT运算也相应分为两组

N

r=0,1,L,2−1

N−1

N

X（k）=DFT[x（n）]=∑x（n）Wkn

n=0

N−1

=∑x（n）Wkn＋

N−1

∑x（n）Wkn

n=0

n为偶数

n=0

n为奇数

N/2−1

∑

（2）

2rk+

N/2−1

∑（2

+1）

（2r+1）k

＝x

r=0

rWN

xrWN

r=0

N/2−1

∑

2rk+

N/2−1

k∑

2rk

＝

r=0

x1（r）WN

WN

k

r=0

x2（r）WN

rk

N/2−1

＝∑x1（r）W

r=0

N/2+WN

N/2−1

∑x2（r）W

r=0

rk

N/2

（因为W

2rk

N

rk

N/2

k

=X1（k）+WNX2（k）

其中X1（k）、X2（k）分别是x1（n）、x2（n）的N/2点的DFT

N/2−1

N/2−1

rkrk

X1（k）＝∑x1（r）WN/2=∑x（2r）WN/2,0≤k≤−1

r=0

r=02

N/2−1

N/2−1

rk

rkN

X2（k）＝∑x2（r）WN/2=∑x（2r+1）WN/2,0≤k≤−1

r=0

r=02

至此，一个N点DFT被分解为两个N/2点的DFT。

上面是否将全部N点的X（k）求解出来了？

分析：

X1（k）和

N

X2（k）只有N/2个点（k=0,1,L,

2

−1），则由

1N2

X（k）=X（k）+WkX（k）只能求出X（k）的前N/2个点的DFT，要求出

全部N点的X（k），需要找出X1（k）、X2（k）和X（k+N/2）的关系，其

N

中k=0,1,L,2−1。

由式子X（k）=X1

N

（k）+WkX

2（k）可得

1N

X（k+N/2）=X（k+N/2）+Wk+N/2

k

X2（k+N/2）化简得

N

X（k+N/2）＝=X1（k）−WNX2（k），k=0,1,L,2−1

这样N点DFT可全部由下式确定出来：

1N2

⎧⎪X（k）=X（k）+WkX（k）

1N2

⎩⎪⎨X（k+N/2）=X（k）−WkX（k）

k=0,1,L

N−1

2

（＊）

上式可用一个专用的碟形符号来表示，这个符号对应一次复乘和两次复加运算。

aa+Wkb

ba−Wkb

-1

N

图蝶形运算符号

2

通过这样的分解以后，每一个N／2点的DFT只需要（N）2=N

次复数乘

24

2

法，两个N/2点的DFT需要2（N）2=N

次复乘，再加上将两个N／2点

22

DFT合并成为N点DFT时有N／2次与W因子相乘，一共需要

2

N+N

22

N2

≈次复乘。

可见，通过这样的分解，运算量节省了近一半。

2

因为N=2M，N/2仍然是偶数，因此可以对两个N/2点的DFT再

分别作进一步的分解，将两个N/2点的DFT分解成两个N/4点的DFT。

例如对x1（r），可以在按其偶数部分及奇数部分进行分解：

⎧x1（2l）=x3（l）

⎨

⎩x1（2l+1）=x4（l）

则的运算可相应分为两组：

N

l=0,1,L,4−1

N/4−1

X1（k）＝∑x1（2l）W

l=0

2lk

N/2

N/4−1

+∑x1（2l+1）W

l=0

（2l+1）k

N/2

N/4−1

∑

3（）/4

k

N/2

N/4−1

∑

4（）

lk

N/4

＝x

l=0

lWNW

k

xlW

l=0

N

=X3（k）+WN/2X4（k）

N/2

将系数统一为以Ｎ为周期，即Wk

k=0,1,L,4−1

N

=W2k，可得

13N

⎧⎪X（k）=X（k）+W2k

⎨

X4（k）

2k

k=0,1,L

N−1

X1（k+N/4）=X3（k）−WN

X4（k）4

同样，对X2（k）也可进行类似的分解。

一直分解下去，最后是２点的

DFT，２点DFT的运算也可用碟形符号来表示。

这样，对于一个N=23=8

的DFT运算，其按时间抽取的分解过程及完整流图如下图所示。

这种方法，由于每一步分解都是按输入序列在时域上的次序是属于偶数

还是奇数来抽取的，故称为“时间抽取法”。

分析上面的流图，N=2M，一共要进行M次分解，构成了从x（n）到

X（k）的M级运算过程。

每一级运算都是由N/2个蝶形运算构成，因此每一级运算都需要N/2次复乘和N次复加，则按时间抽取的M级运算后总共需要

复数乘法次数：

mF

=N⋅M

2

=NlogN

22

复数加法次数：

aF

=N⋅M=Nlog2N

根据上面的流图，分析ＦＦＴ算法的两个特点，它们对ＦＦＴ的软硬件

构成产生很大的影响。

（１）原位运算也称为同址运算，当数据输入到存储器中以后，每一级运算的结果仍然存储在原来的存储器中，直到最后输出，中间无需其它的存储器。

根据运算流图分析原位运算是如何进行的。

原位运算的结构可以节省存储单元，降低设备成本。

（２）变址分析运算流图中的输入输出序列的顺序，输出按顺序，输入是“码位倒置”

的顺序。

见图。

自然顺序

二进制表示

码位倒置

码位倒置顺序

0

000

000

0

1

001

100

4

2

010

010

2

3

011

110

6

4

100

001

1

5

101

101

5

6

110

011

3

7

111

111

7

码位倒置顺序

01234567

X（0）X（4）X

（2）X（6）X

（1）X（5）X（3）X（7）码位倒置的变址处理

在实际运算中，直接将输入数据x（n）按码位倒置的顺序排好输入很不

方便，一般总是先按自然顺序输入存储单元，然后通过变址运算将自然顺序

的存储换成码位倒置顺序的存储，这样就可以进行FFT的原位运算。

变质的功能如图所示。

用软件实现是通用采用雷德（Rader）算法，算出I的倒序Ｊ以后立即将输入数据X（I）和X（J）对换。

尽管变址运算所占运算量的比例很小，但对某些高要求的应用（尤其在实时信号处理中），也可设法用适当的电路结构直接实现变址。

例如单片数字信号处理器TMS320C25就有专用于FFT的二进制码变址模式。

4.２按频率抽取（DIF）的FTT

除时间抽取法外，另外一种普遍使用的FFT结构是频率抽取法。

频率抽取法将输入序列不是按奇、偶分组，而是将Ｎ点DFT写成前后两部分：

N−1

N

X（k）=DFT[x（n）]=∑x（n）Wkn

n=0

（N/2）−1

N−1

=∑x（n）Wkn＋∑x（n）Wkn

n=0

NN

n=N/2

N/2−1

∑（）

nk+

N/2−1

∑（+

nk

/2）

（n+N/2）k

＝

n=0

xnWN

xnNWN

n=0

（N/2）k

N/2−1

＝∑[x（n）+WN

x（n+N/2）]WN

n=0

因为WN/2=−1,W（N/2）k

=（−1）k，k为偶数时（−1）k

=1，k为奇数时

NN

（−1）k

=−1，由此可将X（k）分解为偶数组和奇数组：

N/2−1

∑N

X（k）＝

n=0

[x（n）+（−1）kx（n+N/2）]Wnk

N/2−1

∑N

X（2r）＝

n=0

[x（n）+x（n+N/2）]W2nr

N/2−1

∑

N/2

=

n=0

[x（n）+x（n+N/2）]Wnr

N/2−1

∑N

X（2r+1）＝

N/2−1

n=0

[x（n）−x（n+N/2）]W（2r+1）n

∑NN/2

=

N

n=0

[x（n）−x（n+N/2）]WnWnr

⎧x1（n）=x（n）+x（n+N/2）

令⎨

n=0,1,L,N/2−1

⎩x2

（n）=[x（n）−x（n+N/2）]Wn

这两个序列都是N/2点的序列，对应的是两个N/2点的DFT运算：

N/2−1

∑1

N/2

X（2r）＝

n=0

x（n）]Wnr

N/2−1

∑2

N/2

X（2r+1）＝

x

n=0

（n）Wrn

这样，同样是将一个N点的DFT分解为两个N/2点的DFT了。

频率抽选法

对应的碟形运算关系图如下：

aa+b

（a−b）Wn

bN

n

-1WN

对于N=8时频率抽取法的FFT流图如下：

这种分组的办法由于每次都是按输出X（k）在频域的顺序上是属于偶数还是

奇数来分组的，称为频率抽取法。

与前面按时间抽取的方法相比，相同点问题：

如何利用快速算法计算IDFT？

分析IDFT的公式：

x（n）=IDFT[X（k）]=

1N−1

N

比较DFT的公式：

Nk=0

N−1

N

X（k）=DFT[x（n）]=∑x（n）Wnk,k=0,1,L,N−1

n=0

得知可用两种方法来实现IDFT的快速算法：

（１）只要把DFT运算中的每

一个系数Wnk该为W−nk，并且最后再乘以常数1，就可以用时间抽取法

NNN

或频率抽取的FFT算法来直接计算IDFT。

这种方法需要对FFT的程序和参

数稍加改动才能实现。

（２）因为

x（n）=

1N−1

N

1{DFT[X*（k）]},n=0,1,L,N−1，也就

[∑X*（k）Wnk]*=

Nk=0N

是说，可先将X（k）取共轭变换，即将X（k）的虚部乘以－１，就可直接调用

FFT的程序，最后再对运算结果取一次共轭变换并乘以常数1/N即可得到x（n）的值。

这种方法中，FFT运算和IFFT运算都可以共用一个子程序块，在使用通用计算机或用硬件实现时比较方便。

4.1.3混合基FFT算法

以上讨论的是基2的FFT算法，即N=2M的情况，这种情况实际上使用得最多，这种FFT运算，程序简单，效率很高，用起来很方便。

另

外，在实际应用时，有限长序列的长度N到底是多少在很大程度时是由人

为因素确定的，因此，大多数场合人们可以将N选定为N=2M，从而可以直接调用以２为基数的FFT运算程序。

如果长度N不能认为确定，而N的数值又不是以2为基数的整数次

方，一般可有以下两种处理方法：

（１）将x（n）用补零的方法延长，使N增长到最邻近的一个

N=2M数值。

例如，N=30，可以在序列x（n）中补进

x（30）＝x（31）＝０两个零值点，使N=32。

如果计算FFT的目的是为了了解整个频谱，而不是特定频率点，则此法可行。

因为有限长序列补零以后并不影响其频谱

X（ejw），只是频谱的采样点数增加而已。

（２）如果要求特定频率点的频谱，则N不能改变。

如果N

为复合数，则可以用以任意数为基数的FFT算法来计算。

快速傅里叶变换的基本思想就是要将DFT的运算尽量分小。

例如，N=6时，可以按照N=3×２分解，将６点DFT分解为３组２点DFT。

举例：

N=9时的快速算法。

凡是可以利用傅里叶变换来进行分析、综合、变换的地方，都可以利

用FFT算法及运用数字计算技术加以实现。

FFT在数字通信、语音信号处理、图像处理、匹配滤波以及功率谱估计、仿真、系统分析等各个领域都得到了广泛的应用。

但不管FFT在哪里应用，一般都以卷积积分或相关积分的具体处理为依据，或者以用FFT作为连续傅里叶变换的近似为基础。

4.2.1利用FFT求线性卷积—快速卷积

在实际中常常遇到要求两个序列的线性卷积。

如一个信号序列x（n）

通过FIR滤波器时，其输出y（n）应是x（n）与h（n）的卷积：

y（n）=x（n）∗h（n）=

∞

∑x（m）h（n−m）

m=−∞

有限长序列x（n）与h（n）的卷积的结果y（n）也是一个有限长序列。

假设x（n）与

h（n）的长度分别为N1和N2，则y（n）的长度为N=N1+N2-1。

若通过补零使x（n）与h（n）都加长到N点，就可以用圆周卷积来计算线性卷积。

这样得到用FFT运算来求y（n）值（快速卷积）的步骤如下：

（１）对序列x（n）与h（n）补零至长为N，使N≥N1+N2-1，并且

N=2M（M为整数），即

⎨

x（n）=⎧x（n）,n=0,1,L,N1−1

⎩0,

n=N1,N1+1,L,N−1

⎨

h（n）=⎧h（n）,n=0,1,L,N2−1

⎩0,

n=N2,N2+1,L,N−1

（２）用FFT计算x（n）与h（n）的离散傅里叶变换

x（n）⇔（FFT）⇔X（k）

（N点）

h（n）⇔（FFT）⇔H（k）

（N点）

（３）计算Y（k）=X（k）H（k）

（４）用IFFT计算Y（k）的离散傅里叶反变换得：

y（n）=IFFT[Y（k）]（N点）

互相关及自相关的运算已广泛的应用于信号分析与统计分析，应用于连

续时间系统也用于离散事件系统。

用FFT计算相关函数称为快速相关，它与快速卷积完全类似，不同的是一个应用离散相关定理，另一个应用离散卷积定理。

同样都要注意到离散傅里叶变换固有的周期性，也同样用补零的方法来绕过这个障碍。

设两个离散时间信号x（n）与y（n）为已知，离散互相关函数记作Rxy（n），定义为

Rxy（n）=

∞

∑x（m）y（n+m）

m=−∞

如果x（n）与y（n）的序列长度分别为N1和N2，则用FFT求相关的计算步骤

如下：

（１）对序列x（n）与y（n）补零至长为N，使N≥N1+N2-1，并且

N=2M（M为整数），即

⎨

x（n）=⎧x（n）,n=0,1,L,N1−1

⎩0,

n=N1,N1+1,L,N−1

⎨

y（n）=⎧y（n）,n=0,1,L,N2−1

⎩0,

n=N2,N2+1,L,N−1

（２）用FFT计算x（n）与y（n）的离散傅里叶变换

x（n）⇔（FFT）⇔X（k）

（N点）

y（n）⇔（FFT）⇔Y（k）

（N点）

（３）将X（k）的虚部Im[X（k）]改变符号，求得其共轭X*（k）

（４）计算Rxy（k）=X*（k）Y（k）

（５）用IFFT求得相关序列Rxy（n）

Rxy（n）＝IFFT[Rxy（k）]（N点）

如果x（n）＝y（n），则求得的是自相关序列Rxx（n）

4.2.3Chirp－Z变换

采用FFT算法可以很快的计算出全部DFT值，也即Z变换在单位圆上的全部等间隔采样值。

但是，很多场合下，并非整个单位圆上的频谱都是很有意义的，例如对于窄带信号过程，往往只需要对信号所在的一段频带进行分析，这是，希望采样能密集在这段频带内，而对频带以外的部分，则可以完全不管。

另外，有时也希望采样能不局限于单位圆上，例如，语音信号处理中，往往需要知道其Z变换的极点所在频率，如果极点位置离单位圆较远，，则其单位圆上的频谱就很平滑，如图所示，这是很难从中识别出极点所在的频率。

要是采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则所得的结果将会在极点所在频率上出现明显的尖峰，从而可以较准确的测定极点频率。

螺线采样就是一种适用于这种需要的变换，并且可以采用FFT来快

速计算，这种变换也称为Chirp-Z变换，它是沿Z平面上的一段螺线作等分角的采样，这些采样点可表达为

k

z=AW−k,k=0,1,L,M−1

其中M为采样点的总数，A为起始点位置，这个位置可以进一步用它的半径

0

A0及相角θ0来表示

参数W可表示为

A=Aejθ0

0

W=WejΦ0

其中W0为螺线的伸展率，W0>１时螺线内缩（反时针方向）；W0<１时螺

线外伸。

Φ0为螺线上采样点之间的等分角。

螺线采样点在Z平面上的分布可表示为下图。

下面分析这些点上采样值计算的特点。

假定x（n）是长度为N的有限

长信号序列，则其Z变换在采样点上的值为