http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-21
一、徹底理解傅里葉變換
快速傅里葉變換(Fast Fourier Transform)是離散傅里葉變換的一種快速算法，簡稱FFT，通過FFT可以將一個信號從時域變換到頻域。
模擬信號經過A/D轉換變為數字信號的過程稱為采樣。為保證采樣后信號的頻譜形狀不失真，采樣頻率必須大於信號中最高頻率成分的2倍，這稱之為采樣定理。
| 假設采樣頻率為fs，采樣點數為N，那么FFT結果就是一個N點的復數，每一個點就對應着一個頻率點，某一點n(n從1開始)表示的頻率為：fn=(n-1)*fs/N。 舉例說明：用1kHz的采樣頻率采樣128點，則FFT結果的128個數據即對應的頻率點分別是0，1k/128，2k/128，3k/128，…，127k/128 Hz。 這個頻率點的幅值為：該點復數的模值除以N/2(n=1時是直流分量，其幅值是該點的模值除以N)。 |
二、傅里葉變換的C語言編程
1、碼位倒序。
假設一個N點的輸入序列，那么它的序號二進制數位數就是t=log2N.
碼位倒序要解決兩個問題：①將t位二進制數倒序;②將倒序后的兩個存儲單元進行交換。
如果輸入序列的自然順序號i用二進制數表示，例如若最大序號為15，即用4位就可表示n3n2n1n0，則其倒序后j對應的二進制數就是n0n1n2n3，那么怎樣才能實現倒序呢?利用C語言的移位功能!
/*變址計算，將x(n)碼位倒置*/ void Reverse(complex *IN_X, int n) { int row = log(n) / log(2); //row為FFT的N個輸入對應的最大列數 complex temp; //臨時交換變量 unsigned short i = 0, j = 0, k = 0; double t; for (i = 0; i < row; i++) //從第0個序號到第N-1個序號 { k = i; //當前第i個序號 j = 0; //存儲倒序后的序號，先初始化為0 t = row; //共移位t次 while ((t--) > 0) //利用按位與以及循環實現碼位顛倒 { j = j << 1; j |= (k & 1); //j左移一位然后加上k的最低位 k = k >> 1; //k右移一位，次低位變為最低位 /*j為倒序，k為正序， 從k右向左的值依次移位， j向左依次填入對應的k的右向位*/ } if (j > i) //如果倒序后大於原序數，就將兩個存儲單元進行交換(判斷j>i是為了防止重復交換 { temp = IN_X[i]; IN_X[i] = IN_X[j]; IN_X[j] = temp; } } } //復數加法的定義 complex add(complex a, complex b) { complex c; c.real = a.real + b.real; c.img = a.img + b.img; return c; } //復數乘法的定義 complex mul(complex a, complex b) { complex c; c.real = a.real*b.real - a.img*b.img; c.img = a.real*b.img + a.img*b.real; return c; } //復數減法的定義 complex sub(complex a, complex b) { complex c; c.real = a.real - b.real; c.img = a.img - b.img; return c; }
2、第二個要解決的問題就是蝶形運算
對N個輸入行，m表示第m列蝶形。
①第1級(第1列)每個蝶形的兩節點“距離”為1，第2級每個蝶形的兩節點“距離”為2，第3級每個蝶形的兩節點“距離”為4，第4級每個蝶形的兩節點“距離”為8。由此推得，
第m級蝶形運算，每個蝶形的兩節點“距離”L=2^(m-1)。
②對於16點的FFT，第1級有16組蝶形，每組有1個蝶形;第2級有4組蝶形，每組有2個蝶形;第3級有2組蝶形，每組有4個蝶形;第4級有1組蝶形，每組有8個蝶形。由此可推出，
對於N點的FFT，第m級有N/(2L)組蝶形，每組有L=2^(m-1)個蝶形。
③旋轉因子
的確定
以16點FFT為例，第m級第k個旋轉因子為
，其中k=0～2^(m-1)-1，即第m級共有2^(m-1)個旋轉因子，根據旋轉因子的可約性，
，所以第m級第k個旋轉因子為
，其中k=0～2^(m-1)-1。
生成一個的序列，代碼如下：
/*產生一個周期歐拉形式的Wn的值*/ void InitWn(complex *w,int n) //n為輸入的個數,w為權值Wn { int k; for (k = 0; k < n; k++) { w[k].real = cos(2 * PI / n * k); //用歐拉公式計算旋轉因子 w[k].img = -1 * sin(2 * PI / n * k); } }
3、算法實現
我們已經知道，N點FFT從左到右共有log2N級蝶形，每級有N/2L組，每組有L個。所以FFT的C語言編程只需用3層循環即可實現：最外層循環完成每一級的蝶形運算(整個FFT共log2N級)，中間層循環完成每一組的蝶形運算(每一級有N/2L組)，最內層循環完成單獨1個蝶形運算(每一組有L個)。
/*快速傅里葉變換*/ void fft(complex *IN_X, int n) { int row = log(n) / log(2); //row為FFT的N個輸入對應的最大列數 int i = 0, j = 0, k = 0, L = 0; complex top, bottom, product; Reverse(IN_X,n); //調用變址函數 for (i = 0; i < row ; i++) /*第i級蝶形運算 */ { L = 1 << i; //L等於2的i次方,表示每一級蝶形中組間的距離，即一組中蝶形個數 for (j = 0; j < n; j += 2 * L) /*j為組偏移地址，每L個蝶形是一組，每級有N/2L組*/ { for (k = 0; k < L; k++) /*k為一組中蝶形的偏移地址，一個蝶形運算 每個group內的蝶形運算*/ { product = mul(IN_X[j + k + L], Wn[n*k/2/L]); top = add(IN_X[j + k], product); bottom = sub(IN_X[j + k], product); IN_X[j + k] = top; IN_X[j + k + L] = bottom; } } } }
4.全部代碼
#include <stdio.h> #include <math.h> #include <stdlib.h> #include <malloc.h> #define N 1000 #define PI atan(1)* 4 /*定義復數類型*/ typedef struct { double real; double img; }complex; complex x[N], *Wn; /*輸入序列,復數系數Wn*/ int size_x; /*size_x為采樣的信號個數，必須為2的倍數*/ void fft(complex *IN_X, int n); /*對輸入IN_X，快速傅里葉變換*/ void InitWn(complex *w, int n); /*生成一個長度為n的Wn歐拉形式序列*/ void Reverse(complex *IN_X, int n); /*對輸入IN_X地址*/ complex add(complex, complex); /*復數加法*/ complex mul(complex, complex); /*復數乘法*/ complex sub(complex, complex); /*復數減法*/ void output(); /*輸出快速傅里葉變換的結果*/ int main() { int i; /*輸出結果*/ system("cls"); //調用系統命令，CLS的功能為清除屏幕輸出 printf("輸出DIT方法實現的FFT結果\n"); printf("Please input the size of x:\n");//輸入序列的大小 scanf_s("%d", &size_x); printf("Please input the data in x[N]:\n");//輸入序列的實部和虛部 for (i = 0; i < size_x; i++) { printf("請輸入第%d個序列：", i); scanf_s("%lf%lf", &x[i].real, &x[i].img); } printf("輸出倒序后的序列\n"); Wn = (complex *)malloc(sizeof(complex) * size_x); //分配變換后的值的內存空間 InitWn(Wn , size_x);//調用變換核 fft(x, size_x);//調用快速傅里葉變換 printf("輸出FFT后的結果\n"); output();//調用輸出傅里葉變換結果函數 scanf_s("%d", &size_x); //free(Wn); return 0; } /*快速傅里葉變換*/ void fft(complex *IN_X, int n) { int row = log(n) / log(2); //row為FFT的N個輸入對應的最大列數 int i = 0, j = 0, k = 0, L = 0; complex top, bottom, product; Reverse(IN_X,n); //調用變址函數 for (i = 0; i < row ; i++) /*第i級蝶形運算 */ { L = 1 << i; //L等於2的i次方,表示每一級蝶形中組間的距離，即一組中蝶形個數 for (j = 0; j < n; j += 2 * L) /*j為組偏移地址，每L個蝶形是一組，每級有N/2L組*/ { for (k = 0; k < L; k++) /*k為一組中蝶形的偏移地址，一個蝶形運算 每個group內的蝶形運算*/ { product = mul(IN_X[j + k + L], Wn[n*k/2/L]); top = add(IN_X[j + k], product); bottom = sub(IN_X[j + k], product); IN_X[j + k] = top; IN_X[j + k + L] = bottom; } } } } /*產生一個周期歐拉形式的Wn的值*/ void InitWn(complex *w,int n) //n為輸入的個數,w為權值Wn { int k; for (k = 0; k < n; k++) { w[k].real = cos(2 * PI / n * k); //用歐拉公式計算旋轉因子 w[k].img = -1 * sin(2 * PI / n * k); } } /*變址計算，將x(n)碼位倒置*/ /* 碼位倒序要解決兩個問題： ①將t位二進制數倒序;②將倒序后的兩個存儲單元進行交換。 如果輸入序列的自然順序號i用二進制數表示， 例如若最大序號為15，即用4位就可表示n3n2n1n0， 則其倒序后j對應的二進制數就是n0n1n2n3 */ void Reverse(complex *IN_X, int n) { int row = log(n) / log(2); //row為FFT的N個輸入對應的最大列數 complex temp; //臨時交換變量 unsigned short i = 0, j = 0, k = 0; double t; for (i = 0; i < row; i++) //從第0個序號到第N-1個序號 { k = i; //當前第i個序號 j = 0; //存儲倒序后的序號，先初始化為0 t = row; //共移位t次 while ((t--) > 0) //利用按位與以及循環實現碼位顛倒 { j = j << 1; j |= (k & 1); //j左移一位然后加上k的最低位 k = k >> 1; //k右移一位，次低位變為最低位 /*j為倒序，k為正序， 從k右向左的值依次移位， j向左依次填入對應的k的右向位*/ } if (j > i) //如果倒序后大於原序數，就將兩個存儲單元進行交換(判斷j>i是為了防止重復交換 { temp = IN_X[i]; IN_X[i] = IN_X[j]; IN_X[j] = temp; } } } /*輸出傅里葉變換的結果*/ void output() { int i; printf("The result are as follows：\n"); for (i = 0; i < size_x; i++) { printf("%.4f", x[i].real); if (x[i].img >= 0.0001)printf("+%.4fj\n", x[i].img); else if (fabs(x[i].img) < 0.0001)printf("\n"); else printf("%.4fj\n", x[i].img); } } //復數加法的定義 complex add(complex a, complex b) { complex c; c.real = a.real + b.real; c.img = a.img + b.img; return c; } //復數乘法的定義 complex mul(complex a, complex b) { complex c; c.real = a.real*b.real - a.img*b.img; c.img = a.real*b.img + a.img*b.real; return c; } //復數減法的定義 complex sub(complex a, complex b) { complex c; c.real = a.real - b.real; c.img = a.img - b.img; return c; }