快速傅里叶变换FFT算法及其应用当文网提供

1、快速傅里叶变换FFT算法及其应用本文较为系统地阐述了快速傅里叶变换的算法原理及其在数字信号处理等 工程技术中的应用。根据抽取方法的不同，一维基2 FFT算法分为两种：频域抽 取的FFT算法和时频域抽取的FFT算法。第1节阐述了这两种FFT算法的原理。 第2节给出了两种算法的编程思想和步骤。 第3节阐述了一维非基2 FFT的两种 算法：Cooley-tukey FFT算法和素因子算法(Prime Factor Algorithm)的思想原理， 给出了在把一维非基2 DFT的多层分解式转化为二层分解的过程中，如何综合 运用这两种算法以达到总运算次数最少的方案；并以20点DFT为例描述了非基2 FF

2、T算法实现的一般步骤。第4节介绍了一维FFT算法在计算连续时间信号的 傅里叶变换、离散信号的线性卷积、离散信号压缩和滤波等数字信号处理中的典 型应用。第5节把一维FFT变换推广到二维FFT变换，并在一维FFT算法的基 础上，给出了二维FFT算法的原理和实现过程。最后在附录中给出了一维DFT的基2 FFT算法(包括频域抽取的 FFT和IFFT算法、时域抽取的FFT和IFFT 算法)，一维任意非基2 FFT算法，二维DFT的基2 FFT算法以及二维DFT的 任意非基2 FFT算法的详细的Visual C+程序。本文通过各种流程图和表格，较为深入系统地阐述了FFT的算法原理；运用Matlab编程，通

3、过大量生动的实例，图文并茂地列举出了FFT算法的各种应用，并在每个实例中都附上了完整的 Matlab程序，可供读者参考。由于篇幅所 限，本文未涉及FFT变换以及其应用的数学理论背景知识。关键词：FFT算法的应用，一维基2 FFT算法，频域抽取，时域抽取，非基2 FFT算法，Cooley-Tukey算法，素因子算法，线形卷积，信号压缩和滤波，二 维FFT算法快速傅里叶变换FFT算法及其应用目 录摘要0.目录2.1 一维DFT的快速算法 一FFT31.1频域抽取的基2算法3.1.1.1正变换的计算31.1.2逆变换的计算61.2时域抽取的基2算法 7.2 一维基2 FFT算法编程8.3 一维任意非

4、基2 FFT算法123.1 COOLEY-TUKEY FFT 算法123.2素因子算法(PFA)133.3 一维任意非基2FFT算法154 一维FFT算法的应用184.1利用FFT计算连续时间信号的傅里叶变换 1 84.2利用FFT计算离散信号的线性卷积 214.3利用FFT进行离散信号压缩 234.4利用FFT对离散信号进行滤波 264.5利用FFT提取离散信号中的最强正弦分量 295二维DFT的快速变换算法及应用简介 345.1二维FFT变换及其算法介绍 345.2二维FFT变换算法的应用 35参考文献35附录361 .一维 DFT 的基 2FFT 算法 Visual C+程序36(1)

5、频域抽取的FFT和IFFT算法 36(2) 时域抽取的 FFT和IFFT算法 412. 一维任意非基2FFT算法Visual C+程序463. 二维 DFT 的基 2FFT 算法 VISUAL C+程序5 14. 二维DFT的任意非基2FFT算法VISUAL C+程序591 一维DFT的快速算法一FFT当序列f n的点数不超过N时，它的N点DFT定义为N三k nF k = £ f n e0 k N 1（ 1）n -0反变换ID FT定义为1 N丄i2兀knf nF kNeQ n N 1（2）N k -Q二者形式相似，快速算法的原理一样，这里先就其正变换进行讨论。令WN =6上卅，当k

6、依次取为0,1, 2,，N -1时，可表示为如下的方程组：”F 0 U fQ Q(WQ + f0 1W +0 2f+ f NN _Q (1N WF1 U fqT + f讪1+f 個*+ f N脚1）讣2 U fa硒° + fWh1 +f 個孙+ f N脚1）F N 1 =f oW 一1)Q+ f斷1却+ fN-NW 2）（(3)1 )由上式可见，直接按照定义计算N点序列的N点DFT时，每行含N个复乘 和N个加，从而直接按定义计算点的总计算量为 N2个复乘和N 2个加。当N较 大时，N2很大，计算量过大不仅耗时长，还会因字长有限而产生较大的误差， 甚至造成计算结果不收敛。所谓快速傅里叶

