基礎知識
復數表示
C = R + jI
極坐標:C = |C|(cosθ + jsinθ)
歐拉公式:C = |C|ejθ
有關更多的時域與復頻域的知識可以學習復變函數與積分變換,本篇文章只給出DFT公式,性質,以及實現方法
二維離散傅里葉變換(DFT)
其中f(x,y)為原圖像,F(u,v)為傅里葉變換以后的結果,根據歐拉公式可得,每個F(u,v)值都為復數,由實部和虛部組成
代碼示例
1 void dft(short** in_array, double** re_array, double** im_array, long height, long width) 2 { 3 double re, im, temp; 4 5 for (int i = 0; i < height; i++){ 6 for (int j = 0; j < width; j++){ 7 re = 0; 8 im = 0; 9 10 for (int x = 0; x < height; x++){ 11 for (int y = 0; y < width; y++){ 12 temp = (double)i * x / (double)height + 13 (double)j * y / (double)width; 14 re += in_array[x][y] * cos(-2 * pi * temp); 15 im += in_array[x][y] * sin(-2 * pi * temp); 16 } 17 } 18 19 re_array[i][j] = re; 20 im_array[i][j] = im; 21 } 22 } 23 printf("dft done\n"); 24 }
傅里葉譜
相角
功率譜
傅里葉變換頻譜圖
對於上面得兩幅圖案,在區間[0, M-1]中,變換數據由兩個在點M/2處碰面的背靠背的半個周期組成
針對顯示和濾波的目的,在該區間中有一個完整的變換周期更加方便,因為完整周期中數據是連續的
我們希望得到上圖所示的圖案
傅里葉變換的平移性質
因此對每個f(x, y)項乘以(-1)x+y可達目的
代碼示例
1 void fre_spectrum(short **in_array, short **out_array, long height, long width) 2 { 3 double re, im, temp; 4 int move; 5 6 for (int i = 0; i < height; i++){ 7 for (int j = 0; j < width; j++){ 8 re = 0; 9 im = 0; 10 11 for (int x = 0; x < height; x++){ 12 for (int y = 0; y < width; y++){ 13 temp = (double)i * x / (double)height + 14 (double)j * y / (double)width; 15 move = (x + y) % 2 == 0 ? 1 : -1; 16 re += in_array[x][y] * cos(-2 * pi * temp) * move; 17 im += in_array[x][y] * sin(-2 * pi * temp) * move; 18 } 19 } 20 21 out_array[i][j] = (short)(sqrt(re*re + im*im) / sqrt(width*height)); 22 if (out_array[i][j] > 0xff) 23 out_array[i][j] = 0xff; 24 else if (out_array[i][j] < 0) 25 out_array[i][j] = 0; 26 27 } 28 } 29 }
執行結果
旋轉性質
即f(x, y)旋轉一個角度,F(u, v)旋轉相同的角度
二維離散傅里葉反變換
代碼示例
1 void idft(double** re_array, double** im_array, short** out_array, long height, long width) 2 { 3 double real, temp; 4 5 for (int i = 0; i < height; i++){ 6 for (int j = 0; j < width; j++){ 7 real = 0; 8 9 for (int x = 0; x < height; x++){ 10 for (int y = 0; y < width; y++){ 11 temp = (double)i * x / (double)height + 12 (double)j * y / (double)width; 13 14 real += re_array[x][y] * cos(2 * pi * temp) - 15 im_array[x][y] * sin(2 * pi * temp); 16 } 17 } 18 19 out_array[i][j] = (short)(real / sqrt(width*height)); 20 if (out_array[i][j] > 0xff) 21 out_array[i][j] = 0xff; 22 else if (out_array[i][j] < 0) 23 out_array[i][j] = 0; 24 } 25 } 26 printf("idft done\n"); 27 }
經驗證,圖像經傅里葉變換,然后再反變換以后可恢復原圖
改進
本篇文章只是按照二維離散傅里葉變換公式進行了實現,在測試的過程中發現,執行速度真的是非常慢,算法時間復雜度O(n4),等以后有時間再對這段代碼進行優化。