快速傅里叶变换FFT算法及其应用_————.doc

快速傅里叶变换 FFT 算法及其应用 摘 要 本文较为系统地阐述了快速傅里叶变换的算法原理及其在数字信号处理等 工程技术中的应用 根据抽取方法的不同 一维基 2 FFT 算法分为两种 频域 抽取的 FFT 算法和时频域抽取的 FFT 算法 第 1 节阐述了这两种 FFT 算法的 原理 第 2 节给出了两种算法的编程思想和步骤 第 3 节阐述了一维非基 2 FFT 的两种算法 Cooley tukey FFT 算法和素因子算法 Prime Factor Algorithm 的思想原理 给出了在把一维非基 2 DFT 的多层分解式转化为二层分解的过程 中 如何综合运用这两种算法以达到总运算次数最少的方案 并以 20 点 DFT 为例描述了非基 2 FFT 算法实现的一般步骤 第 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 编程 通过大量生动的实例 图文并茂地列举出了 FFT 算法的各种应 用 并在每个实例中都附上了完整的 Matlab 程序 可供读者参考 由于篇幅所 限 本文未涉及 FFT 变换以及其应用的数学理论背景知识 关键词 FFT 算法的应用 一维基 2 FFT 算法 频域抽取 时域抽取 非 基 2 FFT 算法 Cooley Tukey 算法 素因子算法 线形卷积 信号压缩和滤波 二维 FFT 算法 快速傅里叶变换 FFT 算法及其应用 1 摘 要 2 目 录 摘 要 1 目 录 2 1 一维 DFT 的快速算法 FFT 4 1 1 频域抽取的基 2 算法 4 1 1 1 正变换的计算 4 1 1 2 逆变换的计算 7 1 2 时域抽取的基 2 算法 8 2 一维基 2 FFT 算法编程 10 3 一维任意非基 2 FFT 算法 14 3 1 COOLEY TUKEY FFT 算法 14 3 2 素因子算法 PFA 14 3 3 一维任意非基 2 FFT 算法 16 4 一维 FFT 算法的应用 20 4 1 利用 FFT 计算连续时间信号的傅里叶变换 20 4 2 利用 FFT 计算离散信号的线性卷积 23 4 3 利用 FFT 进行离散信号压缩 25 4 4 利用 FFT 对离散信号进行滤波 28 4 5 利用 FFT 提取离散信号中的最强正弦分量 31 5 二维 DFT 的快速变换算法及应用简介 36 5 1 二维 FFT 变换及其算法介绍 36 5 2 二维 FFT 变换算法的应用 37 参考文献 38 附 录 39 1 一维 DFT 的基 2 FFT 算法 VISUAL C 程序 39 1 频域抽取的 FFT 和 IFFT 算法 39 2 时域抽取的 FFT 和 IFFT 算法 44 2 一维任意非基 2 FFT 算法 VISUAL C 程序 49 快速傅里叶变换 FFT 算法及其应用 3 3 二维 DFT 的基 2 FFT 算法 VISUAL C 程序 54 4 二维 DFT 的任意非基 2 FFT 算法 VISUAL C 程序 62 1 一维 DFT 的快速算法 FFT 当序列的点数不超过时 它的点定义为 f nNNDFT 2 1 0 01 N ikn N n F kf n ekN 1 反变换定义为IDFT 2 1 0 1 01 N ikn N k f nF k enN N 2 二者形式相似 快速算法的原理一样 这里先就其正变换进行讨论 令 当依次取为时 可表示为如下的方程组 2 iN N We k0 1 2 1N 0001020 1 1011121 1 2021222 1 1 0 1 1 1 1 0 0 1 2 1 1 0 1 2 1 2 0 1 2 1 1 0 1 1 N NNNN N NNNN N NNNN NNNN NNN FfWfWfWf NW FfWfWfWf NW FfWfWfWf NW F NfWfWf NW 3 由上式可见 直接按照定义计算点序列的点时 每行含个复NNDFTN 乘和个加 从而直接按定义计算点的总计算量为个复乘和个加 当N 2 N 2 N 较大时 很大 计算量过大不仅耗时长 还会因字长有限而产生较大的N 2 N 误差 甚至造成计算结果不收敛 所谓快速傅里叶变换就是能大大减少计算量 而完成全部点计算的算法 下面介绍两种经典的的快速算法 频域抽取的DFT 算法和时域抽取的算法 FFTFFT 1 一维 DFT 的快速算法 FFT 4 1 1 频域抽取的基 2 算法 1 1 1 正变换的计算 这里仅介绍基 2 算法 即是 2 的整次幂的情况 由定义 1 0 01 N kn N n F kf n WkN 4 把分成两半 即和 代入 4 式得 f n f n 2 f nN 0 2 1 nN 2 1 2 1 2 00 2 01 NN knk n N NN nn F kf n Wf nNWkN 5 由于 2 2 1 k n NknkNkkn NNNN WW WW 5 式两项又可合并为 2 1 0 1 2 01 N kkn N n F kf nf nNWkN 6 当为偶数时 注意到 6 2kr 1 1 k 22 2 knrnirn N NN WWe 2 rn N W 式变为 7 2 1 2 0 2 1 2 0 2 2 0 2 1 N rn N n N rn N n Frf nf nNW g n W G rrN 当为奇数时 6 式变为21kr 21 2 21 2 knrnirn Nnrn NNNN WWeW W 8 2 1 2 0 2 1 2 0 21 2 0 2 1 N nrn NN n N rn N n Frf nf nNWW p n WP rrN 这样就把一个点序列 的点 计算化成了两个N f nNDFT F k 点序列 和 的点 和 计算 由划分 2N g n p n 2NDFT G r P r f n 成 g n 和的计算量为个加 即 p nN 快速傅里叶变换 FFT 算法及其应用 5 和 2 f nf nN 2 0 2 1f nf nNnN 和个乘 即 2N 2 0 2 1 n N f nf nNWnN 由于算出的点 是的点 中为偶数的 g n 2NDFT f nNDFT F kk 那一半 由算出的则是为偶数的那一半 故需要把偶数的抽出来 p nkk F k 放在一起作为的 输出 同时把奇数的抽出来放在一起 g nDFT G rk F k 作为的 输出 由于是频域变量 故这种算法称为频域抽取 p nDFT P rk 的算法 FFT 接着 两个点仍可用上述方法各经个乘个加划分成两 2NDFT 4N 2N 个点 同时还要做相应的频域抽取 从而共划分成 4 个点 4NDFT 2N 总划分计算量仍是个加和个乘 这样的划分可一步步做下去 DFTN 2N 不难看出 每步的总划分计算量都是个加和个乘 N 2N 经过步的划分后就划成了个如下两点的计算问题1M 2NDFT 9 0 001 22 10110 222 AaWbWab BaWbWab W 上式所需计算量是 2 个加和 1 个乘 于是完成个两点的总计算量 2NDFT 仍是个加和个乘 从而完成全部点的总计算量N 2NNDFT 个加和个乘 这比直接按定义计算所 2 logMNNN 2 2 2 logMNNN 需的个乘和加要少得多 例如 用算法计算 2 N 10 21024N 10M FFT 所需的乘法个数为 而直接按定义计算所需的乘法个数为 2MN 5 1024 二者相差倍 若直接计算需半小时 而用 2 1024 1024N 1024 5200 计算只需 9s 即可完成 可见其效率之高 而且越大 的效率提高FFTNFFT 越明显 1 一维 DFT 的快速算法 FFT 6 f 0 F 0 000 F 0 F 0 000 f 1 1 W20 F 4 100 F 2 F 1 001 g n f 2 1 W40 F 2 010 F 4 F 2 010 f 3 1 W41 1 W20 F 6 110 F 6 F 3 011 f 4 1 W80 F 1 001 F 1 F 4 100 f 5 1 W81 1 W20 F 5 101 F 3 F 5 101 p n f 6 1 W82 1 W40 F 3 011 F 5 F 6 110 f 7 1 W83 1 W41 1 W20 F 7 111 F 7 F 7 111 图图 1 频域抽取的频域抽取的 8 点点 FFT 计算流图计算流图 一般情况下 由于做了次分奇偶的抽取 此算法最后的个两点1M 2N 计算出的不是顺序抽取的 次序的变化可用二进码来说明 第一次抽DFT F k 取所分的奇偶是由二进码第 1 位是 1 或 0 来区分的 该位为 0 时为偶数 该位 为 1 时为奇数 第二次抽奇偶是由二进码第 2 位是 1 或 0 来区分的 每次 抽取都是把偶数项放在前 左 边 把奇数项放在后 右 边 从而抽取以后 数的二进码是按照二进制位从左向右依次排列的 和普通二进制数从右向左依 次排列的的规律正好相反 所以称为倒位序 在计算出之后要把倒位序变 F k 成顺序 1 1 2 逆变换的计算 所谓逆变换是指由求的计算 若直接按定义 F k f n 1 0 1 01 N kn N k f nF k WnN N 做计算 则除了求和号和正变换相同的计算量外 每算一个都还需再多做 f n 一个乘的乘法运算 故按定义完成全部点的总计算量是个加和1 NNDFT 2 N 个乘 下面从图导出它的快速算法 先讨论第 3 列的 2 点的逆运 1 N N DFT 算如何完成 由式 7 得 0 2 abA ab WB 快速傅里叶变换 FFT 算法及其应用 7 由上式不难解出 0 2 0 2 1 2 1 2 aABW bABW 10 F 0 000 F 0 F 0 000 1 8 f 0 F 1 001 F 2 F 4 100 1 8 W2 0 1 f 1 F 2 010 F 4 F 2 010 1 8 W4 0 1 f 2 F 3 011 F 6 F 6 110 1 8 W2 0 1 W4 1 1 f 3 F 4 100 F 1 F 1 001 1 8 W8 0 1 f 4 F 5 101 F 3 F 5 101 1 8 W2 0 1 W8 1 1 f 5 F 6 110 F 5 F 3 011 1 8 W4 0 1 W8 2 1 f 6 F 7 111 F 7 F 7 111 1 8 W2 0 1 W4 1 1 W8 3 1 f 7 图图 2 频域抽取的频域抽取的 8 点点 IFFT 计算流图计算流图 此计算过程如图 2 所示 可以看出 左边各列的划分计算也都是由个 2N 碟形运算来完成的 只是各碟形运算所乘的相移因子不同 把每个碟形运算W 都用图的办法变成对应的逆运算 并把它们按输入在左 输出在右重新排列 就得到了全部点的计算流图 给出了的示例 图中先对顺序输入NIDFT8N 的做次的频域抽取 并把 个乘的运算合成了一个乘的运算 F k1M 31 21 8 放在了最前边 然后就开始做求逆的碟形运算 1 2 时域抽取的基 2 算法 比较正变换和反变换的定义式DFTIDFT 1 0 01 N kn N n F kf n WkN 1 0 1 01 N kn N k f nF k WnN N 1 一维 DFT 的快速算法 FFT 8 可见 去掉乘的运算 把换成 交换 和 反1 N 1 W W F k f nkn 变换定义式就变成了正变换的定义式 对图 2 做这些变换 则得到图 3 的流程 图 对图 1 做这些变换 则得到图 4 的流程图 这就是时域抽取的算法流图 进行碟形运算之前 先要对顺序的时域输入序列进行次的奇偶抽取 故称1M 为时域抽取的算法 FFT f 0 000 f 0 f 0 000 F 0 f 1 001 f 2 f 4 100 W20 1 F 1 f 2 010 f 4 f 2 010 W40 1 F 2 f 3 011 f 6 f 6 110 W20 1 W41 1 F 3 f 4 100 f 1 f 1 001 W80 1 F 4 f 5 101 f 3 f 5 101 W20 1 W81 1 F 5 f 6 110 f 5 f 3 011 W40 1 W82 1 F 6 f 7 111 f 7 f 7 111 W20 1 W41 1 W83 1 F 7 图图 3 时域抽取的时域抽取的 8 点点 FFT 计算流图计算流图 比较图 2 和图 3 不难看出 两种算法的计算量是完全一样的 这里先算出 个两点的 2NDFT 0 001 22 2 2 Af n Wf nNWf nf nN 11 1011 22 2 2 Bf n Wf nNWf nf nN 快速傅里叶变换 FFT 算法及其应用 9 f 0 1 8 F 0 000 F 0 F 0 000 f 1 1 8 1 W2 0 F 4 100 F 2 F 1 001 f 2 1 8 1 W4 0 F 2 010 F 4 F 2 010 f 3 1 8 1 W4 1 1 W2 0 F 6 110 F 6 F 3 011 f 4 1 8 1 W8 0 F 1 001 F 1 F 4 100 f 5 1 8 1 W8 1 1 W2 0 F 5 101 F 3 F 5 101 f 6 1 8 1 W8 2 1 W4 0 F 3 011 F 5 F 6 110 f 7 1 8 1 W8 3 1 W4 1 1 W2 0 F 7 111 F 7 F 7 111 图图 4 时域抽取的时域抽取的 8 点点 IFFT 计算流图计算流图 然后把个两点的组合成个4点的 再把个4点的 2NDFT 4NDFT 4N 组合成个8点的 经过次的组合之后 就得到了顺序点DFT 8NDFT1M N 计算结果 算法原理参考文献 4 DFT 2 一维基 2 FFT 算法编程 以频域抽取的基 2 FFT 正变换为例 对 FFT 的信号流图进行讨论 以找到 FFT 算法的规律 1 分级 在进行变换的过程中 从点到两点共分了 M 级 如图 1DFTNDFTDFT 所示 从左到右依次为级 级 级 1m 2m mM 2 倒位序 在频域抽取的基 2 FFT 算法中 输出数据不是按照序列的先后顺序排列的 这是由于变换过程中 输出按奇 偶抽取的缘故 如果将序列中标号用 x nn 二进制值表示 那么在 FFT 信号流图输入端 位于 021 2 MM nnn x n 处 称为倒序 以 8 点 FFT 为例 顺序和倒序的关系如表 1 所 1202 MM nnn 示 2 一维基 2 FFT 算法编程 10 表表 1 顺序和倒序对照表顺序和倒序对照表 顺序倒序 十进制数二进制数二进制数十进制数 0 1 2 3 4 5 6 7 0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 0 4 2 6 1 5 3 7 从表 1 可以看出 一个自然顺序二进制数 是在最低位加 1 逢 2 向左移 位 而倒序数的顺序是在最高位加 1 逢 2 向右移位 用 表示顺序数 表示ij 倒序数 表示位权重 对于一个倒序数来说 下一个倒序数可以按下面的kj 方法求 先对最高位加 1 相当于十进制运算 如果 说明二 2jN 2jN 进制最高位为 0 则直接由得到下一个倒序值 如果 说明二 2jN 2jN 进制最高 位为 1 则将的最高位变为 0 通过实现 同时令 接j 2jjN 2kk 着判断次高位是否为 0 直到位为 0 时 令 归结起来算法流程图jjk 如图 5 所示 快速傅里叶变换 FFT 算法及其应用 11 j N 2 i 1 to N 2 t f i f j f i f i f j k N 2 j j k k k 2 kj j j k i j 是 是 否 否 输入 N 信号 f N M log2N l 1 to lb r l 1 2m 1 n l 1 la N 2 m 1 to M la 2M 1 m lb la 2 lc n lb t f n f lc f lc f n f lc WNr f n t 倒位序重排信号 图图 5 倒位序程序流程图倒位序程序流程图 图图 6 频域抽取频域抽取 FFT 程序流程图程序流程图 3 蝶形运算单元和同址计算 频域抽取的信号流图中 基本的运算结构如图 7 所示 该运算结构的形状 像蝴蝶 故称为 蝶形运算单元 在这一结构中 DFT 和 IDFT 运算关系分别为 12 11 11 mmm r mmmN FpFpFq FpFpFq W 11 11 2 2 r mmmN r mmmN fpfpfq W fpfpfq W 输 出 2 一维基 2 FFT 算法编程 12 Fm 1 p Fm p Fm 1 q 1 WNr Fm q a 两点 DFT 的计算 fm 1 p 1 2 fm p fm 1 q WN r 1 1 2 fm q b 两点 IDFT 的计算 图图 7 频域抽取频域抽取 FFT 算法流图中第算法流图中第 m 级碟形单元级碟形单元 而在时域抽取的信号流图中 基本的运算结构如图 8 所示 在这一结构中 DFT 和 IDFT 运算关系分别为 11 11 r mmmN r mmmN FpFpFq W FpFpFq W 11 11 2 2 mmm r mmmN fpfpfq fpfpfq W 13 Fm 1 p Fm p Fm 1 q WNr 1 Fm q a 两点 DFT 的计算 fm 1 p 1 2 fm p fm 1 q 1 WN r 2 fm q b 两点 IDFT 的计算 图图 8 时域抽取时域抽取 FFT 算法流图中第算法流图中第 m 级碟形单元级碟形单元 其中 分别表示该蝶形运算单元的上下节点的序号 可以看出参与pq 运算的输入序号 输出序号仍为 并且该运算不涉及到其它的点 因此我pq 们可以将输出的结果仍然放在数组中 称这样的操作为同址计算 也就是说 共同占有同一个存储单元 4 寻址和相移因子的计算 r N W 时域抽取基 2 FFT 信号流图中 每一级有个蝶形单元 每一级的个蝶形单 元又可以分为若干组 每一组具有相同的结构和因子的分布 如图 1 所示 第 1 级分为 1 组 第 2 级分为 2 组 第级分为组 m 1 2m 在第级中 相邻组之间的间距 也即每个分组所含节点数 为 每个m 1 2M m 蝶形单元的上下节点之间的距离 也即每个分组所含碟形单元数 为 每2M m 快速傅里叶变换 FFT 算法及其应用 13 组的相移因子 其中 22 cos sin r N Wrir NN m 1 r l 1 2 0 1 l 21 Mm 综合以上各步骤 得到频域抽取 FFT 程序流程图如图 6 所示 采用类似的 步骤可得到频域抽取 IFFT 流程图 时域抽取 FFT 与 IFFT 流程图 略 频域抽取 FFT 算法 时域抽取 FFT 算法的 Visual C 源程序分别见附录 1 1 1 2 在 Matlab 中 傅里叶变换及其逆变换分别用 dwt 和 idwt 函 数实现 3 一维任意非基 2 FFT 算法 3 1 Cooley Tukey FFT 算法 的核心是将一层运算映射为二层运算 作点时 若不是素FFTNFFTN 数 则可分解为 那么由的N 12 NN N f nDFT 1 0 01 N nk N n F kf n WkN 14 通过映射 2121122 1121122 01 01 01 01 nN nnnNnN kkN kkNkN 15 可得到 2 1211 22 1 112 1 22 11 2 2 N nnkN kN n kN N n kn kN n knk NNN WWW 而 可化简为 12 NN N 1 2 N NN WW 2 1 N NN WW 1 12 12 2 12 n kn kn knk NNNN WWWW 16 从而式 14 转化为 17 21 2 22 11 1 21 21 11 1212 00 NN n kn kn k NNN nn F k kWWf n n W 其中 1122 01 01kNkN 以点为例 映射方式为 20FFT 12 20 5 4NNN 12 4nnn 3 一维任意非基 2 FFT 算法 14 则计算流图如图 9 所示 12 5kkk 3 2 素因子算法 Prime Factor Algorithm PFA 当 Cooley Tukey FFT 算法中的 互素的话 相移因子可通过选 1 N 2 N 2 1 n k N W 择前的特殊系数而消去 FFT 的算法进一步的简化 12 n n 12 k k 121122 121122 01 01 01 01 nAnBnnNnN kCkDkkNkN 18 当满足以下条件ABCD 0 mod 0 mod ADN BCN 19 n1 k2 f 0 0 W0 0 F 0 f 4 1 W0 1 F 5 f 8 2 W0 2 F 10 f 12 3 W0 3 F 15 f 16 4 W0 0 F 1 f 1 0 W1 1 F 6 f 5 1 W2 2 F 11 f 9 2 W3 3 F 16 f 13 3 f 17 4 W0 0 F 2 W2 1 F 7 f 2 0 W4 2 F 12 f 6 1 W6 3 F 17 f 10 2 f 14 3 W0 0 F 3 f 18 4 W3 1 F 8 W6 2 F 13 f 3 0 W9 3 F 18 f 7 1 f 11 2 W0 0 F 4 f 15 3 W4 1 F 9 f 19 4 W8 2 F 14 W12 3 F 19 k1 0 n2 0 n2 1 n2 2 n2 3 k1 1 k1 2 k1 3 k1 4 图图 9 Cooley Tukey 20 点点 FFT 算法算法 快速傅里叶变换 FFT 算法及其应用 15 2 1 mod mod ACNN BDNN 19 时 有 1212 1 11 22 12 2 2 1 11 2 2 1 12 2 12 AnBnCkDknk NN ACn kADn kBCn kBDn k N N n kN n k N n kn k NN WW W W WW 20 于是式 14 化简为 21 2 21 1 21 21 11 1212 00 NN n kn k NN nn F k kWf n n W 21 从而消去了相移因子 同样以点为例 修改映射方式为 2 1 n k N W20FFT 20 N 12 5 4NN 1212 45 04 03nnnnn 22 1212 1625 04 03kkkkk 23 则由式 22 得到的映射关系如表 2 由式 23 得到的映射关系如表 3 计算流图如图 10 所示 表表 2 由式由式 22 得到的映射关系得到的映射关系 n1 n20123 0051015 1491419 2813183 3121727 4161611 表表 3 由式由式 23 得到的映射关系得到的映射关系 k1 k20123 0051015 1161611 2121727 n k 3 一维任意非基 2 FFT 算法 16 3813183 4491419 3 3 一维任意非基 2 FFT 算法 对于非素数点 对做因式分解NDFTN 12l NN NN 24 令 则 于是式 24 中多层转化为二层 22ll NNN 12l NN N FFTFFT 如果与互素 那么采用 PFA 算法进行映射 否则采用 Cooley Tukey FFT 1 N 2l N 算法 此时需乘以相移因子 可采用同样的方法分解成个 22ll NNN 2 N 点 其中 依此类推 直到长度为 3l NDFT 33ll NNN DFT l N 可以证明 复乘的总次数不大于 12 l N NNNl 25 例如 若 则复乘总次数至多为 而633 3 7N 63 3 3 73 630 用基 2 的算法计算 64 点 需要的复乘总次数为 192 这说明 将FFTDFT 分解得越细 运算量越少 实际中 常常将输入序列补零 使得成为 2 的NN 幂次数 这样能够减少复乘运算的次数 快速傅里叶变换 FFT 算法及其应用 17 n1 k2 f 0 0 0 F 0 f 4 1 1 F 5 f 8 2 2 F 10 f 12 3 3 F 15 f 16 4 0 F 16 f 5 0 1 F 1 f 9 1 2 F 6 f 13 2 3 F 11 f 17 3 f 1 4 0 F 12 1 F 17 f 10 0 2 F 2 f 14 1 3 F 7 f 18 2 f 2 3 0 F 8 f 6 4 1 F 13 2 F 18 f 15 0 3 F 3 f 19 1 f 3 2 0 F 4 f 7 3 1 F 9 f 11 4 2 F 14 3 F 19 k1 0 n2 0 n2 1 n2 2 n2 3 k1 1 k1 2 k1 3 k1 4 图图 10 素因子素因子 PFA 20 点点 FFT 算法算法 再以点为例 进行如下三层因式分解20FFT 123 5 2 2NN N N 即 由于 5 与 4 互素 所以第一层采用 PFA 算法 相应 1 5N 23 2 24N 的二层映射为 123123 45 04 03nnnnn 123123 1625 04 03kkkkk 由于 2 与 2 不互素 所以第二层采用 Cooley Tukey FFT 算法 相应的二层 映射为 232323 2 01 01nnnnn 3 一维任意非基 2 FFT 算法 18 232323 2 01 01kkkkk 总的变换公式如式 26 计算流图如图 11 所示 FFT n1 k1 n2 k2 n3 k3 f 0 0 0 0 0 W40 0 0 F 0 f 4 1 1 1 1 W40 1 1 F 10 f 8 2 2 f 12 3 3 0 0 W40 0 0 F 5 f 16 4 4 1 1 W41 1 1 F 15 0 0 W40 0 0 F 16 1 1 W40 1 1 F 6 f 5 0 0 f 9 1 1 0 0 W40 0 0 F 1 f 13 2 2 1 1 W41 1 1 F 11 f 17 3 3 f 1 4 4 0 0 W40 0 0 F 12 1 1 1 1 W40 1 1 F 2 0 0 W40 0 0 F 17 f 10 0 0 1 1 W41 1 1 F 7 f 14 1 1 f 18 2 2 0 0 W40 0 0 F 8 f 2 3 3 1 1 W40 1 1 F 18 f 6 4 4 0 0 W40 0 0 F 13 1 1 W41 1 1 F 3 f 15 0 0 0 0 W40 0 0 F 4 f 19 1 1 1 1 W40 1 1 F 14 f 3 2 2 f 7 3 3 0 0 W40 0 0 F 9 f 11 4 4 1 1 W41 1 1 F 19 n23 0 n23 1 n23 2 n23 3 n3 0 n3 1 n3 0 n3 1 n3 0 n3 1 n3 0 n3 1 n3 0 n3 1 图图 11 多层分解时多层分解时 20 点点 FFT 算法算法 321 3 33 22 21 1 32321 321 111 123123 000 NNN n kn kn kn k NNNN nnn F k k kWWWf n n n W 快速傅里叶变换 FFT 算法及其应用 19 3 33 22 21 1 321 114 2421235 000 n kn kn kn k nnn WWWf n n n W 26 比较正变换和反变换的定义式可知 在正变换前加上乘的DFTIDFT1 N 运算 把换成 并交换 和 就得到反变换的算法 W 1 W f n F knk 一维任意非基 2 FFT 算法的 Visual C 源程序见附录 2 在 Matlab 中 傅 里叶变换及其逆变换也分别用 dwt 和 idwt 函数实现 4 一维 FFT 算法的应用 4 1 利用 FFT 计算连续时间信号的傅里叶变换 设是连续时间信号 并假设时 则其傅里叶变换由下式给 x t0t 0 x t 出 0 i t Xx t edt 令是一个固定的正实数 是一个固定的正整数 当 N 时 利用 FFT 算法可计算 0 1 2 1kkN X 已知一个固定的时间间隔 选择足够小 使得每一个秒的间隔TTT 内 的变化很小 则式中积分可近似为 1 nTtnT x t 1 0 nT iwt nT n Xedt x nT 1 0 1 i t tnT t nT n ex nT i 0 1 i T i nT n e ex nT i 27 假设足够大 对于所有的整数 幅值很小 则式 27 变NnN x nT 为 1 0 1 i TN i nT n e Xex nT i 28 当时 式 28 两边的值为2 k NT 4 一维 FFT 算法的应用 20 2 2 1 2 0 211 2 2 ik Nik N N ink N n kee Xex nTX k NTik NTik NT 29 其中代表抽样信号的点 最后令 则上式 X k x nx nT NDFT2 NT 变为 2 1 0 1 2 1 2 ik N e X kX kkN ik NT 30 首先用 FFT 算法求出 然后可用上式求出时的 X k0 1 2 1kN X k 应该强调的是 式 28 只是一个近似表示 计算得到的只是一个近 X 似值 通过取更小的抽样间隔 或者增加点数 可以得到更精确的值 如TN 果时 幅度谱很小 对应于奈奎斯特抽样频率 抽样间隔B X 2 s B 选择比较合适 如果已知信号只在时间区间内存在 可以通过T B 1 0tt 对时的抽样信号补零 使足够大 1 nTt x nx nT N 例例 1 利用 FFT 计算傅里叶变换 如图 12 所示的信号 102 0 tt x t 其它 其傅里叶变换为 2 cos sin 2 i Xei 利用下面的命令 可得到的 x t 近似值和准确值 图图 12 连续时间信号连续时间信号 x t N input Input N T input Input T 计算 X w 近似值 t 0 T 2 x t 1 zeros 1 N length t X fft x gamma 2 pi N T k 0 10 gamma Xapp 1 exp i k gamma T i k gamma X 计算真实值 X w w 0 05 0 05 10 快速傅里叶变换 FFT 算法及其应用 21 Xact exp i w 2 i w cos w sin w w w plot k gamma abs Xapp 1 length k o w abs Xact legend 近似值 真实值 xlabel 频率 rad s ylabel X 运行程序后输入 N 128 T 0 1 此时 得到实际的和近似的傅0 4909 里叶变换的幅度谱如图 13 所示 此时近似值已经相当准确 通过增加 NT 可以 增加更多的细节 减少 T 使得到的值更精确 再次运行程序后输入 N 512 T 0 05 此时 得到实际的和近似的傅里叶变换的幅度谱如0 2454 图 14 所示 图图 13 N 128 T 0 1 时的幅度谱时的幅度谱 4 一维 FFT 算法的应用 22 图图 14 N 512 T 0 05 时的幅度谱时的幅度谱 4 2 利用 FFT 计算离散信号的线性卷积 已知两个离散时间信号与 0 1 2 1 x nnM 0 1 2 1 y nnN 取 对和右端补零 使得L 1MN x n y n 0 1 2 1x nnMML 0 1 2 1y nnNNL 31 利用 FFT 算法可以求得和的 L 点 DFT 分别是和 利 x n y n X k Y k 用 DTFT 卷积性质 卷积等于乘积的 L 点 DFT 反变换 这 x ny n X k Y k 也可以通过 FFT 算法得到 例例 2 利用 FFT 计算线性卷积 已知 其中为单位阶跃序列 信号如图 15 所示 0 8 n x nu n u n y n 由于当时 很小 故可以取为 17 N 取 10 16n x nM126LMN 利用下面的 Matlab 命令 可得到 的卷积图形如图 15 所示 x n y n subplot 3 1 1 n 0 16 快速傅里叶变换 FFT 算法及其应用 23 x 0 8 n stem 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 信号信号 x n y n 及其卷积及其卷积 z n x n y n 利用下面的 Matlab 命令 可得到信号 x n y n 的幅度谱与相位谱如图 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 4 一维 FFT 算法的应用 24 xlabel k ylabel X k subplot 2 2 2 stem k angle X axis 0 25 1 1 xlabel k ylabel Angle X k 弧度 subplot 2 2 3 y ones 1 10 Y fft y L stem k abs Y axis 0 25 0 10 xlabel k ylabel Y k subplot 2 2 4 stem k angle Y axis 0 25 3 3 xlabel k ylabel Angle Y k 弧度 图图 16 信号信号 x n y n 的幅度谱与相位谱的幅度谱与相位谱 4 3 利用 FFT 进行离散信号压缩 利用 FFT 算法对离散信号进行压缩的步骤如下 1 通过采样将信号离散 化 2 对离散化信号进行傅里叶变换 3 对变换后的系数进行处理 将绝对 值小于某一阈值的系数置为 0 保留剩余的系数 4 利用 IFFT 算法对处理后 快速傅里叶变换 FFT 算法及其应用 25 的信号进行逆傅里叶变换 例例 3对单位区间上的下列连续信号 1 cos 4 sin 8 2 f tttt 以采样频率进行采样 将其离散化为个采样值256 s fHz 8 2 256 t nT f nf tf nTf n 0 1 2 255n 用 FFT 分解信号 对信号进行小波压缩 然后重构信号 令绝对值最小的 80 系数为 0 得到重构信号图形如图 17 a 所示 均方差为 0 0429 相对误差 为 0 0449 令绝对值最小的 90 系数为 0 得到重构信号图形如图 17 b 所示 均方差为 0 0610 相对误差为 0 0638 a 绝对值最小的 80 系数为 0 的重构信号 FFT b 绝对值最小的 90 系数为 0 的重构信号 FFT 图图 17 用用 FFT 压缩后的重构信号压缩后的重构信号 相关 Matlab 程序如下 function wc compress w r 压缩函数 compress m 输入信号数据 w 压缩率 r 输出压缩后的信号数据 if r1 error r 应该介于 0 和 1 之间 end N length w Nr floor N r ww sort abs w tol abs ww Nr 1 wc abs w tol w function unbiased variance error fftcomp t y r 利用 FFT 做离散信号压缩 输入时间 t 原信号 y 以及压缩率 r 4 一维 FFT 算法的应用 26 输出原信号和压缩后重构信号的图像 以及重构均方差和相对 l 2 误差 if r1 error r 应该介于 0 和 1 之间 end fy fft y fyc compress fy r 调用压缩函数 compress m yc ifft fyc plot t y r t yc b legend 原信号 重构信号 unbiased variance norm y yc sqrt length t error norm y yc norm y 输入以下 Matlab 命令 t 0 255 256 f t cos 4 pi t 1 2 sin 8 pi t unbiased variance error fftcomp t f 0 8 unbiased variance 0 0429 error 0 0449 如果用 Harr 尺度函数和 Harr 小波分解信号 对信号进行小波压缩 然后 重构信号 令绝对值最小的 80 系数为 0 得到重构信号图形如图 18 a 所示 均方差为 0 0584 相对误差为 0 0611 令绝对值最小的 90 系数为 0 得到重 构信号图形如图 18 b 所示 均方差为 0 1136 相对误差为 0 1190 a 绝对值最小的 80 系数为 0 的重构信号 Harr b 绝对值最小的 90 系数为 0 的重构信号 Harr 图图 18 用用 Harr 小波压缩后的重构信号小波压缩后的重构信号 相关 Matlab 程序如下 快速傅里叶变换 FFT 算法及其应用 27 function unbiased variance error daubcomp t y n r 利用 Daubechies 系列小波做离散信号压缩 输入时间 t 原信号 y 分解层数 n 以及压缩率 r 输出原信号和压缩后重构信号的图像 以及重构均方差和相对 l 2 误差 if r1 error r 应该介于 0 和 1 之间 end c l wavedec y n db1 cc compress c r 调用压缩函数 compress m yc waverec cc l db1 plot t y r t yc b legend 原信号 重构信号 unbiased variance norm y yc sqrt length t error norm y yc norm y 输入以下 Matlab 命令 t 0 255 256 f t cos 4 pi t 1 2 sin 8 pi t unbiased variance error daubcomp t f 8 0 8 unbiased variance 0 0584 error 0 0611 结论 在信号没有突变 快变化或者大致上具有周期性的信号 用 FFT 可 以处理得很好 甚至比小波还要好 4 4 利用 FFT 对离散信号进行滤波 利用 FFT 算法对信号进行滤波的步骤如下 1 通过采样将信号离散化 2 对离散化信号进行傅里叶变换 3 对变换后的系数进行处理 将绝对值大 于某一阈值的系数置为 0 保留剩余的系数 4 利用 IFFT 算法对处理后的信 号进行逆傅里叶变换 例例 4 股票价格分析 首先进入网址 Download To Spreadsheet 按钮 即可把以 Excel 表格格式存储的价格数据下载 4 一维 FFT 算法的应用 28 到本地计算机 表格从 1 列至第 6 列分别给出了从 1996 年 4 月 12 日至 2007 年 5 月 30 日的交易期里每天的开盘价 最高价 最低价 收盘价 成交量以及 趋势 数据下载完成后 需要颠倒顺序 使得最早时间的数据首先显示 然后 另存到 Matlab 所在的目录中 并重新命名为 yhoodata csv 本次分析选择开 盘价 时间是从 2007 年 1 月 1 日至 2007 年 5 月 30 日的 102 个交易日期 N 令代表一支股票的开盘价 为了便于分析 可以先从 0 1 1o nnN 中减去跃变 得到 即 o n 0 1 1x nnN 32 1 0 0 0 1 1 1 o No x no nonnN N 输入以下命令 得到的频谱如图 19 所示 x n o csvread yhoodata csv 2700 1 2700 1 2801 1 N 102 for n 1 N x n o n o 1 o N o 1 N 1 n 1 end X fft x k 0 N 1 stem k abs X axis 0 101 0 300 xlabel k ylabel X k 快速傅里叶变换 FFT 算法及其应用 29 图图 19 x n 的幅度谱的幅度谱 可以看出上图中有 5 个较强谱分量 频率 分别0 1 4 98 100k 2 k N 对应和 保留这 5 个频率分量的系数 将其他频率分量的系0 2 102 8 102 数置为 0 然后再进行逆傅里叶变换 得到滤波后的近似值 输入如下 x n Matlab 程序 得到真实值与滤波后的近似值 如图 20 所示 x n x n plot x hold on fliterX X 1 2 0 0 X 5 zeros 1 102 9 X 99 0 0 X 102 fliterx ifft fliterX plot real fliterx r axis 1 102 0 7 xlabel n ylabel x n 的值和它的近似值 legend x n 真实值 x n 近似值 4 一维 FFT 算法的应用 30 图图 20 x n 的真实值与滤波后的近似值的真实值与滤波后的近似值 从上图可以看出 滤波后的近似值既大致上保留了真实值的变化趋势 x n 而且与其十分接近 与滤波前比较 滤波后的图形要比滤波前平滑得多 再由 式 32 即可求得 33 1 0 0 0 0 1 1 1 o No ox nonnN N 输入如下 Matlab 程序 画出真实开盘价与近似开盘价的图形 如 o n o n 图 21 所示 可以看出是近似基础上的平滑值 o n o n plot o hold on for n 1 N oapp n fliterx n o 1 o N o 1 N 1 n 1 end plot oapp r axis 1 102 25 34 xlabel n ylabel o n 的值和它的近似值 legend o n 真实值 o n 近似值 快速傅里叶变换 FFT 算法及其应用 31 图图 21 o