7、变换就是能大大减少计算量而完成全 部点计算的算法。下面介绍两种经典的 DFT的快速算法：频域抽取的FFT算法和时域抽取的FFT算法1.1频域抽取的基2算法1.1.1正变换的计算这里仅介绍基2算法，即是2的整次幕的情况。由定义N丄(4)F k = 7 f n WNknHk 岂 N1n -Q把f n分成两半，即f n和f n N / 2 （0空n乞N / 2 _1），代入（4）式得N / 2_ 1N L 21_kn_k(n N / 2 )F k - v f nWn v f n N / 2Wn0 乞 k 空 N _1(5)n 羽n _0由于k (n -N /2)N=Wn"wkN / 2=(

8、_1)kwknN9(5)式两项又可合并为N / 2_ 1 Fk=v fn -(k 1 )f n N 人2切< 0k N 1( 6)n £当 k =2r 为偶数时，注意到(_1)k =1，W,n =Wn2" =e"rn/N =W(n/2，( 6)式 变为N / 2 1T .rnF2 门二' (f n f nN/ 2) Wn/2nN /2 A二g(n)W/2(7)n =Q=G (r)0 乞 r 乞 N / 2 -1当k =2r 1为奇数时,（6）式变为(8)kn( 2r 1) n i2 :(2r 彳)n/ Nn rnWN =WN=e=WNWN/2，N /

9、空1F2r ( fn - f n N / 2) WnW； 2n卫N /2 -1=p(n )WNrn/2=P(r) 0_r_N/2-1n zQ这样就把一个N点序列（f n）的N点DFT （ Fk）计算化成了两个N /2 点序列（gn和pn）的N / 2点DFT （ Gr和Pr）计算。由f n划分成gn 和pn的计算量为N个加，即f n - fn N / 2和 f n - f n N / 2, 0 < n 乞 N / 2 1和N /2个乘，即（f n - f n N / 2） W,0 乞 n 乞 N /2 -1由于gn算出的N /2点D FT，是fn的N点DFT （ Fk）中k为偶数的那一半

10、，由pn算出的则是k为偶数的那一半，故需要把偶数k的Fk抽出来放 在一起作为gn的DFT （ G（r）输出，同时把奇数k的Fk抽出来放在一起作 为pn的DFT （ P（r）输出。由于k是频域变量，故这种算法称为频域抽取的 FFT算法。接着，两个N / 2点DFT仍可用上述方法各经N / 4个乘N / 2个加划分成两个 N/4点DFT (同时还要做相应的频域抽取)，从而共划分成4个N /2点D FT， 总划分计算量仍是N个加和N/2个乘。这样的划分可一步步做下去，不难看出, 每步的总划分计算量都是N个加和N /2个乘。经过M _1步的划分后就划成了 N/2个如下两点DFT的计算问题0 .00 1

11、飞A = aW2 bW2 a b22J( 9)B 二 aW2 - bW2(a -b) W2上式所需计算量是2个加和1个乘，于是完成N /2个两点DFT的总计算量仍是N个加和N / 2个乘。从而完成全部N点DFT的总计算量M N二N log 2 N 个加和M N / 2 =( N / 2) log 2 N个乘，这比直接按定义计算所需的N 2个乘和加要 少得多。例如，N = 210 =1024，M =10，用FFT算法计算所需的乘法个数为M N/ 2 =5 1024，而直接按定义计算所需的乘法个数为 N 2 =1024 1024，二 者相差1024 1 5 ：、200倍。若直接计算需半小时，而用F

12、FT计算只需9s即可完成,可见其效率之高，而且N越大，FFT的效率提高越明显。f0f1f2f3f4f5f6f7138*W8°一W4=1V-W21W81/ 乂 0V/-1八/20亠W401>PnW4W40W2°4F0000F0F0000F4100F2F1001F2010F4F2010F6 110F6F3011F1001F1F4100F5 101F3F5101F3011F5F6110F7111F7F71111W20 图1频域抽取的8点FFT计算流图一般情况下，由于做了 M -1次分奇偶的抽取，此算法最后的N/2个两点 DFT计算出的Fk不是顺序抽取的。次序的变化可用二进码

13、来说明：第一次抽 取所分的奇偶是由二进码第1位是1或0来区分的，该位为0时为偶数，该位为 1时为奇数，第二次抽奇偶是由二进码第 2位是1或0来区分的”，每次抽取都是把偶数项放在前（左）边，把奇数项放在后（右）边，从而抽取以后数的二 进码是按照二进制位从左向右依次排列的，和普通二进制数从右向左依次排列的的规律正好相反，所以称为倒位序。在计算出Fk之后要把倒位序变成顺序。1.1.2逆变换的计算所谓逆变换是指由F k求f n的计算，若直接按定义/ N 11knf n=瓦 F kWN_0 兰 n 兰 N -1N k -Q做计算，则除了求和号和正变换相同的计算量外， 每算一个fn都还需再多做一 个乘1

14、/ N的乘法运算。故按定义完成全部 N点DFT的总计算量是N _0b (A -BW2 ) 个加和N（N 1）个乘。下面从图导出它的快速算法，先讨论第3列的2点DFT的逆运算如何完成。由式（7）得，"a + b = A 0（a b） W2= B由上式不难解出1 q a (A BW2 )(10)F0000F0F1001F2F2010F4F3011F6F4100F1F5101F3F6110F5F7111F7F0000 1/8 F4100 1d-W2-0F2010 1/8 - yW4F6110 1/8 .-W2-0 -1 W4-1F1001F5101 1/8，W2F3011 1/8 <

15、W4-0F7111 1/8 农 Wa-% W4-1W80W8-2W3f0f1 f2 f3 f4 f5 f6 f72图2频域抽取的8点IFFT计算流图快速傅里叶变换FFT算法及其应用此计算过程如图2所示，可以看出：左边各列的划分计算也都是由N/2个碟形运算来完成的，只是各碟形运算所乘的相移因子W不同。把每个碟形运算都用图的办法变成对应的逆运算， 并把它们按输入在左、输出在右重新排列，就 得到了全部N点ID FT的计算流图。给出了 N =8的示例，图中先对顺序输入的 Fk做M _1次的频域抽取，并把3个乘1 /2的运算合成了一个乘1/8的运算放在 了最前边，然后就开始做求逆的碟形运算。1.2时域抽

16、取的基2算法比较正变换DFT和反变换ID FT的定义式N丄knF k - a fnWN 0 _ k _ N -1n -01 N丄fn= 瓦 FkwfnN k z0可见，去掉乘1/N的运算，把w 换成W，交换Fk、fn和k、n，反变 换定义式就变成了正变换的定义式。对图 2做这些变换，则得到图3的流程图。 对图1做这些变换，则得到图4的流程图。这就是时域抽取的算法流图。 进行碟 形运算之前，先要对顺序的时域输入序列进行 M -1次的奇偶抽取，故称为时域抽取的FFT算法f0000f0f0000f1001f2f4100f2010f4f2010f3011f6f6110f4100f1f1001f5101

17、f3f5101f6110f5f3011f7111f7f7111W20-1.V /W40yw20-1wF-1wWs°w2?亠VWs!w402Ws2-1W20-14W41Ws3-1F0F1F2F3F4F5F6F715图3时域抽取的8点FFT计算流图.F0000-1 W2-° F4100F2010二-1 亠 W2-0 F6110_F1001-1_Wt°. F5101W4-0 F3011W4-1-1 W2-0 F7111F0F0000F2F1001F4F2010F6F3011F1F4100F3F5101F5F6110F7F7111比较图2和图3不难看出，两种算法的计算量是

18、完全一样的。这里先算出 N /2个两点的DFT0 0 0 1A = f (n)W2f (n N /2)W2 1 = f (n) f (n N / 2)1 ,01 1/ A AB = f (n)W2 f (n N /2)W2 = f (n) _ f (n N / 2)( 11)f01d f11/8. f21/8 f31/8 . f41/8 * f1d f61/8 f71/8 图4时域抽取的8点IFFT计算流图然后把N /2个两点的DFT组合成N /4个4点的D FT，再把N / 4个4点的DFT组合成N /8个8点的DFT，经过M -1次的组合之后，就得到了顺序N点D FT计算结果。算法原理参考

19、文献【4】2 一维基2 FFT算法编程以频域抽取的基2 FFT正变换为例，对FFT的信号流图进行讨论，以找到 FFT算法的规律。1) 分级在进行DFT变换的过程中，从N点DFT到两点D FT共分了 M级，如图1 所示，从左到右依次为m =1级，m = 2级，m = M级。2) 倒位序在频域抽取的基2 FFT算法中，输出数据不是按照序列的先后顺序排列的，这是由于变换过程中，输出按奇、偶抽取的缘故。如果将序列xn中标号n用二进制值（n°“m訥m）2表示，那么在FFT信号流图输入端，xn位于（Rm Anm<n。"处，称为倒序。以8点FFT为例，顺序和倒序的关系如表1所示 表

