二維快速傅里葉變換的java實現

圖像處理與模式識別課作業，沒學過信號與系統（哭暈）。
  惡補了兩天岡薩里斯的書，寫一下實現原理及過程
  看了網絡上很多版本的概念講解，道理我都懂，但是在將算法遷移到傅里葉變換的實現上時，出現了一些問題。接下來我會簡要介紹快速傅里葉變換是如何加速的，着重寫如何將加速算法應用到傅里葉變換上。

二維快速傅里葉變換原理介紹

1.1普通的二維傅里葉變換

​ 二維傅里葉變換的公式如下：

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

對上述表達式調整和分順序，先對x求和，再對y求和。我們可以將二維離散傅里葉變換公式轉化為兩次一維離散傅里葉變換。一維離散傅里葉變換公式如下：

\[F(u)=\sum^{N-1}_{x=0}f(x)e^{{{-j2\pi}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_n\quad k\in(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{2\pi}k/n}=cos(2\pi k/n)+i*sin(2\pi k/n) \]

從虛數相乘等於在圓上的旋轉，我們可以看出它確實是1的n次單位根的集合。它有幾個特殊的性質

\[w^{k+n/2}_n=e^{i2\pi (k+n/2)/n}=e^{i2\pi k/n}*e^{i\pi}=e^{i2\pi k/n}*(cos(\pi)+i*sin(\pi))=-e^{i2\pi k/n}=-w^k_n\\ w^{k+n}_n=e^{i2\pi k/n}*e^{i2\pi}=e^{i2\pi 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^{{{-j2\pi}ux}/N}\quad u = 0,1,2...N-1 \]

是不是看到了\(e^{{{-j2\pi}ux}/N}\)? 注意：我們接下來使用的\(w^k_n\)=\(e^{-j2\pi 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)] \]

第三次到遞歸終點，我們看一下上訴編號的二進制碼

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

將二進制反序

\[000\quad 001\quad 010\quad 011\quad 100\quad 101\quad 110\quad 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";
    }
    
}

填坑完畢，轉發需注明出處。