/* 新手上路還望見諒。 *
- #include <stdio.h>
- #include <math.h>
- #include <stdlib.h>
- #define N 8 //64 輸入樣本總數
- #define M 3 //DFT運算層數 //2^m=N
- #define PI 3.1415926
- float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
- float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0}; //輸入數據，此處設為8個
- float x_i[N]; //N=8
- /**
- * 初始化輸出虛部
- */
- static void fft_init( void )
- {
- int i;
- for(i=0; i<N; i++) x_i[i] = 0.0;
- }
-
- /**
- * 反轉算法.將時域信號重新排序.
- * 這個算法有改進的空間
- */
- static void bitrev( void )
- {
- int p=1, q, i;
- int bit_rev[ N ];
- float xx_r[ N ];
-
- bit_rev[ 0 ] = 0;
- while( p < N )
- {
- for(q=0; q<p; q++)
- {
- bit_rev[ q ] = bit_rev[ q ] * 2;
- bit_rev[ q + p ] = bit_rev[ q ] + 1;
- }
- p *= 2;
- }
-
- for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];
-
- for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];
- }
- void fft( void )
- { fp = fopen("log2.txt", "a+");//此處
- int cur_layer, gr_num, i, k, p; //cur_layer代表正要計算的當前層，gr_num代表當前層的顆粒數
- float tmp_real, tmp_imag, temp; // 臨時變量, 記錄實部
- float tw1, tw2;// 旋轉因子,tw1為旋轉因子的實部cos部分, tw2為旋轉因子的虛部sin部分.
-
- int step; // 步進
- int sample_num; // 顆粒的樣本總數(各層不同, 因為各層顆粒的輸入不同)
-
- /* 對層循環 */
- for(cur_layer=1; cur_layer<=M; cur_layer++)
- {
- /* 求當前層擁有多少個顆粒(gr_num) */
- gr_num = 1;
- i = M - cur_layer;
- while(i > 0)
- {
- i--;
- gr_num *= 2;
- }
-
- /* 每個顆粒的輸入樣本數N' */
- sample_num = (int)pow(2, cur_layer);
- /* 步進. 步進是N'/2 */
- step = sample_num/2;
-
- /* */
- k = 0;
-
- /* 對顆粒進行循環 */
- for(i=0; i<gr_num; i++)
- {
- /*
- * 對樣本點進行循環, 注意上限和步進
- */
- for(p=0; p<sample_num/2; p++)
- {
- // 旋轉因子, 需要優化...
- tw1 = cos(2*PI*p/pow(2, cur_layer));
- tw2 = -sin(2*PI*p/pow(2, cur_layer));
-
- tmp_real = x_r[k+p];
- tmp_imag = x_i[k+p];
- temp = x_r[k+p+step];
-
- /* 蝶形算法 */
- x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
- x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
- /* X[k] = A(k)+WB(k)
- * X[k+N/2] = A(k)-WB(k) 的性質可以優化這里*/
- /*旋轉因子, 需要優化...
- tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
- tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
- x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
- x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );*/
- x_r[k+p+step] = tmp_real - ( tw1* temp - tw2*x_i[k+p+step] );
- x_i[k+p+step] = tmp_imag - ( tw2* temp + tw1*x_i[k+p+step] );
-
- printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
- printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);
- }
- /* 開跳!:) */
- k += 2*step;
- }
- }
- }
- void display( void )
- {
- printf("\n\n");
- int i;
- for(i=0; i<N; i++)
- printf("%f\t%f\n", x_r[i], x_i[i]);
- }
- int main( void )
- {
- fft_init( ); //初始化
- bitrev( ); //將輸入直接按FFT計算要求排序，如8點FFT計算，排序為x[0]、x[4]、x[2]、x[6]、x[1]、x[5]、x[3]、x[7]
- fft( ); //進行FFT計算
- display( ); //顯示計算結果
-
- system( "pause" );
- return 1;
- }
|
QQ好友和群
QQ空間
騰訊微博
騰訊朋友
收藏
淘帖
頂
踩
