基本原理在这里就不多讲了，可以看看其他高浏览量的博文，这篇文章针对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);
};旋转因子
教科书上的旋转因子一般长成这样把它变成可以运算的形式
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); //这里的PI是宏定义
return w;
}反位序操作
我用的方法是基于时间抽取的基2FFT （DIT-FFT）所以开始进入FFT之前要先把序列顺序取反位序
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;
}解释一下这个代码含义，用取余操作取出最低位，把最低位放在m中最高位置，取出后除2操作把每一位右移，直到为0
FFT
我们首先拿出祖传的蝶形运算信号流图，
在其中抽取一个蝶形详细分析，注意这其中上半部分是相加【+】下半部分是相减【-】，而且下半部分都会带有一个旋转因子，先记住这个特征，接下来编程会利用这个特征
为了方便编程，把它表达成数学式子，这其中前一级与后一级的关系用等式表达
现在问题就是k和j的确定，重新观察完整的蝶形图在第【1】级两两运算的节点序号是【0、1】【2、3】【4、5】【6、7】；在第【2】级两两运算的节点序号是【0、2】【1、3】【4、6】【5、7】；第【3】级运算的节点是【0、4】【1、5】【2、6】【3、7】。注意这里不要用图上的反位置序，而是反位后的顺序，也就是下图红色序号
首先【j】与【k】的间隔可以容易观察出规律，每增加一级，间隔就会扩大两倍，换句话说【j】与【k】的间搁可以表达为下面的式子
好了，这里我们破解了第一条规律，急需确定的是这个旋转因子当中的，这里我没有观察到规律，只好直接引用教材上的解释
我们可以简单验证一下：
比如第【3】级 【k=2】的情况
第【2】级【k=5】的情况
两次验证都符合实际FFT蝶形图的结果。
接下来看似没有急需破解的规律了，但是真正编程时候还是存在一个问题，【k】的取值不是从上到下按顺序的例如第一级的【k】取值为【0，2，4，6】第二级的取值为【0，1，4，5】第三级的取值为【0，1，2，3】。
这样如果按找从头到尾的顺序确实找不到什么规律，但是整体没规律，部分总会有规律吧，注意看我笔迹画的线为啥颜色是这样分的呢，没有发现很有规律吗？我们按找相同颜色分组，第一级的第一组两个蝶形运算的【k】之间相差距离刚好是【dist】的两倍,第二级的第一组两个蝶形运算的【k】之间相差的距离也是【dist】的两倍，但是到了第三级， 【dist】的两倍是【8】如果【k=0】 那么【k+8】就会超出最大序号【7】，那么考虑第三级分成【4】组，每一组只有一个蝶形运算，有别于其他两个，可以用for循环这样写
for(k=0;k<len;k+=(dist<<1)) //在二进制左移一位相当于乘2 ，这里len是FFT的点数N同时考虑每一级内又有很多组，可以用for循环嵌套
for(j=0;j<dist;j++)
for(k=0;k<len;k+=(dist<<1)) 最终把分析的结果整合到一起，写成单独的函数 ，注意一下，程序中的一些标号和前面分析的会有点不同
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)); //拷贝到输出的数组
}结果测试
这里我用序列输入到FFT，输出应该是
文章来源：https://www.toymoban.com/news/detail-414009.html
对照结果，和计算一样，程序成功文章来源地址https://www.toymoban.com/news/detail-414009.html
完整程序
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define PI 3.1415926
struct 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));
}
到了这里，关于快速傅里叶变换(FFT)c语言实现的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
