直接上代码了
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...