C語言數據結構算法之實現快速傅立葉變換
本實例將實現二維快速傅立葉變換，同時也將借此實例學習用c語言實現矩陣的基本操作、復數的基本掾作，復習所學過的動態內存分配、文件操作、結構指針的函數調用等內容。
很久以來，傅立葉變換一直是許多領域，如線性系統、光學、概率論、量子物理、天線、數字圖像處理和信號分析等的一個基本分析工具，但是即便使用計算速度驚人的計算機計算離散傅立葉變換所花費的時間也常常是難以接受的，因此導致了快速傅立葉變換（FFT)的產生。
本實例將對一個二維數組進行正、反快速傅立葉變換。正傅立葉變換時dfft()函數先調用fft()按行對數組進行變換，再對結果調用fft()按列進行變換，此時完成了快速傅立葉變換，再調用rdfft()函數進行傅立葉逆變換。如果程序設計正確的話，變換的結果應與原數組相同。
實例代碼：
#include <stdio.h>#include <stdlib.h>#include <math.h>#define PI 3.14159265358979323846 struct COMPLEX{ float re; float im;} cplx , * Hfield , * S , * R , * w; int n , m;int ln , lm; void initiate ();void dfft ();void rdfft ();void showresult (); void fft (int l , int k);int reverse (int t , int k);void W (int l);int loop (int l);void conjugate (); void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);struct COMPLEX * Hread(int i , int j);void Hwrite (int i , int j , struct COMPLEX x); void main (){ initiate (); printf("/n原始數據:/n"); showresult(); getchar (); dfft (); printf("/n快速復利葉變換后的結果:/n"); showresult (); getchar (); rdfft (); printf("/n快速復利葉逆變換后的結果:/n"); showresult (); getchar (); free (Hfield);} void initiate (){//程序初始化操作，包括分配內存、讀入要處理的數據、進行顯示等 FILE * df; df = fopen ("data.txt" , "r"); fscanf (df , "%5d" , &n); fscanf (df , "%5d" , &m); if ((ln = loop (n)) == -1) { printf (" 列數不是2的整數次冪 "); exit (1); } if ((lm = loop (m)) == -1) { printf (" 行數不是2的整數次冪 "); exit (1); } Hfield = (struct COMPLEX *) malloc (n * m * sizeof (cplx)); if (fread (Hfield , sizeof (cplx) , m * n , df) != (unsigned) (m * n)) { if (feof (df)) printf (" Premature end of file "); else printf (" File read error "); } fclose (df);} void dfft (){//進行二維快速復利葉變換 int i , j; int l , k; l = n; k = ln; w = (struct COMPLEX *) calloc (l , sizeof (cplx)); R = (struct COMPLEX *) calloc (l , sizeof (cplx)); S = (struct COMPLEX *) calloc (l , sizeof(cplx)); W (l); for ( i = 0 ; i < m ; i++ ) {//按行進行快速復利葉變換 for (j = 0 ; j < n ; j++) { S[j].re = Hread (i , j)->re; S[j].im = Hread (i , j)->im; } fft(l , k); for (j = 0 ; j < n ; j++) Hwrite (i , j , R[j]); } free (R); free (S); free (w); l = m; k = lm; w = (struct COMPLEX *) calloc (l , sizeof (cplx)); R = (struct COMPLEX *) calloc (l , sizeof (cplx)); S = (struct COMPLEX *) calloc (l , sizeof (cplx)); W (l); for (i = 0 ; i < n ; i++) {//按列進行快速復利葉變換 for(j = 0 ; j < m ; j++) { S[j].re = Hread(j , i)->re; S[j].im = Hread(j , i)->im; } fft(l , k); for (j = 0 ; j < m ; j++) Hwrite (j , i , R[j]); } free (R); free (S); free (w);} void rdfft (){ conjugate (); dfft (); conjugate ();} void showresult (){ int i , j; for (i = 0 ; i < m ; i++) { printf ( " /n第%d行/n " , i); for (j = 0 ; j < n ; j++) { if (j % 4 == 0) printf (" /n "); printf(" (%5.2f,%5.2fi) " , Hread (i , j)->re , Hread (i , j)->im); } }} void fft (int l , int k){ int i , j , s , nv , t; float c; struct COMPLEX mp , r; nv = l; c = (float) l; c = pow (c , 0.5); for (i = 0 ; i < k ; i++) { for (t = 0 ; t < l ; t += nv) { for (j = 0 ; j < nv / 2 ; j++) { s = (t + j) >> (k - i -1); s = reverse(s , k); r.re = S[t + j].re; r.im = S[t + j].im; mul (&w[s] , &S[t + j + nv / 2] , &mp);/////////講解傳遞結構指針和結構本身的區別 add (&r , &mp , &S[t + j]); sub (&r , &mp , &S[t + j + nv / 2]); } } nv = nv >> 1; } for (i = 0 ; i < l ; i++) { j = reverse(i , k); R[j].re = S[i].re / c; R[j].im = S[i].im / c; }} int reverse (int t , int k){ int i , x , y; y = 0; for (i = 0 ; i < k ; i++) { x = t & 1; t = t >> 1; y = (y << 1) + x; } return y;} void W (int l){ int i; float c , a; c = (float) l; c = 2 * PI / c; for (i = 0 ; i < l ; i++) { a = (float) i; w[i].re = (float) cos(a * c); w[i].im = -(float) sin(a * c); }} int loop (int l){//檢驗輸入數據是否為2的整數次冪，如果是返回用2進制表示時的位數 int i , m; if (l != 0) { for (i = 1 ; i < 32 ; i++) { m = l >> i; if (m == 0) break; } if (l == (1 << (i - 1))) return i - 1; } return -1;} void conjugate (){//求復數矩陣的共軛矩陣 int i , j; for (i = 0 ; i < m ; i++) { for (j = 0 ; j < n ; j++) { Hread (i , j)->im *= -1; } }} struct COMPLEX * Hread (int i , int j){//按讀矩陣方式返回Hfield中指定位置的指針 return (Hfield + i * n + j);} void Hwrite (int i , int j , struct COMPLEX x){//按寫矩陣方式將復數結構x寫到指定的Hfield位置上 (Hfield + i * n + j)->re = x.re; (Hfield + i * n + j)->im = x.im;} void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定義復數加法 z->re = x->re + y->re; z->im = x->im + y->im; } void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定義復數減法 z->re = x->re - y->re; z->im = x->im - y->im;} void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定義復數乘法 z->re = (x->re) * (y->re) - (x->im) * (y->im); z->im = (x->im) * (y->re) + (x->re) * (y->im);}感謝閱讀，希望能幫助到大家，謝謝大家對本站的支持！
新聞熱點
疑難解答
圖片精選
