c2java 第9篇 FFT

傅里叶分析是什么？

设多项式A(x) = a_0 + a_1*x^1 + a_2*x^2 + ... + a_{n-1}*x^{n-1},
系数向量记为a = (a_0, a_1, ..., a_{n-1}),
多项式求值向量记为y=(y_0, y_1, ..., y_{n-1}),
y_{j} = A(w_n^j),
w_n = e^{2*\pi*i/n},
称y = FT(a)为离散傅里叶变换
反过来， 跟据y求系数a也很容易，称a=IFT(y)为逆傅里叶变换。
fft-test.c/fft_ref就是按定义来求FT的，复杂度为O(N^2)

基本原理
====
设n是2的幂。把A(x)拆分为
B(x) = a_0 + a_2 * x + ... + a_{n-2} * x ^{n/2-1}
C(x) = a_1 + a_3 * x + ... + a_{n-1} * x ^{n/2-1}
则A(x) = B(x^2) + x*C(x^2); (I)
直观看B(x), B(x)的次数为n/2, 问题规模减半，那计算量f(n)是否减半呢？
我们需要计算n点：(w_n^0)^2, (w_n^1)^2, ..., (w_n^n)^2,
由折半引理:(w_n^{k+n/2})^2 = (w_n^k)^2， 第k个点等于第n/2+k个点。
得到f(n) = 2*f(n/2) + O(n).

如果按递归写法，需要给子问题开辟系数存储所需的空间，某一时刻
最大为：n + n/2 + n/4 + ... + 2 + 1 = 2*n.
实际内存要求比这个值大，因为内存管理还有额外开销。
一定需要额外空间么，可不可以原地操作呢？

迭代版本
====
递归并不一定能转化为迭代的，有时需要我们重新定义问题和输入数据。
递归中，计算顺序是假设子问题已经算好的，因此不需要特别关；而迭代中却是首要解决的问题。
观察下面这个由递归调用形成的树。
从底到上，第二层左边第一个节点(b0, b1)：
b0 = a0 + a4; b1 = a0 - a4;
因为在同一层中，没有其他地方修改或引用原来的a0, a4, 因此
我们可以直接把存储空间b0,b1 放到a0,a4。
因为我们首先是按a0,a4,a2,a6,a1,a5,a3,a7顺序计算的，因此开始时需要
把原始数据A摆成这样的顺序。
下面是把这个过程一般化：
FFT(A, bits)
arrange(A);
n = 2^bits;
for(s = 1; s <= bits; ++s){//s表示在哪一层
 m = 2^s;//当前问题的规模
 for(k = 0; k < n; k += m){//k表示从左边开始数第几个节点
  A[k,...,k+m-1] = 应用等式(I)合并
   A[k,...,k+m/2-1] 和 A[k+m/2, ..., k+m-1]
 }
}

ffmpeg 中的FFT
====
fft_template.c/fft_calc_c 的实部和虚部是放在同一个数组的两端，
比且对于n=4,n=8直接内循环展开。

SSE 3DNow向量指令优化：
x86/fft.asm
 fft_calc

应用
====
1. 多项式乘法，即卷积， 为什么卷积重要？
2. ffmpeg:wma*.c 等，DCT似乎在每个codec都有。
3. 为什么要有MDCT ？

 

代码：

/*FFT.java -- 快速傅里页变换*/
import java.util.*;

class Complex
{
	double re, im;
	public Complex(double a, double b)
	{
		re = a; im = b;
	}
	public String toString()
	{return re + "+" + im+"i";}

	public void exp(int n)
	{//e^{i*2\pi/n}
		re = Math.cos(2*Math.PI/n);
		im = Math.sin(2*Math.PI/n);
	}
	public void multiply(Complex b)
	{//a *= b
		double bre = b.re, bim = b.im,
			   are = re, aim = im;
		re = are * bre - aim * bim;
		im = are * bim + aim * bre;
	}
	public void add(Complex b)
	{
		double bre = b.re, bim = b.im,
			   are = re, aim = im;
		re = are + bre;	
		im = aim + bim;
	}
	public void sub(Complex b)
	{
		double bre = b.re, bim = b.im,
			   are = re, aim = im;
		re = are - bre;	
		im = aim - bim;
	}
	public void assign(Complex b)
	{
		re = b.re;
		im = b.im;
	}	
}

public class FFT
{
	static void fft(Complex[] a, int bits)
	{
		int s, k, j, m, n = 1 << bits;
		Complex w0 = new Complex(.0, .0),
				w  = new Complex(.0, .0),
				u  = new Complex(.0, .0),
				t  = new Complex(.0, .0),
				one = new Complex(1., .0);

		for(s = 1; s <= bits; ++s){
			m = 1 << s;
			w0.exp(m);
			//System.out.printf("w0=%s ", w0);
			for(k = 0; k < n; k += m){
				w.assign(one);
				for(j = 0; j < m/2; ++j){
					u.assign(a[k+j]); //u = a[k+j];
					t.assign(a[k+m/2+j]); t.multiply(w); //t = w * a[k+m/2+j];
					a[k+j].assign(u); a[k+j].add(t); //a[k+j] = u + t;
					a[k+m/2+j].assign(u); a[k+m/2+j].sub(t); //a[k+m/2+j] = u - t;
					w.multiply(w0);
				}
			}
		}	
	}

	static int bitReverse(int n, int bits)
	{
		if(bits <= 1)return n&1;
		return ((n&1)<<(bits-1)) | bitReverse(n>>1, bits-1);
	}
	
	static void init(Complex[] a, Complex[] b, int bits)
	{
		int i, t, n = 1 << bits;	
		for(i = 0; i < n; ++i){
			t = bitReverse(i, bits);
			b[t].assign(a[i]); //a[rev[i]] = a[i];
			//System.out.printf("%d  ", t);
		}
	}

	static void dump(String prefix, Complex[] a)
	{
		int i;
		System.out.printf("%s", prefix);
		for(i = 0; i < a.length; ++i){
			System.out.printf("%.16f+%.16fi, ", a[i].re, a[i].im);
		}
		System.out.printf("%n");
	}

	static void fftRef(Complex[] b, Complex[] a, int bits)
	{
		int i, j, n = 1 << bits;
		Complex x = new Complex(.0, .0),
				sum = new Complex(.0, .0),
				w = new Complex(.0, .0),
				one = new Complex(1.0, .0);	
		x.assign(one);
		w.exp(n);
		for(j = 0; j < n; ++j){
		sum.assign(a[n-1]);
		//下面使用horner方法求多项式值，因为w^j可以直接求，也可以
		//直接从a0开始计算。
		for(i = n - 2; i >= 0; --i){
			sum.multiply(x); sum.add(a[i]); //sum = sum * x +  a[i];
		}
		b[j].assign(sum);
		x.multiply(w);
		//System.out.printf("x=%s ", x);
		}
	}

	static void checkDiff(Complex[] a, Complex[] b)
	{
		int i, n = a.length;
		double max = 0, err = 0., sum = 0.;
		for(i = 0; i < n; ++i){
			err = Math.abs(a[i].re - b[i].re);
			sum += err * err;
			if(err > max)max = err;
			err = Math.abs(a[i].im - b[i].im);
			sum += err * err;	
			if(err > max)max = err;
		}
		System.out.printf("err max=%g, avg=%g%n", max,
				Math.sqrt(sum/n));
	}	

	public static void main(String[] arg){
		int i, j, bits = 16; //3;
		int n = 1 << bits;
		System.out.printf("bits = %d%n", bits);
		/*for(i = 0; i < (1<
   
   
    
    ludi@ubun:~/java