20、1顺序和倒序对照表顺序倒序十进制数二进制数二进制数十进制数00 0 00 0 0010 0 11 0 0420 1 00 1 0230 1 11 1 0641 0 00 0 1151 0 11 0 1561 1 00 1 1371 1 11 1 17从表1可以看出，一个自然顺序二进制数，是在最低位加1,逢2向左移位; 而倒序数的顺序是在最高位加1,逢2向右移位。用i表示顺序数，j表示倒序 数，k表示位权重。对于一个倒序数j来说，下一个倒序数可以按下面的方法求： 先对最咼位加1，相当于十进制运算j N / 2。如果j ：: N / 2，说明二进制最咼 位为0,则直接由j N /2得到下一个倒序值

21、；如果j _ N / 2，说明二进制最高位为1，则将j的最高位变为0，通过j 二j 一 N / 2实现，同时令k 二k /2，接 着判断次高位是否为0,直到位为0时，令j二j "。归结起来算法流程图 如图5所示。j=N/2i=1 to N-2t=fi fj=fi fi= fj1k=N/2b>j=j-k k=k/2j=j+k图5倒位序程序流程图出图6频域抽取FFT程序流程图3）蝶形运算单元和同址计算频域抽取的信号流图中，基本的运算结构如图7所示，该运算结构的形状像 蝴蝶，故称为“蝶形运算单元”。Fm P二 FmP FmqrFm P =（FmP Fmq）WN在这一结构中，DFT和I

