youtome 发表于 2024-1-25 11:22

C语言写的FFT代码

#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 = {1.0, 0.707, 0.0, -0.707};
float   x_r = {1, 1, 1, 1, 0, 0, 0, 0};//输入数据，此处设为8个
float   x_i;                         //N=8


/**
* 初始化输出虚部
*/
static void fft_init( void )
{
    int i;
    for(i=0; i<N; i++)   x_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 = xx_r[ bit_rev ];
}

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;
            tmp_imag = x_i;
            temp = x_r;
            
            /* 蝶形算法 */
            x_r   = tmp_real + ( tw1*x_r - tw2*x_i );
            x_i   = tmp_imag + ( tw2*x_r + tw1*x_i );
            /* X = A(k)+WB(k)
               * X = 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 = tmp_real + ( tw1*temp - tw2*x_i );
            x_i = tmp_imag + ( tw2*temp + tw1*x_i );*/
      x_r   = tmp_real - ( tw1* temp - tw2*x_i );
            x_i   = tmp_imag - ( tw2* temp + tw1*x_i );
            
            printf("k=%d, x_r=%f, x_i=%f\n", k+p, x_r, x_i);
            printf("k=%d, x_r=%f, x_i=%f\n", k+p+step, x_r, x_i);
         }
         /* 开跳!:) */
         k += 2*step;
       }   
    }
}

void display( void )
{
    printf("\n\n");
    int   i;
    for(i=0; i<N; i++)
       printf("%f\t%f\n", x_r, x_i);
}

int main( void )
{
    fft_init( );                //初始化
    bitrev( );                //将输入直接按FFT计算要求排序，如8点FFT计算，排序为x、x、x、x、x、x、x、x
    fft( );                        //进行FFT计算
    display( );                //显示计算结果
   
    system( "pause" );
    return 1;
}

rickluo 发表于 2024-1-26 17:09



点赞支持！

MT6701QT增量式磁编码控制器

feifeifeichang 发表于 2024-4-25 09:43

这个得点赞
页: [1]
查看完整版本: C语言写的FFT代码