fft程序：快速傅里叶变换(FFT)c语言实现

基本原理在这里就不多讲了，可以看看其他高浏览量的博文，这篇文章针对c语言的实现

复数运算算子

        我们都知道C语言本身是没有复数运算的，很多DSP、单片机要用到也没有开源库可以使用复数运算，针对FFT在硬件上运行只能手动从底层开始

定义复数类型

        这里用最简单高效的方法——结构体

struct complex{ double real; double image;};

复数加法

struct complex complex_add(struct complex c1,struct complex c2) //复数加法{ struct complex p; p.real = c1.real + c2.real; p.image = c1.image + c2.image; return p;}

复数减法

struct complex complex_sub(struct complex c1,struct complex c2) //复数减{ struct complex p; p.real = c1.real - c2.real; p.image = c1.image - c2.image; return p;}

复数乘法

struct complex complex_multi(struct complex c1,struct complex c2) //复数乘法{ struct complex c3; c3.real=c1.real*c2.real - c1.image *c2.image; c3.image = c2.real* c1.image + c1.real*c2.image; return c3;};

复数取模

double mold_length(struct complex c) //幅度{ return sqrt(c.real * c.real + c.image * c.image);};

旋转因子

    教科书上的旋转因子一般长成这样

 在其中抽取一个蝶形详细分析，注意这其中上半部分是相加【+】下半部分是相减【-】，而且下半部分都会带有一个旋转因子

 为了方便编程，把它表达成数学式子，这其中前一级与后一级的关系用等式表达

首先【j】与【k】的间隔可以容易观察出规律，每增加一级，间隔就会扩大两倍，换句话说【j】与【k】的间搁可以表达为下面的式子

 我们可以简单验证一下：

比如第【3】级 【k=2】的情况

 对照结果，和计算一样，程序成功

完整程序

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#define PI 3.1415926struct complex{ double real; double image;};struct complex complex_add(struct complex c1,struct complex c2);struct complex complex_sub(struct complex c1,struct complex c2);struct complex complex_multi(struct complex c1,struct complex c2);struct complex rotation_factor(int N,int n,int k);double mold_length(struct complex c);void fft(int len, struct complex in_x[],struct complex out_y[]);int main(){ int sam[8] = {1,-1,1,-1,1,-1,1,-1}; int jhg[8] = {0,0,0,0,0,0,0,0}; struct complex x[8]; struct complex y[8]; for(int i=0;i<8;i++) { x[i].real = sam[i]; x[i].image = jhg[i]; } printf("时域序列\n"); for(int i=0;i<8;i++) { printf("(%.2f, %.2fi) \n",x[i].real,x[i].image); } fft(8,x,y); printf("频域序列\n"); for(int i=0;i<8;i++) { printf("(%.2f, %.2fi)\n",y[i].real,y[i].image); } return 0;}struct complex complex_add(struct complex c1,struct complex c2) //复数加法{ struct complex p; p.real = c1.real + c2.real; p.image = c1.image + c2.image; return p;}struct complex complex_sub(struct complex c1,struct complex c2) //复数减{ struct complex p; p.real = c1.real - c2.real; p.image = c1.image - c2.image; return p;}struct complex complex_multi(struct complex c1,struct complex c2) //复数乘法{ struct complex c3; c3.real=c1.real*c2.real - c1.image *c2.image; c3.image = c2.real* c1.image + c1.real*c2.image; return c3;};struct complex rotation_factor(int N,int n,int k) //旋转因子{ struct complex w; w.real = cos(2* PI * n * k /N); w.image = - sin(2* PI * n * k /N); return w;}double mold_length(struct complex c) //幅度{ return sqrt(c.real * c.real + c.image * c.image);};int reverse_num(int l,int oringin_num) //反位序{ int q=0,m=0; for(int k=l-1;k>=0;k--) { q = oringin_num % 2; m += q*(1<<k); oringin_num = oringin_num / 2; } return m;}void fft(int len, struct complex in_x[],struct complex out_y[]){ /* param len 序列长度，目前只能是2的指数 param in_x输入的序列 param out_y输出的序列 */ int l,k,r,z,dist,q,j; //l是级 struct complex w,tmp; struct complex in_x_mem[len]; l = log2(len); for(k=0;k<len;k++) { in_x_mem[k] = in_x[reverse_num(l,k)]; //反位序号操作 } for(r = 0;r<l;r++) //遍历每一级蝶形运算 { dist = 1<<r; //提前计算每一级的间隔距离 for( j=0;j<dist;j++ ) //计算策略是拆成上下两组，先上计算，后下计算，j是计算的起始序号 { for(k=j;k<len;k+=(dist<<1)) //不好解释，得画图理解 { q = k+dist; //q同一组蝶形运算第二个序号 z = k << (l - r -1); //确定旋转因子的指数 w = rotation_factor(len,1,z); //由于不是并行计算，必须先缓存 tmp = in_x_mem[k]; in_x_mem[k] = complex_add( in_x_mem[k] , complex_multi(w, in_x_mem[q]) ); in_x_mem[q] = complex_sub(tmp , complex_multi(w, in_x_mem[q]) ); } } } memcpy(out_y,in_x_mem,len*sizeof(struct complex));}