22、DFT运算关系分别为fmP"E（ 12）.fm【P =（fmP fmqWN_）/2快速傅里叶变换FFT算法及其应用(b)两点IDFT的计算图7频域抽取FFT算法流图中第 m级碟形单元而在时域抽取的信号流图中，基本的运算结构如图 8所示。在这一结构中,DFT和IDFT运算关系分别为.rFm【P二 Fm i P - Fm igW 2 - -rFm【P二 FmP- FmqWNt【P =（fm【P - fmjq） /2_rfm【P =（fmP - fmq）WN /2（13）Fm-1 prFm-1 qW N F” Fmp” F mq(b)两点IDFT的计算图8时域抽取FFT算法流图中第m级碟形

23、单元23其中，P、q分别表示该蝶形运算单元的上下节点的序号。可以看出参与运 算的输入序号P，输出序号仍为q，并且该运算不涉及到其它的点，因此我们可 以将输出的结果仍然放在数组中，称这样的操作为同址计算。也就是说，共同占有同一个存储单元。4) 寻址和相移因子WN的计算时域抽取基2 FFT信号流图中，每一级有个蝶形单元。每一级的个蝶形单元 又可以分为若干组，每一组具有相同的结构和因子的分布。如图1所示，第1级分为1组，第2级分为2组，第m级分为2m"组。在第m级中，相邻组之间的间距（也即每个分组所含节点数）为 2M，每个蝶形单元的上下节点之间的距离（也即每个分组所含碟形单元数）为 2M。

