分治法：快速傅里叶变换算法

分治法：快速傅里叶变换算法

学院：网研院 姓名：xxx 学号：xxx

一、 分治法原理

分治法是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题……直到最后子问题可以简单的直接求解，原问题的解即子问题的解的合并。

分治法的精髓：

◆ 分--将问题分解为规模更小的子问题； ◆ 治--将这些规模更小的子问题逐个击破；

◆ 合--将已解决的子问题合并，最终得出“母”问题的解； 二、 快速傅里叶变换（FFT ）简介

快速傅里叶变换（Fast Fourier Transform, FFT ），是离散傅里叶变换的快速算法，也可用于计算离散傅里叶变换的逆变换。它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。快速傅里叶变换有广泛的应用，如数字信号处理、计算大整数乘法、求解偏微分方程等等。

序列

的离散傅里叶变换公式为：

令 则上式可写为：

从算法分析角度：

于是：分别考虑对其奇数项和偶数项作傅氏变换： 则

从而，可以将对N 个量的傅氏变换变成为对两个规模更小的序列的变换。这

样，将变换的量大大减小。

快速傅里叶变换是分治法的一种具体实现。

[][]()

∑-=-=10

2N n N jnk e n x k X π

1,,0]},[{-=N i i x N

j N e W π

2-=∑-==10

][][N n kn

N

W n x k X k

n N

N

n N n nk N N n kn N W n x W n x W n x k X )12(120

12

210

]12[]2[][][+-=-=-=∑∑∑?++?==nk

N N n N n nk N E W n x W n x k X 2

/12

12

2]2[]2[][∑∑-=-=?=?=kn

N N n N n nk N O W n x W n x k X 2

/120

120

2]12[]12[][∑∑-=-=?+=?+=]

[][][][1

k X W k X W n x k X O k N E N n kn N ?+==∑-=

三、 快速傅里叶变换算法FFT

1. 算法

● 输入采样值；

● 对采用值进行补0操作，使采样值的个数是2的幂次；

● 对补0后的序列进行重排，重排的原则是按照序列的下标奇偶性排

序，先偶后奇，分成两个子序列，然后对子序列继续重排。如

1,2,3,4,5,6,7->1,3,5,7;2,4,6,8->1,5;3,7;2,6;4,8->1,5,3,7,2,6,4,8； ● 对补0重排后的序列用三层嵌套的循环来完成M=log 2N （阶数）次迭

代，三层循环的功能是：最里的一层循环完成相同W N P 的蝶形运算，

中间的一层循环完成因子W N P

的变化，而最外的一层循环则是完成M 次迭代过程。 2. 算法复杂度

令

，则得 从而快速傅氏变换的复杂性度为

3. 可能的改进

本次实验使用的是最基本的基2FFT 算法，根据维基百科，有如下比较出名的FFT 算法，在性能上应该都优于本实现。

Cooley-Tukey 算法是最常见的FFT 算法。这一方法以分治法为策略递归地将长度为的DFT 分解为长度为

的

个较短序列的DFT ，

以及与

个旋转因子的复数乘法。

在Cooley-Tukey 算法之外还有其他DFT 的快速算法。对于长度

且

与

互质的序列，可以采用基于中国剩余定理的互质因

子算法将N 长序列的DFT 分解为两个子序列的DFT 。与 Cooley-Tukey 算法不同的是，互质因子算法不需要旋转因子。

Rader-Brenner 算法是类似于 Cooley-Tukey 算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey 所取代，与

?????

=+

=0

)2(2

)2(2)(T N N T N T m N 2=12)1()2(--=m m m T )log (N N O

Rader-Brenner算法相比，有一样多的乘法量，却有较少的加法量，且不牺牲数值的准确性。

Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。

(Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字)。

尤其是Bruun算法，把FFT看成是，并把它分解成与

的形式。

另一个从多项式观点的快速傅立叶变换法是Winograd 算法。此算法把

分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。

这样只需要很少的乘法量(如果有需要的话)，所以winograd是可以得到最少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。

更精确地说，winograd算法让DFT可以用点的DFT来简化，但减少乘法

量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。

Rader算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋折积来表示原本的DFT，如此就可利用折积用一对基本的FFT来计算 DFT。

另一个prime-size的FFT算法为chirp-Z算法。此法也是将DFT用折积来表示，此法与Rader算法相比，能运用在更一般的转换上，其转换的基础为Z转换。

四、算法实现框架

本次实验使用的语言是java，只有一个类FFT.java，类有一个属性xConv

用来存储补0和重排后的序列；FFT()构造方法对序列进行补0操作，然后调用

i2Sort()方法对补0后序列进行重排，最后调用myFFT()方法进行快速傅里叶变换；i2Sort()法对补0后序列进行重排；myFFT()方法对补0重排后序列进行快速傅里叶变换操作；main()方法是程序的入口，用于输入数据并调用构造方法FFT()。如下图所示。

i2Sort()方法对补0后序列进行重排，如1,2,3,4,5,6,7->1,5,3,7,2,6,4,8。

myFFT()方法对补0重排后序列进行快速傅里叶变换操作，方法的原理是用三层嵌套的循环来完成M=log

2

N（阶数）次迭代，三层循环的功能是：最里的一

层循环完成相同W

N P的蝶形运算，中间的一层循环完成因子W

N

P的变化，而最外

的一层循环则是完成M次迭代过程。

FFT()构造方法对序列进行补0操作，然后调用i2Sort()方法对补0后序列进行重排，最后调用myFFT()方法进行快速傅里叶变换；

main()方法是程序的入口，用于输入数据并调用构造方法FFT()。

五、总结

本科的时候曾经系统学习过傅里叶变换的知识，考研的时候也复习了该方面知识，在理解题目上没有遇到太大的困难。但由于从未尝试编程实现FFT，因此在本次实验的过程中还是遇到了许多问题。经过查看课件及网上查找资料，我对如何编码实现FFT有了一个较好的理解，知道了算法主要应该包括补0、排序和蝶形迭代3步操作。通过查找资料参考别人的代码，并结合自己的理解，我使用了java语言完成了该实验。由于时间较为紧张，加上期末考临近等原因，该实验完成的比较仓促，仍有许多不足的地方，但总的来说，通过本次试验，我再次复习了FFT的原理及了解了分治法的原理，同时锻炼了自己动手编写代码的能力，收获多多。

附录

程序运行示例：

源程序：

import java.util.Scanner;

public class FFT {

private double[] xConv;// 对x[n]进行二进制倒序排列的结果

public FFT(double[] x) throws Exception {

int n = x.length;

if (n < 1)

throw new Exception("采样序列的个数不能小于1！");

int m = (int) Math.ceil(Math.log(n) / Math.log(2));

int j = (int) Math.pow(2, m);

xConv = new double[j]; // 初始化xConv

for (int i = 0; i < n; i++)

xConv[i] = x[i];

for (int i = n; i < j; i++)

xConv[i] = 0;

System.out.println("x[n]共有" + x.length + "个采样值" + '(' + m + "阶)"

+ ",补零个数为" + (j - n));

System.out.println("原x[n]数组：");

for (int i = 0; i < x.length; i++) {

System.out.print(String.format("%6.2f", x[i]) + " ");

}

System.out.println("");

System.out.println("补零之后的x[n]数组：");

for (int i = 0; i < xConv.length; i++) {

System.out.print(String.format("%6.2f", xConv[i]) + " ");

}

System.out.println("");

i2Sort(xConv, m);

myFFT(xConv, m);

}

private void i2Sort(double[] xConv2, int m) {

int[] index = new int[xConv2.length]; // index数组用于，倒序索引

int[] bits = new int[m];

double[] temp = new double[xConv2.length];

for (int i = 0; i < xConv2.length; i++)

// xConv2的原序映像

temp[i] = xConv2[i];

for (int i = 0; i < index.length; i++) {

index[i] = i; // 第i个位置，倒序前的值为i

for (int j = 0; j < m; j++) {

bits[j] = index[i] - index[i] / 2 * 2; // 提取index[i]的第j位二进制的值

index[i] /= 2;

}

index[i] = 0; // 清零第i个位置的值

for (int j = m, power = 1; j > 0; j--) {

index[i] += bits[j - 1] * power; // 第i个位置，倒序后的位置

power *= 2;

}

}

System.out.println("二进制排序之后的x[n]数组：");

for (int i = 0; i < xConv2.length; i++) {

// 倒序实现

xConv2[i] = temp[index[i]];

System.out.print(String.format("%6.2f", xConv2[i]) + " ");

}

System.out.println("");

}

private void myFFT(double[] xConv2, int m) {

int divBy; // divBy等分

double[] Xr, Xi, Wr, Wi; // 分别表示：FFT结果的实部和虚部、旋转因子的实部和虚部

double[] tempXr, tempXi; // 蝶形结果暂存器

int n = xConv2.length;

double pi = Math.PI;

divBy = 1;

Xr = new double[n];

Xi = new double[n];

tempXr = new double[n];

tempXi = new double[n];

Wr = new double[n / 2];

Wi = new double[n / 2];

for (int i = 0; i < n; i++) { // 初始化Xr、Xi，之所以这样初始化，是为了方便下面的蝶形结果暂存

Xr[i] = xConv2[i];

Xi[i] = 0;

}

for (int i = 0; i < m; i++) { // 共需要进行m次蝶形计算

divBy *= 2;

for (int k = 0; k < divBy / 2; k++) { // 旋转因子赋值

Wr[k] = Math.cos(k * 2 * pi / divBy);

Wi[k] = -Math.sin(k * 2 * pi / divBy);

}

for (int j = 0; j < n; j++) { // 蝶形结果暂存

tempXr[j] = Xr[j];

tempXi[j] = Xi[j];

}

for (int k = 0; k < n / divBy; k++) { // 蝶形运算：每一轮蝶形运算，都有n/2对的蝴蝶参与；n/2分为n/divBy组，每组divBy/2个。

int wIndex = 0; // 旋转因子下标索引

for (int j = k * divBy; j < k * divBy + divBy / 2; j++) {

double X1 = tempXr[j + divBy / 2] * Wr[wIndex]

- tempXi[j + divBy / 2] * Wi[wIndex];

double X2 = tempXi[j + divBy / 2] * Wr[wIndex]

+ tempXr[j + divBy / 2] * Wi[wIndex];

Xr[j] = tempXr[j] + X1;

Xi[j] = tempXi[j] + X2;

Xr[j + divBy / 2] = tempXr[j] - X1; // 蝶形对两成员距离相差divBy/2

Xi[j + divBy / 2] = tempXi[j] - X2;

wIndex++;

}

}

}

System.out.println("FFT结果：");

for (int i = 0; i < n; i++)

// FFT结果显示

System.out.println(String.format("%6.2f", Xr[i])

+ (Xi[i] < 0 ? " - " : " + ")

+ String.format("%6.2f", Math.abs(Xi[i])) + " * j");

}

public static void main(String[] args) throws Exception {

Scanner in = new Scanner(System.in);

while (true) {

System.out.println("******基2时域抽取FFT算法开始******");

System.out.println("---请输入采样点个数(正整数),回车结束---");

int n = in.nextInt();

System.out.println("---请输入" + n + "个采样值(支持浮点型，空格隔开),回车结束---");

double[] x = new double[n];

for (int i = 0; i < n; i++) {

x[i] = in.nextDouble();

}

System.out.println("---基2时域抽取FFT算法处理结果---");

new FFT(x);

System.out.println("");

}

}

}