C語(yǔ)言寫(xiě)的FFT代碼

ID:996634 · 發(fā)表于 2021-12-26 17:19

/* 新手上路還望見(jiàn)諒。 *

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 8 //64 輸入樣本總數(shù)
#define M 3 //DFT運(yùn)算層數(shù) //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}; //輸入數(shù)據(jù)，此處設(shè)為8個(gè)
float x_i[N]; //N=8
/**
* 初始化輸出虛部
*/
static void fft_init( void )
{
int i;
for(i=0; i<N; i++) x_i[i] = 0.0;
}
/**
* 反轉(zhuǎn)算法.將時(shí)域信號(hào)重新排序.
* 這個(gè)算法有改進(jìn)的空間
*/
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代表正要計(jì)算的當(dāng)前層，gr_num代表當(dāng)前層的顆粒數(shù)
float tmp_real, tmp_imag, temp; // 臨時(shí)變量, 記錄實(shí)部
float tw1, tw2;// 旋轉(zhuǎn)因子,tw1為旋轉(zhuǎn)因子的實(shí)部cos部分, tw2為旋轉(zhuǎn)因子的虛部sin部分.
int step; // 步進(jìn)
int sample_num; // 顆粒的樣本總數(shù)(各層不同, 因?yàn)楦鲗宇w粒的輸入不同)
/* 對(duì)層循環(huán) */
for(cur_layer=1; cur_layer<=M; cur_layer++)
{
/* 求當(dāng)前層擁有多少個(gè)顆粒(gr_num) */
gr_num = 1;
i = M - cur_layer;
while(i > 0)
{
i--;
gr_num *= 2;
}
/* 每個(gè)顆粒的輸入樣本數(shù)N' */
sample_num = (int)pow(2, cur_layer);
/* 步進(jìn). 步進(jìn)是N'/2 */
step = sample_num/2;
/* */
k = 0;
/* 對(duì)顆粒進(jìn)行循環(huán) */
for(i=0; i<gr_num; i++)
{
/*
* 對(duì)樣本點(diǎn)進(jìn)行循環(huán), 注意上限和步進(jìn)
*/
for(p=0; p<sample_num/2; p++)
{
// 旋轉(zhuǎn)因子, 需要優(yōu)化...
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) 的性質(zhì)可以優(yōu)化這里*/
/*旋轉(zhuǎn)因子, 需要優(yōu)化...
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ā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[i], x_i[i]);
}
int main( void )
{
fft_init( ); //初始化
bitrev( ); //將輸入直接按FFT計(jì)算要求排序，如8點(diǎn)FFT計(jì)算，排序?yàn)閤[0]、x[4]、x[2]、x[6]、x[1]、x[5]、x[3]、x[7]
fft( ); //進(jìn)行FFT計(jì)算
display( ); //顯示計(jì)算結(jié)果
system( "pause" );
return 1;
}

復(fù)制代碼

帳號(hào)		自動(dòng)登錄	找回密碼
密碼			立即注冊(cè)

C語(yǔ)言寫(xiě)的FFT代碼

評(píng)分