FFT算法的原理详解

/*

功能：将input里的数据进行快速傅里叶变换 并且输出

*/

#include

#include

#define FFT_LENGTH 8

double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};

struct complex1{ //定义一个复数结构体

double real; //实部

double image; //虚部

};

//将input的实数结果存放为复数

struct complex1 result_dat[8];

/*

虚数的乘法

*/

struct complex1 con_complex(struct complex1 a,struct complex1 b){

struct complex1 temp;

temp.real=(a.real*b.real)-(a.image*b.image);

temp.image=(a.image*b.real)+(a.real*b.image);

return temp;

}

/*

简单的a的b次方

*/

int mypow(int a,int b){

int i,sum=a;

if(b==0)return 1;

for(i=1;i

sum*=a;

}

return sum;

}

/*

简单的求以2为底的正整数

*/

int log2(int n){

unsigned i=1;

int sum=1;

for(i;;i++){

sum*=2;

if(sum>=n)break;

}

return i;

}

/*

简单的交换数据的函数

*/

void swap(struct complex1 *a,struct complex1 *b){

struct complex1 temp;

temp=*a;

*a=*b;

*b=temp;

}

/*

dat为输入数据的数组

N为抽样次数 也代表周期 必须是2^N次方

*/

void fft(struct complex1 dat[],unsigned char N){

/*最终 dat_buf计算出 当前蝶形运算奇数项与W 乘积

dat_org存放上一个偶数项的值

*/

struct complex1 dat_buf,dat_org;

/* L为几级蝶形运算 也代表了2进制的位数

n为当前级蝶形的需要次数 n最初为N/2 每级蝶形运算后都要/2

i j为倒位时要用到的自增符号 同时 i也用到了L碟级数 j是计算当前碟级的计算次数

re_i i_copy均是倒位时用到的变量

k为当前碟级 cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k

*/

unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;

//经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N 次

unsigned char fft_counter=0;

//在此要进行补2 N必须是2^n 在此略

//蝶形级数 (L级)

L=log2(N);

//计算每级蝶形计算的次数(这里只是一个初始值) 之后每次要/2

//n=N/2;

//对dat的顺序进行倒位

for(i=1;i

i_copy=i;

re_i=0;

for(j=L-1;j>0;j--){

//判断i的副本最低位的数字 并且移动到最高位 次高位 ..

//re_i为交换的数 每次它的数字是不能移动的 并且循环之后要清0

re_i|=((i_copy&0x01)<

i_copy>>=1;

}

swap(&dat,&dat[re_i]);

}

//进行fft计算

for(i=0;i

fft_flag=1;

fft_counter=0;

for(j=0;j

if(fft_counter==mypow(2,i)){ //控制隔几次，运算几次

fft_flag=0;

}else if(fft_counter==0){ //休止结束，继续运算

fft_flag=1;

}

//当不判断这个语句的时候 fft_flag保持 这样就可以持续运算了

if(fft_flag){

dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));

dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));

dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);

//计算 当前蝶形运算奇数项与W 乘积

dat_org.real=dat[j].real;

dat_org.image=dat[j].image; //暂存

dat[j].real=dat_org.real+dat_buf.real;

dat[j].image=dat_org.image+dat_buf.image;

//实部加实部 虚部加虚部

dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;

dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;

//实部减实部 虚部减虚部

k++;

fft_counter++;

}else{

fft_counter--; //运算几次，就休止几次

k=0;

}

}

}

}

void main{

int i;

//先将输入信号转换成复数

for(i=0;i

result_dat.image=0;

//输入信号是二维的，暂时不存在复数

result_dat.real=input;

//result_dat.real=10;

//输入信号都为实数

}

fft(result_dat,FFT_LENGTH);

for(i=0;i

input=sqrt(result_dat.real*result_dat.real+result_dat.image*result_dat.image);

//取模

printf("%lfn",input);

}

while(1);

}

这个程序中input这个数组是输入信号，在这里只模拟抽样了8次，输出的数据也是input，如果想看其它序列的话，可以改变FFT_LENGTH 的值以及 input里的内容，程序输出的是实部和虚部的模