# 【C++】Image FFT and IFFT

## 環境

* OS: Windows 10
* IDE: Visual Studio 2017
* 語言與框架： C++/CLI，Windows form

## 結果

* 可開啟Windows File Dialog選擇圖片。
* 
* 快速傅立葉轉換(FFT)。
* 
* 逆向快速傅立葉轉換(IFFT)。
* 

## 實作講解

### FFT(Fast Fourier Transform)

- <u>詳細的原理、證明這邊不會細講，只會提概念。</u>

#### 1D

* 我們先來想**一維**的FFT怎麼完成。
* FFT和DFT相比，用了兩個概念來加速：
* DP動態規劃。
* bit-reversal重新排序。
* 我們知道DFT必須得到整個傅立葉轉換矩陣，和原輸入相乘之後才能得到結果；FFT也是一樣，只是我們加速得到矩陣的過程。
* 每個一維都要一個矩陣，展成二維之後N*M個矩陣(行和列都需要做轉換，後面會提到)，因此這部分的化簡會提升相當大的速度。
* FFT可以利用divide and conquer，將運算從小推到大(NlogN)。
* 然而這樣會讓陣列的索引變很亂，程式碼不易撰寫。此時利用bit-reversal預先排列，就可以讓DP需要的部分排好，不需要額外做處理。

//bit-reversal重新排序
for (int i = 1, j = 0; i < N; ++i)
{
    for (int k = N >> 1; !((j ^= k)&k); k >>= 1);
    if (i > j) swap(x[i], x[j]);
}

* FFT作完之後還要將圖像做一個大小的縮放，不然值會太大。

//scale
for (int i = 0; i < N; i++)
{
    x[i] /= sqrt(N);
}

* 合併之後就是以下程式碼。

//N為Size大小，x是輸入的一維資料。
void FFT(int N,std::complex<double>* x)
{
    /* bit-reversal permutation */
    for (int i = 1, j = 0; i < N; ++i)
    {
        for (int k = N >> 1; !((j ^= k)&k); k >>= 1);
        if (i > j) swap(x[i], x[j]);
    }

    /* dynamic programming */
    for (int k = 2; k <= N; k <<= 1)
    {
        float theta = -2.0 * 3.14159 / k;
        std::complex<float> delta_w(cos(theta), sin(theta));

        // 每k個做一次FFT
        for (int j = 0; j < N; j += k)
        {
            // 前k/2個與後k/2的三角函數值恰好對稱，
            // 因此兩兩對稱的一起做。
            std::complex<double> w(1, 0);
            for (int i = j; i < j + k / 2; i++)
            {
                std::complex<double> a = x[i];
                std::complex<double> b = x[i + k / 2] * w;
                x[i] = a + b;
                x[i + k / 2] = a - b;
                w *= delta_w;
            }
        }
    }
    //scale
    for (int i = 0; i < N; i++)
    {
        x[i] /= sqrt(N);
    }
}

#### 2D

* 學會1D之後，我們就可以用1D擴展至2D。
* 概念是將每一橫排都做1D的FFT，再將**作完的結果**的每一直排做1D的FFT。
* 特別注意，不是原圖的橫排結果+原圖的直排結果，直排的資料來源要用橫排的結果。

/*///////////////////////////
// Varables Info
// M = image height
// N = image width
// InputImage** = image source
// FreqReal** = 頻譜的實數
// FreqImag** = 頻譜的虛數
///////////////////////////*/

//each col
std::complex<double> **temp = new std::complex<double>*[M];
for (int i = 0; i < M; i++)
{
    temp[i] = new std::complex<double>[N];
}

for (int i = 0; i < M; i++)
{
    //將一維的資料抽出來
    std::complex<double> *x = new std::complex<double>[N];
    for (int j = 0; j < N; j++)
    {
        //資料來源用原圖
        x[j].real(InputImage[i][j]);
        x[j].imag(0);
    }
    FFT(N,x);
    for (int j = 0; j < N; j++)
    {
        //將資料暫存到一個空間
        temp[i][j] = x[j];
    }
    delete[] x;
}
//each row
for (int i = 0; i < N; i++)
{
    //將一維的資料抽出來
    std::complex<double> *x = new std::complex<double>[M];
    for (int j = 0; j < M; j++)
    {
        //資料來源用剛剛的結果
        x[j] = temp[j][i];
    }
    FFT(M,x);
    for (int j = 0; j < M; j++)
    {
        FreqReal[i][j] = x[j].real();
        FreqImag[i][j] = x[j].imag();
    }
    delete[] x;
}

* 得到實數和虛數之後，再將兩者混合到輸出即可。

//實數虛數結合
for (int i = 0; i < M; i++)
{
    for (int j = 0; j < N; j++)
    {
        OutputImage[i][j] = sqrt(pow(FreqReal[i][j], 2) + pow(FreqImag[i][j], 2));
    }
}

### IFFT(Inverse Fast Fourier Transform)

* 反向傅立葉和正向的非常像，這代表我們幾乎不用重寫任何code。
* 如果去觀察正向和反向所需要的傅立葉矩陣的話，你會發現它們之間只差了一個共軛的關係，也就是說，你在虛數部份加上負號即可。

float theta = 2.0 * 3.14159 / k;

## 其他資源
* [演算法筆記 Wave](http://www.csie.ntnu.edu.tw/~u91029/Wave.html#4)
* 我code主要參考的網頁，裡面還有許多圖片可以更容易理解概念。
* [Paul Bourke.net](http://paulbourke.net/miscellaneous/dft/)
* 查資料的時候發現的網頁，1993年的歷史物。雖然Code比較雜亂，不過對1D轉2D的部分比較詳細。
* [我的GitHub](https://github.com/Zero871015/FourierTransform)
* 內含該專案完整程式碼。
* 如要服用最好看一下裡面的ReadMe，不然會踩雷。

## Info

* [我的網頁](https://zero871015.github.io/)
* [name=Zero871015][time=Sat, Jun 8, 2019 11:12 AM]

###### tags: `FFT` `C++`