最近做一個東西，要用到快速傅里葉變換，抱着蛋疼的心態，自己嘗試寫了一下，遇到一些問題。
首先看一下什么叫做快速傅里葉變換（FFT）（來自Wiki）：
快速傅里葉變換（英語：Fast Fourier Transform, FFT），是離散傅里葉變換的快速算法，也可用於計算離散傅里葉變換的逆變換。快速傅里葉變換有廣泛的應用，如數字信號處理、計算大整數乘法、求解偏微分方程等等。
對於復數串行
，離散傅里葉變換公式為：
直接變換的計算復雜度是O(n^2).快速傅里葉變換可以計算出與直接計算相同的結果，但只需要O(nlogn)的計算復雜度。
參考：快速傅里葉變換Wiki
首先，看一個概念叫做N次單位復根：復數WN^n = 1。N次單位復根一共有N個，分別為：W(k,N) = cis(2*PI*k/N),k=0,1,...,n-1
其中cis(u) = e^i*u = cos(u) + isin(u)（歐拉公式）
參考：歐拉公式Wiki
N次單位復根有如下性質：
- W(k+N,N) = W(k,N)
- W(k+N/2,N) = -W(k,N)
- W(dk,dN) = W(k,N) (d>0)
庫利-圖基算法：
將長度為N的序列分成兩個N/2的子序列N1和N2，其中N1是原來序列的奇數項，N2為偶數項。
將離散傅里葉變換寫成：
N(x) = ∑(j=0,n-1)aj*xj x = W(k,N)
分解后有：
N1(x) = a0 + a2*x + a4*x^2 + ... + a(n-2)*x^(n/2-1)……1
N2(x) = a1 + a3*x + a5*x^2 + ... + a(n-1)*x^(n/2-1)……2
N(x) = N1(x^2) + x*N2(x^2)……3
就這樣將(W(0,N))^2,(W(1,N))^2,...,(W(N-1,N))^2一個一個往1，2式代，然后根據3式組合就可以得到結果了。
按照相同的辦法，將此式遞歸下去，最終就可以得到結果。這就是庫利-圖基算法的大概思想。
語言描述起來很抽象，使用蝴蝶網絡來表示，就可以一目了然：
如圖，展示了一個N=8的快速傅里葉變換執行流程。
按照這種算法，可以很快寫出遞歸算法類C偽代碼：
Array recursive_fft(a) { int n = a.lengh(); if(n=1) return a; W wn = e^2*PI*i/n; W w = 1; Array a0 = {a0,a2,a3……a(n-2)}; Array a1 = {a1,a3,a5……a(n-1)}; Array y0 = recursive(a[0]); Array y1 = recursive(a[1]); Array y; for(int k = 0;k<=n/2-1;++k) { y[k] = y0[k]+wy1[k]; y[k+n/2] = y0[k]-wy1[0]; w *= wn } return y; }
然而，FFT在實際應用中往往需求很高的效率，雖然都是O(n*logn)d的時間復雜度，但是在遞歸的實現算法中，隱含的常數更大（這點可以參考《算法導論》的相關內容，這里不做具體討論）。並且，遞歸的方案也不好用GPU實現。
所以需要實現迭代的FFT，參考算法導論相代碼實現如下：
#pragma once #include "RBVector2" class RBFFTTools { public: static void iterative_fft(RBVector2 *v,RBVector2 *out,int len) { bit_reverse_copy(v,out,len); int n = len; for (int s = 1; s <= log(n)/log(2);s++) { int m = pow(2, s); RBVector2 wm(cos(2 * PI / m), sin(2 * PI / m)); for (int k = 0; k < n ;k+=m) { RBVector2 w = RBVector2(1,0); for (int j = 0; j < m / 2;++j) { RBVector2 vv = out[k + j + m / 2]; RBVector2 t = RBVector2(w.x*vv.x - w.y*vv.y,w.x*vv.y+w.y*vv.x); RBVector2 u = out[k + j]; out[k + j] = RBVector2(u.x+t.x,u.y+t.y); out[k + j + m / 2] = RBVector2(u.x-t.x,u.y-t.y); w = RBVector2(w.x*wm.x - w.y*wm.y, w.x*wm.y + w.y*wm.x); } } } } //逆變換 static void iterative_ifft(RBVector2 *v, RBVector2 *out, int len) { bit_reverse_copy(v, out, len); int n = len; for (int s = 1; s <= log(n)/log(2); s++) { int m = pow(2, s); RBVector2 wm(cos(2 * PI / m), sin(2 * PI / m)); for (int k = 0; k < n; k += m) { RBVector2 w = RBVector2(1, 0); for (int j = 0; j < m / 2; ++j) { RBVector2 vv = out[k + j + m / 2]; RBVector2 t = RBVector2(w.x*vv.x - w.y*vv.y, w.x*vv.y + w.y*vv.x); RBVector2 u = out[k + j]; out[k + j] = RBVector2(u.x + t.x, u.y + t.y); out[k + j + m / 2] = RBVector2(u.x - t.x, u.y - t.y); w = RBVector2(w.x*wm.x + w.y*wm.y, -w.x*wm.y + w.y*wm.x); } } } for (int i = 0; i < n;++i) { out[i].x /= n; out[i].y /= n; } } private: static void bit_reverse_copy(RBVector2 src[], RBVector2 des[], int len) { for (int i = 0; i < len;i++) { des[bit_rev(i, len-1)] = src[i]; } } public: static unsigned int bit_rev(unsigned int v, unsigned int maxv) { unsigned int t = log(maxv + 1)/log(2); unsigned int ret = 0; unsigned int s = 0x80000000>>(31); for (unsigned int i = 0; i < t; ++i) { unsigned int r = v&(s << i); ret |= (r << (t-i-1)) >> (i); } return ret; } };
上面的代碼出現了兩個函數：bit_reverse_copy()和bit_rev()，這兩個函數用於將初始數據按照bit反轉順序進行排序。
bit反轉看下面例子就可了然：
| 000 | 000 |
| 001 | 100 |
| 010 | 010 |
| 011 | 110 |
| 100 | 001 |
| 101 | 101 |
| 110 | 011 |
| 111 | 111 |
即把所有bit位沖頭到尾顛倒一次。
上面iterative_ifft()是逆變換，逆變換很簡單就可以修改出來，只需要根據W的性質進行修改W的更新，以及除以一個N就可以了。
運行結果如下：
參考：算法導論 第30章
基於GPU的FFT
這里使用OpenGL的Fragment Shader來實現FFT，我使用了兩個浮點紋理作為輸入，一個紋理表示蝴蝶網絡，一個表示當前的輸入數據：
一個N（2的冪）的FFT需要迭代log_2(N)次，蝴蝶網絡和輸入數據每次迭代都會更新，更新的數據用於下次迭代輸入，輸出的數據使用FBO來暫存，效率還可以。
蝴蝶網絡的單元結構如下：
輸入數據僅用了四個通道中的前兩個，用來存放數據的實部和虛部。
GPU算法可以通過蝴蝶網絡圖很容易的得出：
觀察此圖，可以發現：
每次迭代都需要輸出數據，輸出的數據由下圖計算方案決定：
W(r,N)就儲存在對應的蝴蝶網絡紋理中，取出來就可以了。
另外，每次迭代后都要對於蝴蝶網絡更新，更新的內容包括：
更新W和更新對應項的索引x，y值，由上圖可以知道，W只需下標乘以2得到Nnew，上標依次從零到Nnew，以此2循環即可。因為轉換一下就可以知道，W的變換其實是下面這個規律：
W(0,2) W(0,4) W(0,8) W(0,16) ……
W(1,2) W(1,4) W(1,8) W(1,16) ……
W(2,4) W(2,8) W(2,16) ……
W(3,4) W(3,8) W(3,16) ……
W(4,8) W(4,16) ……
W(5,8) .................
W(6,8) .................
W(7,8) .................
.................
而索引的變換就更加簡單了，加上當前的迭代次數即可。
思路很簡單，但是調試的時候很蛋疼，給出的代碼使用了OpenFrameworks提供的OpenGL框架：
#include "ofApp.h" #include "time.h" time_t t_s1,t_e1; time_t t1,t2; void ofApp::setup() { ofDisableArbTex(); imag_test.loadImage("1.jpg"); first = true; N = 1024; cur = 0; source = new RBVector2[N*N]; rev_source = new RBVector2[N*N]; out = new RBVector2[N*N]; //禁用Alpha混合，否則繪制到FBO會混合Alpha,造成數據丟失 ofDisableAlphaBlending(); for(int i=0;i<2;i++) { weight_index[i].allocate(N,N,GL_RGBA32F); res_back[i].allocate(N,N,GL_RGBA32F); data[i].allocate(N,N,ofImageType::OF_IMAGE_COLOR_ALPHA); } //初始化蝴蝶網絡 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { if(j%2==0) data[0].setColor(j,i,ofFloatColor(0.f,2.f,(float)((i*N+j+1)%N),(float)((i*N+j+1)/N))); else data[0].setColor(j,i,ofFloatColor(1.f,2.f,(float)((i*N+j-1)%N),(float)((i*N+j-1)/N))); } } //初始化原始數據，這里使用隨機數 for(int i=0;i<N*N;i++) { float s1 = ofRandomf(); float s2 = ofRandomf(); source[i] = RBVector2(s1,s2); } t_s1 = time(0); bit_reverse_copy(source,rev_source,N*N); t_e1= time(0); t2 = t_e1 - t_s1; //初始化到紋理數據 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { data[1].setColor(j,i,ofFloatColor(rev_source[i*N+j].x,rev_source[i*N+j].y,0,1)); } } data[0].update(); data[1].update(); weight_index[0].begin(); data[0].draw(0,0,N,N); weight_index[0].end(); res_back[0].begin(); data[1].draw(0,0,N,N); res_back[0].end(); ofSetWindowShape(1280,720); weight_index_shader.load("normalv.c","generate_weight_and_index_f.c"); gpu_fft_shader.load("normalv.c","gpu_fft_f.c"); t_s1 = time(0); //CPU變換 RBFFTTools::iterative_fft(source,out,N*N); t_e1 = time(0); t1 = t_e1 - t_s1; } //-------------------------------------------------------------- void ofApp::draw() { if(first) { t_s1 = time(0); for(int i = 1;i<N*N;i*=2) { res_back[1-cur%2].begin(); //更新輸出 gpu_fft_shader.begin(); gpu_fft_shader.setUniform1i("size",N); gpu_fft_shader.setUniform1i("n_step",i); gpu_fft_shader.setUniform3f("iResolution",N,N,0); gpu_fft_shader.setUniformTexture("tex_index_weight",weight_index[cur].getTextureReference(),1); gpu_fft_shader.setUniformTexture("tex_res_back",res_back[cur].getTextureReference(),2); gpu_fft_shader.setUniformTexture("test",imag_test.getTextureReference(),4); ofRect(0,0,N,N); gpu_fft_shader.end(); res_back[1-cur%2].end(); //更新蝴蝶網絡數據 weight_index[1-cur%2].begin(); weight_index_shader.begin(); weight_index_shader.setUniform1i("size",N); weight_index_shader.setUniform1i("n_step",i); weight_index_shader.setUniform3f("iResolution",N,N,0); weight_index_shader.setUniformTexture("tex_input",weight_index[cur].getTextureReference(),1); ofRect(0,0,N,N); weight_index_shader.end(); weight_index[1-cur%2].end(); cur = 1 - cur%2; } t_e1= time(0); t2 += (t_e1 - t_s1); cout<<"==================↓gpu================"<<endl; //輸出GPU計算結果 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { //cout<<s.getColor(j,i).r<<","<<s.getColor(j,i).g<<endl; } } cout<<"==================↓cpu================"<<endl; //輸出CPU計算結果 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { int index = i*N+j; //cout<<out[index].x<<","<<out[index].y<<endl; } } cout<<"===================================="<<endl; //對比CPU和GPU計算結果 for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { float rx = s.getColor(j,i).r; float ry = s.getColor(j,i).g; int index = i*N+j; float delta = 0.1; if(abs(rx-out[index].x)>delta||abs(ry-out[index].y)>delta) { //數據在誤差范圍外，出錯 goto wrong; //cout<<s.getColor(j,i).r<<","<<s.getColor(j,i).g<<"|"<<out[index].x<<","<<out[index].y<<endl; } else { //cout<<s.getColor(j,i).r<<","<<s.getColor(j,i).g<<","<<s.getColor(j,i).b<<","<<s.getColor(j,i).a<<endl; } } } cout<<"CPU:"<<t1<<endl<<"GPU:"<<t2<<endl; first = !first; } if(false) { wrong: printf("wrong!!!\n"); } //繪制數據輸入紋理和數據輸出紋理 res_back[cur].draw(N,0,N,N); data[1].draw(0,0,N,N); }
更新數據Shader：
#version 130 uniform sampler2D tex_index_weight; uniform sampler2D tex_res_back; uniform sampler2D test; uniform int size; uniform int n_step; uniform vec3 iResolution; out vec4 outColor; void main() { vec2 tex_coord = gl_FragCoord.xy/iResolution.xy; float cur_x = gl_FragCoord.x - 0.5; float cur_y = gl_FragCoord.y - 0.5; vec2 outv; vec4 temp = texture(tex_index_weight,tex_coord); vec2 weight = vec2(cos(temp.r/temp.g*2*3.141592653),sin(temp.r/temp.g*2*3.141592653)); vec2 _param2_index = temp.ba; vec2 temp_param2_index = _param2_index; temp_param2_index.x += 0.5; temp_param2_index.y += 0.5; vec2 param2_index = temp_param2_index/iResolution.xy; vec2 param1 = texture(tex_res_back,tex_coord).rg; vec2 param2 = texture(tex_res_back,param2_index).rg; int tex_coord_n1 = int((cur_y*size)+cur_x); int tex_coord_n2 = int((_param2_index.y*size)+_param2_index.x); if(tex_coord_n1<tex_coord_n2) { outv.r = param1.r + param2.r*weight.r-weight.g*param2.g; outv.g = param1.g +weight.r*param2.g + weight.g*param2.r; } else { outv.r = param2.r + param1.r*weight.r-weight.g*param1.g; outv.g = param2.g +weight.r*param1.g + weight.g*param1.r; } outColor = vec4(outv,0,1); }
更新蝴蝶網絡Shader：
#version 130 uniform sampler2D tex_input; uniform int size; uniform int n_total; uniform int n_step; in vec3 normal; uniform vec3 iResolution; out vec4 outColor; void main() { vec2 tex_coord = gl_FragCoord.xy/iResolution.xy; float cur_x = gl_FragCoord.x - 0.5; float cur_y = gl_FragCoord.y - 0.5; vec4 outv; int tex_coord_n = int((cur_y*size)+cur_x); //updata weight vec2 pre_w = texture(tex_input,tex_coord).rg; float i = pre_w.r; float n = pre_w.g; float new_i; float new_n; new_i = i; new_n = n*2; if(tex_coord_n%(n_step*4) > n_step*2-1) { new_i += n_step*2; } outv.r = new_i; outv.g = new_n; //updata index vec2 pre_index = texture(tex_input,tex_coord).ba; int x = int(pre_index.x); int y = int(pre_index.y); int ni = n_step*2; int new_tex_coord_n = tex_coord_n; if((tex_coord_n/ni)%2==0) { new_tex_coord_n += ni; } else { new_tex_coord_n -= ni; } outv.b = new_tex_coord_n%size; outv.a = new_tex_coord_n/size; outColor = outv; }
用time大概對比了一下GPU和CPU的計算時間，結果如下：
N = 512*512的情況下，誤差在0.1以內，在我的應用情況下是可以接受的，更高的數據意味着更多的迭代次數，會進一步降低數據精度，目前還沒有解決誤差問題。
運行50次，大概對比如下，不是很准確，僅供參考：
