標(biāo)題: C語言寫的FFT代碼 [打印本頁]
作者: 略知一二    時間: 2021-12-26 17:19
標(biāo)題: C語言寫的FFT代碼
/*   新手上路還望見諒。  *
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>

  5. #define   N     8     //64  輸入樣本總數(shù)
  6. #define    M     3   //DFT運(yùn)算層數(shù)     //2^m=N  
  7. #define    PI    3.1415926

  9. float   twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
  10. float   x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};  //輸入數(shù)據(jù)，此處設(shè)為8個
  11. float   x_i[N];                         //N=8


  14. /**
  15. * 初始化輸出虛部
  16. */
  17. static void fft_init( void )
  18. {
  19.     int i;
  20.     for(i=0; i<N; i++)   x_i[i] = 0.0;
  21. }

  23. /**
  24. * 反轉(zhuǎn)算法.將時域信號重新排序.
  25. * 這個算法有改進(jìn)的空間
  26. */
  27. static void bitrev( void )
  28. {
  29.     int    p=1, q, i;
  30.     int    bit_rev[ N ];  
  31.     float   xx_r[ N ];   
  32.    
  33.     bit_rev[ 0 ] = 0;
  34.     while( p < N )
  35.     {
  36.        for(q=0; q<p; q++)  
  37.        {
  38.            bit_rev[ q ]     = bit_rev[ q ] * 2;
  39.            bit_rev[ q + p ] = bit_rev[ q ] + 1;
  40.        }
  41.        p *= 2;
  42.     }
  43.    
  44.     for(i=0; i<N; i++)   xx_r[ i ] = x_r[ i ];   
  45.    
  46.     for(i=0; i<N; i++)   x_r[i] = xx_r[ bit_rev[i] ];
  47. }

  49. void fft( void )
  50. {   fp = fopen("log2.txt", "a+");//此處
  51.     int     cur_layer, gr_num, i, k, p;        //cur_layer代表正要計算的當(dāng)前層，gr_num代表當(dāng)前層的顆粒數(shù)
  52.     float   tmp_real, tmp_imag, temp;   // 臨時變量, 記錄實(shí)部
  53.     float   tw1, tw2;// 旋轉(zhuǎn)因子,tw1為旋轉(zhuǎn)因子的實(shí)部cos部分, tw2為旋轉(zhuǎn)因子的虛部sin部分.
  54.       
  55.     int    step;      // 步進(jìn)
  56.     int    sample_num;   // 顆粒的樣本總數(shù)(各層不同, 因?yàn)楦鲗宇w粒的輸入不同)
  57.    
  58.     /* 對層循環(huán) */
  59.     for(cur_layer=1; cur_layer<=M; cur_layer++)
  60.     {      
  61.        /* 求當(dāng)前層擁有多少個顆粒(gr_num) */
  62.        gr_num = 1;
  63.        i = M - cur_layer;
  64.        while(i > 0)
  65.        {
  66.            i--;
  67.            gr_num *= 2;
  68.        }
  69.       
  70.        /* 每個顆粒的輸入樣本數(shù)N' */
  71.        sample_num    = (int)pow(2, cur_layer);
  72.        /* 步進(jìn). 步進(jìn)是N'/2 */
  73.        step       = sample_num/2;
  74.       
  75.        /*  */
  76.        k = 0;
  77.       
  78.        /* 對顆粒進(jìn)行循環(huán) */
  79.        for(i=0; i<gr_num; i++)
  80.        {
  81.            /*
  82.             * 對樣本點(diǎn)進(jìn)行循環(huán), 注意上限和步進(jìn)
  83.             */
  84.            for(p=0; p<sample_num/2; p++)
  85.            {   
  86.               // 旋轉(zhuǎn)因子, 需要優(yōu)化...   
  87.               tw1 = cos(2*PI*p/pow(2, cur_layer));
  88.               tw2 = -sin(2*PI*p/pow(2, cur_layer));
  89.               
  90.               tmp_real = x_r[k+p];
  91.               tmp_imag = x_i[k+p];
  92.               temp = x_r[k+p+step];
  93.               
  94.               /* 蝶形算法 */
  95.               x_r[k+p]   = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
  96.               x_i[k+p]   = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
  97.               /* X[k] = A(k)+WB(k)
  98.                * X[k+N/2] = A(k)-WB(k) 的性質(zhì)可以優(yōu)化這里*/
  99.               /*旋轉(zhuǎn)因子, 需要優(yōu)化...
  100.               tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
  101.               tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
  102.               x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
  103.               x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );*/
  104.         x_r[k+p+step]   = tmp_real - ( tw1* temp - tw2*x_i[k+p+step] );
  105.               x_i[k+p+step]   = tmp_imag - ( tw2* temp + tw1*x_i[k+p+step] );
  106.               
  107.               printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
  108.               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]);
  109.            }
  110.            /* 開跳!:) */
  111.            k += 2*step;
  112.        }   
  113.     }
  114. }

  116. void display( void )
  117. {
  118.     printf("\n\n");
  119.     int   i;  
  120.     for(i=0; i<N; i++)
  121.        printf("%f\t%f\n", x_r[i], x_i[i]);
  122. }

  124. int main( void )
  125. {
  126.     fft_init( );                //初始化
  127.     bitrev( );                //將輸入直接按FFT計算要求排序，如8點(diǎn)FFT計算，排序?yàn)閤[0]、x[4]、x[2]、x[6]、x[1]、x[5]、x[3]、x[7]
  128.     fft( );                        //進(jìn)行FFT計算
  129.     display( );                //顯示計算結(jié)果
  130.    
  131.     system( "pause" );
  132.     return 1;
  133. }

復(fù)制代碼




歡迎光臨 (http://www.torrancerestoration.com/bbs/) Powered by Discuz! X3.1