DSP

网上找的纯C实现的FFT,与matlab计算结果完全一样

2019-07-13 18:44发布

直接上代码了 fft.c #include "math.h" #include "fft.h" //精度0.0001弧度 void conjugate_complex(int n,complex in[],complex out[]) { int i = 0; for(i=0;ireal = a.real + b.real; c->imag = a.imag + b.imag; } void c_sub(complex a,complex b,complex *c) { c->real = a.real - b.real; c->imag = a.imag - b.imag; } void c_mul(complex a,complex b,complex *c) { c->real = a.real * b.real - a.imag * b.imag; c->imag = a.real * b.imag + a.imag * b.real; } void c_div(complex a,complex b,complex *c) { c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag); c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag); } #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void Wn_i(int n,int i,complex *Wn,char flag) { Wn->real = cos(2*PI*i/n); if(flag == 1) Wn->imag = -sin(2*PI*i/n); else if(flag == 0) Wn->imag = -sin(2*PI*i/n); } //傅里叶变化 void fft(int N,complex f[]) { complex t,wn;//中间变量 int i,j,k,m,n,l,r,M; int la,lb,lc; /*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i



fft.h #ifndef __FFT_H__ #define __FFT_H__ typedef struct complex //复数类型 { float real; //实部 float imag; //虚部 }complex; #define PI 3.1415926535897932384626433832795028841971 /////////////////////////////////////////// void conjugate_complex(int n,complex in[],complex out[]); void c_plus(complex a,complex b,complex *c);//复数加 void c_mul(complex a,complex b,complex *c) ;//复数乘 void c_sub(complex a,complex b,complex *c); //复数减法 void c_div(complex a,complex b,complex *c); //复数除法 void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中 void ifft(int N,complex f[]); // 傅里叶逆变换 void c_abs(complex f[],float out[],int n);//复数数组取模 //////////////////////////////////////////// #endif

使用 fft(FFT_NPT, fft_buff); //进行FFT处理

点数必须为8,16,32,64,128,256...