二维快速傅里叶变换的java实现

图像处理与模式识别课作业，没学过信号与系统（哭晕）。
  恶补了两天冈萨里斯的书，写一下实现原理及过程
  看了网络上很多版本的概念讲解，道理我都懂，但是在将算法迁移到傅里叶变换的实现上时，出现了一些问题。接下来我会简要介绍快速傅里叶变换是如何加速的，着重写如何将加速算法应用到傅里叶变换上。

二维快速傅里叶变换原理介绍

1.1普通的二维傅里叶变换

​ 二维傅里叶变换的公式如下：

[F(u, v) = sum^{M-1}_{x=0}sum^{N-1}_{y=0}f(x, y)e^{-j2pi({ux}/M+{vy}/N)} ]

对上述表达式调整和分顺序，先对x求和，再对y求和。我们可以将二维离散傅里叶变换公式转化为两次一维离散傅里叶变换。一维离散傅里叶变换公式如下：

[F(u)=sum^{N-1}_{x=0}f(x)e^{{{-j2pi}ux}/N}quad u = 0,1,2...N-1 ]

​ 从变换公式中我们可以看出来，由于我们要求N个(F(u))，每个(F(u))需要进行N次运算。使用这样算法的时间复杂度是(O(n^2))的。这样的算法在N特别大时，会比较慢。在这里引入快速傅里叶变换的加速算法。

1.2快速傅里叶变换原理

​ 相信大家在看过别的博客内容后已经掌握了快速傅里叶变换的原理。我在这里就简要再提一下。快速傅里叶变换其中一个重要的依据就是函数的点值对表示。其实我们对离散函数求得的傅里叶变换的结果就是一个点值对表示（毕竟我们没有求出傅里叶变换函数的系数）。

​ 一个n次多项式函数，注意：这里的n是2的整数幂，如果一个函数的最高次没有达到2的整数幂，我们可以认为此项系数为0，不影响函数的值

[A(x) = a_0 + a_1x+a_2x^2+...+a_{n-1}x^{n-1} ]

在奇偶项合并后，我们可以得到

[A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})\ A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})\ A(x)=A_{even}(x^2)+xA_{odd}(x^2) ]

​ 我们将(w^k_nquad kin(0, 1, 2,...,n-1))带入A(x)。可以得到A(x)的点值对表示。但是这种点值对的计算方法，仍然是(O(n^2))的。如何将计算加速呢？这里就用到我们取巧的(w^k_n)上了。

​ 简要介绍一下(w^k_n)，它是正整数1的n次单位根。具体值为：

[w^k_n=e^{i{2pi}k/n}=cos(2pi k/n)+i*sin(2pi k/n) ]

从虚数相乘等于在圆上的旋转，我们可以看出它确实是1的n次单位根的集合。它有几个特殊的性质

[w^{k+n/2}_n=e^{i2pi (k+n/2)/n}=e^{i2pi k/n}*e^{ipi}=e^{i2pi k/n}*(cos(pi)+i*sin(pi))=-e^{i2pi k/n}=-w^k_n\ w^{k+n}_n=e^{i2pi k/n}*e^{i2pi}=e^{i2pi k/n}=w^k_n ]

所以我们在计算点值对的时候，有这样一种关系：

[A(w^k_n)=A_{even}((w^k_n)^2)+w^k_n*A_{odd}((w^k_n)^2) -->A(w^k_n)=A_{even}(w^{2k}_n)+w^k_n*A_{odd}(w^{2k}_n)\ A(w^{k+n/2}_n)=A_{even}((w^{k+n/2}_n)^2)+w^{k+n/2}_n*A_{odd}((w^{k+n/2}_n)^2)--> A(w^{k+n/2}_n)=A_{even}(w^{2k+n}_n)-w^k_n*A_{odd}(w^{2k+n}_n)-->A(w^k_n)=A_{even}(w^{2k}_n)-w^k_n*A_{odd}(w^{2k}_n) ]

整理一下，结果如下

[A(w^k_n)=A_{even}(w^{2k}_n)+w^k_n*A_{odd}(w^{2k}_n)\ A(w^{k+n/2}_n)=A_{even}(w^{2k}_n)-w^k_n*A_{odd}(w^{2k}_n) ]

则我们求出(A_{even}(w^{2k}_n))和(A_{odd}(w^{2k}_n))就可以得到两个值。根据我们算法课学的分治+递归的思想。递归层数为logn次，每层计算时间复杂度为(O(n))，则总的时间复杂度为(O(nlogn))。这样我们就做到了快速计算。

1.3将算法应用到傅里叶变换

​ 之前我们大致了解了傅里叶变换的原理，那么如何将这个算法加速我们傅里叶变换的计算呢？一维离散傅里叶变换的公式为：

[F(u)=sum^{N-1}_{x=0}f(x)e^{{{-j2pi}ux}/N}quad u = 0,1,2...N-1 ]

是不是看到了(e^{{{-j2pi}ux}/N})? 注意：我们接下来使用的(w^k_n)=(e^{-j2pi k/n})。可以证明，这并不影响我们之前介绍的性质。

​ 则我们的傅里叶变换公式可以写成

[F(u)=sum^{N-1}_{x=0}f(x)w^{ux}_n ]

那么我们应用加速的地方就出现了

