EEPW论坛

浅谈FFT以及FFT算法的基本实现（下）

       对于FFT算法C语言的实现，网上的方法层出不穷，介于本人比较懒（懒得看别人的程序），再加上自给自足丰衣足食的原则，我自己也写了一个个人认为比较通俗易懂的程序。并且为了帮助读者理解，我特意尽量减少了库函数的使用，一些基本的函数都是自己写的（难免有很多BUG），但是作为FFT算法已经够用了，目前这个程序只能处理2^N的数据，理论上来讲如果不够2^N的话可以在后面数列补0这种操作，当然为了简约文中也就没有实现，但是主要的功能还是具备的，下面将代码开源如下：

/*
功能：将input里的数据进行快速傅里叶变换  并且输出
*/
 
#include <stdio.h>
#include <math.h>
#define PI 3.1415926
#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<b;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<N/2;i++){
i_copy=i;
re_i=0;
for(j=L-1;j>0;j--){
//判断i的副本最低位的数字  并且移动到最高位  次高位  ..
//re_i为交换的数   每次它的数字是不能移动的 并且循环之后要清0
re_i|=((i_copy&0x01)<<j);  
i_copy>>=1;
}
swap(&dat[i],&dat[re_i]);
}
//进行fft计算
for(i=0;i<L;i++){
fft_flag=1;
fft_counter=0;
for(j=0;j<N;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<FFT_LENGTH;i++){
result_dat[i].image=0; 
//输入信号是二维的，暂时不存在复数
result_dat[i].real=input[i];
//result_dat[i].real=10;
//输入信号都为实数
}
fft(result_dat,FFT_LENGTH);
for(i=0;i<FFT_LENGTH;i++){
input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
//取模
printf("%lf\n",input[i]);
}
while(1);
}

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

本次浅谈FFT以及FFT算法的基本实现的全部内容就到这里，欢迎大家多多交流~坐等板砖