C语言写的FFT代码
FFT代码
#include
#include
#include
#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 } /** * 反转算法.将时域信号重新排序. * 这个算法有改进的空间 */ 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 { bit_rev[ q ] = bit_rev[ q ] * 2; bit_rev[ q + p ] = bit_rev[ q ] + 1; } p *= 2; } for(i=0; i for(i=0; 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 { /* * 对样本点进行循环, 注意上限和步进 */ for(p=0; 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 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; }