[F(u)=sum_{t=0}^{frac{n}{2}-1}f_{even}(t)w^{2ut}_n+w^u_n * sum_{t=0}^{frac{n}{2}-1}f_{odd}(t)w^{2ut}_n\ F(u+n/2)=sum_{t=0}^{frac{n}{2}-1}f_{even}(t)w^{2ut}_n-w^u_n * sum_{t=0}^{frac{n}{2}-1}f_{odd}(t)w^{2ut}_n ]

是不是和我们之前看到的算法原理的地方很像？是的，我们快速傅里叶变换就是利用这个公式做的。利用这个公式，一次就可以计算出两个值，加速了傅里叶变换计算的层数。同样的，这个算法共有logn层，每层的时间复杂度为(O(n))。整体算法的时间复杂度为(O(nlogn))。

这里还有一个trick，即蝶形算法。我们在这里举个蝶形算法的例子。假设有下列一组数

[f(0)quad f(1)quad f(2)quad f(3)quad f(4)quad f(5)quad f(6)quad f(7) ]

第一次划分奇偶后，我们得到

[[f(0)quad f(2)quad f(4)quad f(6)]quad [f(1)quad f(3)quad f(5)quad f(7)] ]

第二次划分奇偶后，结果

[[f(0)quad f(4)]quad [f(2)quad f(6)]quad [f(1)quad f(5)]quad [f(3)quad f(7)] ]

第三次到递归终点，我们看一下上诉编号的二进制码

[000quad 100quad 010quad 110quad 001quad 101quad 011quad 111 ]

将二进制反序

[000quad 001quad 010quad 011quad 100quad 101quad 110quad 111 ]

是不是就成了一个自然数列？我们可以根据这个蝶形的方式，先将(f(x))排列好，然后以循环的方式计算最终的结果。

核心java代码

import java.util.Comparator;
import java.util.ArrayList;

public class FFT {
    private ArrayList<sortData> data;
    private ArrayList<Complex> buff;

    public FFT(ArrayList<Complex> list, int pow){
        data = new ArrayList<>();
        for(int i = 0; i < list.size(); i++){
            sortData tmp = new sortData(list.get(i), reverse(i, pow));
            data.add(tmp);
        }
        Comparator<sortData> comparator = (o1, o2) -> Integer.compare(o2.getRevNum(), o1.getRevNum());
        data.sort(comparator);
        buff = new ArrayList<>();
    }
    
    public ArrayList<Complex> trans(){
        recursion(data.size());
        return buff;
    }
    
    private void recursion(int half){
        for (sortData tmpData : data) {
            // init the array
            Complex tmp = new Complex(1.0);
            tmp.multiply(tmpData.getValue());
            buff.add(tmp);
        }
        int count = 2;
        while(half >= count) {
            for (int i = 0; i < data.size(); i = i + count) {
                for (int j = i; j < i + count / 2; j++) {
                    Complex tmp1 = buff.get(j);
                    Complex tmp2 = buff.get(j + count / 2);
                    Complex tmpComplex = tmp1.complexClone();
                    Complex w = root(count, 1, j - i);
                    tmp2.multiply(w);
                    tmpComplex.add(tmp2);
                    tmp1.subtract(tmp2);
                    buff.set(j, tmpComplex);
                    buff.set(j + count / 2, tmp1);
                }
            }
            count = count * 2;
        }
    }

    private int reverse(int value, int pow){
        int res = Integer.reverse(value);
        res = res >> (Integer.SIZE - pow);
        res = res & ((1 << pow) - 1);
        return res;
    }

    private Complex root(int k, int x, int u){
        double real = Math.cos(Math.PI * 2 * u * x / k);
        double imaginary = -1 * Math.sin(Math.PI * 2 * u * x / k);
        return new Complex(real, imaginary);
    }
}

class sortData{
    private Complex value;
    private int revNum;
    public sortData(Complex o, int revNum){
        value = o;
        this.revNum = revNum;
    }
    
    public int getRevNum(){
        return revNum;
    }
    
    
    public Complex getValue(){
        return value;
    }
    
}

recursion函数是计算构造方法中传入的ArrayList的傅里叶变换，计算结果存在buff中。

reverse函数是利用蝶形算法的原理，计算出每个采样点二进制的反转结果。

root函数是计算1的n次单位结果(w^{ux}_k)。

Complex是一个复数类，具体内容如下。

public class Complex {
    private double real;
    private double imaginary;
    
    public Complex(){
        real = 0;
        imaginary = 0;
    }
    
    public Complex(double num){
        real = num;
        imaginary = 0;
    }
    
    public Complex(double num1, double num2){
        real = num1;
        imaginary = num2;
    }
    
    public Complex complexClone(){
        return new Complex(real, imaginary);
    }

    public void add(Complex o){
        real = real + o.real;
        imaginary = imaginary + o.imaginary;
    }
    
    public void subtract(Complex o){
        real = real - o.real;
        imaginary = imaginary - o.imaginary;
    }
    
    public void multiply(Complex o){
        double tmp = real * o.real - imaginary * o.imaginary;
        imaginary = real * o.imaginary + imaginary * o.real;
        real = tmp;
    }
    
    public void absLog(){
        real = Math.sqrt(real * real + imaginary * imaginary)/50;
        if(real > 255){
            real = 255;
        }
        imaginary = 0;
    }
    
    public int getReal(){
        return (int)real;
    }
    
    @Override
    public String toString(){
        if (imaginary == 0) return real + "";
        if (real == 0) return imaginary + "i";
        if (imaginary < 0) return real + " - " + (-imaginary) + "i";
        return real + " + " + imaginary + "i";
    }
    
}

填坑完毕，转发需注明出处。