24、每组 的相移因子Wn"=cos（2 r） is in （2 r），其中 r=（l-1）2m-1 ,丨=0,1,，2“一1NN综合以上各步骤，得到频域抽取 FFT程序流程图如图6所示。采用类似的 步骤可得到频域抽取IFFT流程图、时域抽取FFT与IFFT流程图（略）。频域抽取FFT算法、时域抽取FFT算法的Visual C+源程序分别见附录1.（1），1.（2）。在Matlab中，傅里叶变换及其逆变换分别用dwt（）和idwt（）函数实现。3 一维任意非基2 FFT算法3.1 Cooley-Tukey FFT 算法FFT的核心是将一层运算映射为二层运算。 作N点FFT时，若N不是素数,

25、则N可分解为N =：心“2，那么由fn的DFTN A.(14)nkF【k =» f nWn z0通过映射:< nW NJ 0 Wn en2l兰kWMJ0 兰&兰n2匸(15)可得到nkNn 2）（峪四飞？）(N 2 nk1 -'N! N 2门1 k2 n 2$ M n 2X2)NWnN1 =Wn2，WnN2 -Wn1，可化简为nkn2k1Nn2k2N2(16)从而式（14）转化为N 2 -1(W,2k1 瓦 fgn 2WN：k1)(17)n?k2n?出其中 0乞匕乞2 1,0乞k2乞N2 T。以20点FFT为例，N =20,心=5, N2 =4，映射方式为：n=

26、4nJ n?，k = & 5k2，则计算流图如图9所示3.2 素因子算法(Prime Factor Algorithm, PFA)当Cooley-Tukey FFT算法中的、n?互素的话，相移因子WNn2k1可通过选择 m,n2, ki,k2前的特殊系数而消去，FFT的算法进一步的简化。n =A nB2n<1 n 空1N1,0乞 2n空N1-12_1（18）k =C kD2k 0乞飞乞八1一，0乞2k乞N1当A、B、C、D满足以下条件AD 三 0 m od N（19）BC 三 0 m od NF0F5F10F15F1F6F11F16F2F7F12F17F3F8F13F18F4图

27、9 Cooley-Tukey 20 点 FFT 算法F9F14F19时，有AC 三 NmodNBD三 Nmod N(19)nkN(An i -|Bn 2(Ck -p k ) N(AC nik lAD n ki “BC n k -2BDin k ) 22(20)(21)二 WN-WN：*22于是式(14)化简为“2 丄Ni J卩凶也=£ Wn?2(送 fn 1, n2】WN：k1)n2 =0n1 =0从而消去了相移因子WNn2k1。同样以20点FFT为例，修改映射方式为：N = 20, N <, = 5, N 2 = 4n =4n5 n2, 0 _n _ 4, 0 _n 2 _3

28、（ 22）k =16匕 25k2,0 乞匕乞4, 0< k2 乞 3（23）则由式（22）得到的映射关系如表2，由式（23）得到的映射关系如表3, 计算流图如图10所示。表2由式（22）得到的映射关系012300510151491419281318331217274161611表3由式（23）得到的映射关系0123005101511616112121727381318344914193.3 维任意非基2 FFT算法对于非素数N点DFT，对N做因式分解N = NN Nl（ 24）令N21二N2Nl,则N = N 1 N 21。于是式（24 ）中多层FFT转化为二层FFT。如 果Ni与N21

29、互素，那么采用PFA算法进行映射，否则采用 Cooley-Tukey FFT算 法（此时需乘以相移因子）。N2l二N22可采用同样的方法分解成N2个N3点 DFT，其中N31 = N3N|，依此类推，直到DFT长度为N |。可以证明，复乘的总次数不大于n（nn2 亠 亠-丨）（25）例如，若N =63 =3 3 7，则复乘总次数至多为63 （33 7 _ 3） = 630。而用基2的FFT算法计算64点DFT，需要的复乘总次数为192。这说明，将N分 解得越细，运算量越少。实际中，常常将输入序列补零，使得N成为2的幕次数，这样能够减少复乘运算的次数。再以20点FFT为例，进行如下三层因式分解即

30、Ni =5，N23=2 2=4，由于5与4互素，所以第一层采用PFA算法，相应的 二层映射为n = 4nt 5n 23, 0 _n 勺 _4,0 _n 23 _3k =16 k125 k23,0 乞 k 4, 0 < k23 乞 3由于2与2不互素，所以第二层采用 Cooley-Tukey FFT算法，相应的二层映射为n23 =2n2 n3, 0 _n2 -1,0 _n3_1快速傅里叶变换FFT算法及其应用45k23 = k2 - 2k3, 0 乞 k2 乞 1,0 乞 k3 乞 1总的FFT变换公式如式（26）,计算流图如图11所示f0 f4 f8 f12 f16n1-2_J亠4 _k

31、1 0 1234f5 f9 f13 f17 f1f10 f14 f18 f2 f60*1*234f15 f19 f3 f7 f11n23=0n3=10Jn3=01n3=1n3=1n3=03n3=0n3=11*3=4 ”3401On3=01n2-0*0： 1一23 +0*40120*0*1*n3=01k2n30 W40 0 -W40W/411*k301F0F10W401F5F150w4°0*W401 wJ 1*010101W4°0*W40W4W4J 1J010 一.F171 F7F16F6F1F11F12F20W4°0*W400 F81-, F18n 3=1-0*

32、W40*0*1W411*0r1 rF13F301W41W41UW40 0*W400 F41，F140 = F91 F19图11多层分解时20点FFT算法N3 -1I=WnN 2 -1k 3n 3k 2 n3Wn：2、WnN1 二"(瓦 fn!,n2,n 3叭丁)114八 Wk3W4n3k2 7 Wk2(y fni，n2，n3Wk1)(26)门3 -0门2 £ni _0比较正变换DFT和反变换ID FT的定义式可知，在正变换前加上乘1/ N的运算，把W换成w -，并交换fn、Fk和n、k，就得到反变换的算法。一维任意非基2 FFT算法的Visual C+源程序见附录2。在Ma

33、tlab中，傅里叶变换及其逆变换也分别用dwt()和idwt()函数实现。4 一维FFT算法的应用4.1利用FFT计算连续时间信号的傅里叶变换设x(t)是连续时间信号，并假设t : 0时x(t) =0，则其傅里叶变换由下式给 出X( J 二 0 x(t)e，dt令】是一个固定的正实数，N是一个固定的正整数。当：，k =0,1, 2,N -1时，利用FFT算法可计算X C )。已知一个固定的时间间隔T，选择T足够小，使得每一个T秒的间隔 nT <t <(n 1)T内，x(t)的变化很小，则式中积分可近似为近(n+)TX ()- 7 (e dt)x(nT ) nTn z0n i'

34、; 上 nr ) n TTX (n T)_eT 二=；' e" nTx(nT )(27)i 1n假设N足够大，对于所有n AN的整数，幅值x(nT )很小，则式(27)变为工T N 4X() =ZnT )( 28)in T当川=2二k/ NT时，式(28)两边的值为i_2 飞/N N-12 ；k /N2兀 k1 e _12刑小1 eX( )e - x(nT)Xk( 29)NT i2二 k / NT n 占i2二 k/ NT其中Xk代表抽样信号xn =x( nT)的N点DFT。最后令】=2二/NT，则上式变为i2 卄 /N(30)1 _e 71X (kl)Xk k =0,1,

35、2, N -1i2jrk/ NT首先用FFT算法求出Xk,然后可用上式求出k =0,1, 2,N _1时的X （k ：）。应该强调的是，式（28）只是一个近似表示，计算得到的 x（j只是一个近 似值。通过取更小的抽样间隔T，或者增加点数N ,可以得到更精确的值。如果 B时，幅度谱X C .）很小，对应于奈奎斯特抽样频率.2B，抽样间隔T选 择二/ B比较合适。如果已知信号只在时间区间0学 汀 内存在，可以通过对nT . t时的抽样信号xn = x（nT ）补零，使N足够大N=i nput('l nput N:');T=in put('I nput T:');%计

36、算X(w)近似值t=0:T:2;x=t-1 zeros(1,N-le ngth(t);X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xap p=(1-exp(-i*k*gamma*T)/(i*k*gamma)*X;%计算真实值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-si n(w)./(w.*w); plot(k*gamma,abs(Xapp(1:le ngth(k),'o',w,abs(Xact); lege nd('近似值','真实值');xlabel(

37、'频率(rad/s)');ylabel('|X|')运行程序后输入N=128，T=0.1,此时丨= 0.4909，得到实际的和近似的傅里叶变换的幅度谱如图13所示，此时近似值已经相当准确。通过增加 NT可以增 加更多的细节，减少T使得到的值更精确。再次运行程序后输入N=512, T=0.05, 此时卜=0.2454，得到实际的和近似的傅里叶变换的幅度谱如图 14所示。4.2利用FFT计算离散信号的线性卷积已知两个离散时间信号 xn (n =0,1, 2，M -1)与y n (n =0,1, 2，N _1)，取 L = M N -1，对xn和yn右端补零，使得x

38、n =0, n =M 1, M 2,L -1y n二 0 ,n 二 N 1 ,N 2 , U-( 31)利用FFT算法可以求得xn和yn的L点DFT，分别是X k和Yk，利用 DTFT卷积性质，卷积xn* yn等于乘积XkYk的L点DFT反变换，这也可 以通过FFT算法得到。例2利用FFT计算线性卷积已知xn =0.8n，其中u n为单位阶跃序列，信号y n如图15所示。由 于当n 16时，xn很小，故M可以取为17； N取10， L =M N -26。利用下面的Matlab命令，可得到xn、yn的卷积图形如图15所示。subplot(3,1,1);n=0:16;x=0.8.A n;stem(

39、 n,x); xlabel(' n');ylabel('x n'); subplot(3,1,2);n=0:15;y=ones(1,10) zeros(1,6); stem( n, y);xlabel(' n');ylabel('y n') subplot(3,1,3);L=26 ;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y ;z=ifft(Z,L); stem( n,z);xlabel(' n');ylabel('z n')图 15 信号 xn、yn及其卷积 zn=x

40、n*yn利用下面的Matlab命令，可得到信号xn、yn的幅度谱与相位谱如图16 所示。subplot(2,2,1);L=26;k=0:L-1; n=0:16;x=0.8。n;X=fft(x,L); stem(k,abs(X);axis(0 25 0 5); xlabel('k');ylabel('|Xk|')subplot(2,2,2); stem(k,a ngle(X);axis(0 25 -1 1);xlabel('k');ylabel('A ngle(Xk)(弧度)') subplot(2,2,3);y=o nes(1,1

41、0);Y=fft(y,L); stem(k,abs(Y);axis(0 25 0 10);xlabel('k');ylabel('|Yk|') subplot(2,2,4);stem(k,a ngle(Y); axis(0 25 -3 3);xlabel('k');ylabel('A ngle(Yk)(弧度)')图16信号xn、yn的幅度谱与相位谱4.3利用FFT进行离散信号压缩利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化;2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于 某一阈值的

42、系数置为0,保留剩余的系数；4）利用IFFT算法对处理后的信号进 行逆傅里叶变换。例3 对单位区间上的下列连续信号1f (t) =t ；：；'cos(4 二t) sin(8 二t)2以fs = 256 Hz采样频率进行采样，将其离散化为28个采样值fn二 f(t)|twT 二 f(nT )二 f ( n / 256), n 二 0,1, 2,，255用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0,得到重构信号图形如图17 a)所示，均方差为0.0429,相对误差为0.0449;令绝对值最小的90%系数为0,得到重构信号图形如图17 b)所示，均方差为0

43、.0610,相对误差为0.0638。a)绝对值最小的80%系数为0的重构信号(FFT) b)绝对值最小的90%系数为0的重构信号(FFT)图17用FFT压缩后的重构信号相关Matlab程序如下function wc=compress(w,r)%压缩函数compress.m%输入信号数据w,压缩率r%输出压缩后的信号数据if(r<0)|(r>1)error('r应该介于0和1之间!');en d;N=le ngth(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)>=tol).*w;funct

44、ion un biased_varia nce,error=fftcomp(t,y,r)%利用FFT做离散信号压缩%输入时间t,原信号y,以及压缩率r%输出原信号和压缩后重构信号的图像，以及重构均方差和相对lA2误差 if(r<0)|(r>1)error('r应该介于0和1之间!')；en d; fy=fft(y);fyc=compress(fy,r); %调用压缩函数 compress.m yc=ifft(fyc);plot(t,y,'r',t,yc,'b');lege nd('原信号','重构信号'

45、);un biased_varia nce=norm(y-yc)/sqrt(le ngth(t);error= norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*si n(8*pi*t);un biased_varia nce,error=fftcomp(t,f,0.8)un biased_varia nee =0.0429error =0.0449如果用Harr尺度函数和Harr小波分解信号，对信号进行小波压缩，然后重 构信号。令绝对值最小的80%系数为0,得到重构信号图形如图18 a)所示，均方 差为0.05

46、84，相对误差为0.0611 ;令绝对值最小的90%系数为0,得到重构信号图形如图18 b)所示，均方差为0.1136，相对误差为0.1190a)绝对值最小的80%系数为0的重构信号(Harr)原信号 畫构曲号b)绝对值最小的90%系数为0的重构信号(Harr)图18用Harr小波压缩后的重构信号相关Matlab程序如下function unbiased variance,error=daubcomp(t,y ,n,r)%利用Daubechies系列小波做离散信号压缩%输入时间t,原信号y,分解层数n,以及压缩率r%输出原信号和压缩后重构信号的图像，以及重构均方差和相对lA2误差if(r<

47、;0)|(r>1)error('r应该介于0和1之间!');en d;c,l=wavedec( y,n ,'db1');cc=compress(c,r);%调用压缩函数 compress.myc=waverec(cc,l,'db1');plot(t,y,'r',t,yc,'b');lege nd('原信号','重构信号');un biased_varia nce=norm(y-yc)/sqrt(le ngth(t);error =norm(y yc)/no rm(y);输入以下

48、Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*si n(8*pi*t);un biased_varia nce,error=daubcomp(t,f,8,0.8)un biased_varia nee =0.0584error =0.0611结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好(甚至比小波还要好)。4.4利用FFT对离散信号进行滤波利用FFT算法对信号进行滤波的步骤如下：1)通过采样将信号离散化；2) 对离散化信号进行傅里叶变换；3)对变换后的系数进行处理，将绝对值大于某 一阈值的系数置为0,保留剩余的系数；4)

49、利用IFFT算法对处理后的信号进行 逆傅里叶变换。例4股票价格分析首先进入网址 Download To Spreadsheet按钮，即可把以Excel表格格式存储的价格数据下载到 本地计算机。表格从1列至第6列分别给出了从1996年4月12日至2007年5月30日的交易期里每天的开盘价、最高价、最低价、收盘价、成交量以及趋势。 数据下载完成后，需要颠倒顺序，使得最早时间的数据首先显示。然后另存到 Matlab所在的目录中，并重新命名为“yhoodata.csV'。本次分析选择开盘价，时间是从2007年1月1日至2007年5月30日的N =102个交易日期。令on, n =0,1,，N

50、_1代表一支股票的开盘价。为了便于分析，可以先从on中减去跃变，得到xn, n = 0,1，,N -1，即o N 1 匸 o 0 /“、xn= onj on ,= n0,1,_N, 1( 32)N -1输入以下命令，得到xn的频谱如图19所示。o=csvread('yhoodata.csv',2700,1,2700 1 2801 1)'N=102;for n=1:Nx(n )=o( n)-o(1)-(o(N)-o(1)/(N-1)*( n-1);endX=fft(x);k=0:N-1;stem(k,abs(X);axis(0 101 0 300);xlabel('k');ylabel('|Xk|')图佃 xn的幅度谱可以看出上图中有5个较强谱分量k =0,1, 4, 98,100，频率(2二k/N )分别 对应0, 2二/102和8二/102。保留这5个频率分量的系数，将其他频率分量的系数 置为0,然后再进行逆傅里叶变换，得到滤波后的近似值Xn。输入如下Matlab程序，得到真实值xn与滤波后的近似值 刘n，如图20所示。plot(x);hold on;fliterX=X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0 X(102);fliterx