c語言數字圖像處理（六）：二維離散傅里葉變換

基礎知識

復數表示

C = R + jI

極坐標：C = |C|(cosθ + jsinθ)

歐拉公式：C = |C|e^jθ

有關更多的時域與復頻域的知識可以學習復變函數與積分變換，本篇文章只給出DFT公式，性質，以及實現方法

二維離散傅里葉變換(DFT)

其中f(x,y)為原圖像，F(u,v)為傅里葉變換以后的結果，根據歐拉公式可得，每個F(u,v)值都為復數，由實部和虛部組成

代碼示例

 

void dft(short** in_array, double** re_array, double** im_array, long height, long width)
{
    double re, im, temp;

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            re = 0;
            im = 0;

            for (int x = 0; x < height; x++){
                for (int y = 0; y < width; y++){
                    temp = (double)i * x / (double)height + 
                           (double)j * y / (double)width;
                    re += in_array[x][y] * cos(-2 * pi * temp);
                    im += in_array[x][y] * sin(-2 * pi * temp);
                }
            }
            
            re_array[i][j] = re;
            im_array[i][j] = im;
        }
    }
    printf("dft done\n");
}

 

傅里葉譜

相角

功率譜

傅里葉變換頻譜圖

 

對於上面得兩幅圖案，在區間[0, M-1]中，變換數據由兩個在點M/2處碰面的背靠背的半個周期組成

針對顯示和濾波的目的，在該區間中有一個完整的變換周期更加方便，因為完整周期中數據是連續的

我們希望得到上圖所示的圖案

傅里葉變換的平移性質

因此對每個f(x, y)項乘以(-1)^x+y可達目的

代碼示例

void fre_spectrum(short **in_array, short **out_array, long height, long width)
{
    double re, im, temp;
    int move;

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            re = 0;
            im = 0;

            for (int x = 0; x < height; x++){
                for (int y = 0; y < width; y++){
                    temp = (double)i * x / (double)height + 
                           (double)j * y / (double)width;
                    move = (x + y) % 2 == 0 ? 1 : -1;
                    re += in_array[x][y] * cos(-2 * pi * temp) * move;
                    im += in_array[x][y] * sin(-2 * pi * temp) * move;
                }
            }
            
            out_array[i][j] = (short)(sqrt(re*re + im*im) / sqrt(width*height));
            if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
            else if (out_array[i][j] < 0)
                out_array[i][j] = 0;
 27         }
    }
}

執行結果

 

旋轉性質

 

即f(x, y)旋轉一個角度，F(u, v)旋轉相同的角度

 二維離散傅里葉反變換

代碼示例

void idft(double** re_array, double** im_array, short** out_array, long height, long width)
{
    double real, temp;

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            real = 0;

            for (int x = 0; x < height; x++){
                for (int y = 0; y < width; y++){
                    temp = (double)i * x / (double)height + 
                           (double)j * y / (double)width;

                    real += re_array[x][y] * cos(2 * pi * temp) -
                            im_array[x][y] * sin(2 * pi * temp);
                }
            }
            
            out_array[i][j] = (short)(real / sqrt(width*height));
            if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
            else if (out_array[i][j] < 0)
                out_array[i][j] = 0;
        }
    }
    printf("idft done\n");
}

經驗證，圖像經傅里葉變換，然后再反變換以后可恢復原圖

 

 

改進

本篇文章只是按照二維離散傅里葉變換公式進行了實現，在測試的過程中發現，執行速度真的是非常慢，算法時間復雜度O(n⁴)，等以后有時間再對這段代碼進行優化。