C語言寫的FFT代碼

[復制鏈接]

電梯直達

樓主

ID:996634 發表于 2021-12-26 17:19 | 只看該作者回帖獎勵

|倒序瀏覽 |閱讀模式

/* 新手上路還望見諒。 *

#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;
}