快速傅里叶变换FFT的C程序代码实现
一、彻底理解傅里叶变换
快速傅里叶变换(fast fourier transform)是离散傅里叶变换的一种快速算法，简称fft，通过fft可以将一个信号从时域变换到频域。
模拟信号经过a/d转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。
假设采样频率为fs，采样点数为n，那么fft结果就是一个n点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/n。
举例说明：用1khz的采样频率采样128点，则fft结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 hz。
这个频率点的幅值为：该点复数的模值除以n/2(n=1时是直流分量，其幅值是该点的模值除以n)。
二、傅里叶变换的c语言编程
1、对于快速傅里叶变换fft，第一个要解决的问题就是码位倒序。
假设一个n点的输入序列，那么它的序号二进制数位数就是t=log2n.
码位倒序要解决两个问题：①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。
如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢?利用c语言的移位功能!
程序如下，我不多说，看不懂者智商一定在180以下!
复数类型定义及其运算
#define n 64 //64点
#define log2n 6 //log2n=6
/*复数类型*/
typedef struct
{
float real;
float img;
}complex;
complex xdata x[n]; //输入序列
/*复数加法*/
complex add(complex a,complex b)
{
complex c;
c.real=a.real+b.real;
c.img=a.img+b.img;
return c;
}
/*复数减法*/
complex sub(complex a,complex b)
{
complex c;
c.real=a.real-b.real;
c.img=a.img-b.img;
return c;
}
/*复数乘法*/
complex mul(complex a,complex b)
{
complex c;
c.real=a.real*b.real - a.img*b.img;
c.img=a.real*b.img + a.img*b.real;
return c;
}
/***码位倒序函数***/
void reverse(void)
{
unsigned int i,j,k;
unsigned int t;
complex temp;//临时交换变量
for(i=0;i
{
k=i;//当前第i个序号
j=0;//存储倒序后的序号，先初始化为0
for(t=0;t
{
j=1;//k右移一位，次低位变为最低位
}
if(j>i)//如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换)
{
temp=x;
x=x[j];
x[j]=temp;
}
}
}
2、第二个要解决的问题就是蝶形运算
①第1级(第1列)每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，
第m级蝶形运算，每个蝶形的两节点“距离”l=2m-1。
②对于16点的fft，第1级有16组蝶形，每组有1个蝶形;第2级有4组蝶形，每组有2个蝶形;第3级有2组蝶形，每组有4个蝶形;第4级有1组蝶形，每组有8个蝶形。由此可推出，
对于n点的fft，第m级有n/2l组蝶形，每组有l=2m-1个蝶形。
③旋转因子的确定
以16点fft为例，第m级第k个旋转因子为，其中k=0～2m-1-1，即第m级共有2m-1个旋转因子，根据旋转因子的可约性，，所以第m级第k个旋转因子为，其中k=0～2m-1-1。
为提高fft的运算速度，我们可以事先建立一个旋转因子数组，然后通过查表法来实现。
complex code wn[n]=//旋转因子数组
{ //为节省cpu计算时间，旋转因子采用查表处理
//★根据实际fft的点数n，该表数据需自行修改
//以下结果通过excel自动生成
// wn[k].real=cos(2*pi/n*k);
// wn[k].img=-sin(2*pi/n*k);
{1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
{0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
{0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
{0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
{0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
{-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
{-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
{-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
{-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
{-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
{-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
{-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
{0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
{0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
{0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
{0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}
};
3、算法实现
我们已经知道，n点fft从左到右共有log2n级蝶形，每级有n/2l组，每组有l个。所以fft的c语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算(整个fft共log2n级)，中间层循环完成每一组的蝶形运算(每一级有n/2l组)，最内层循环完成单独1个蝶形运算(每一组有l个)。
/***【快速傅里叶变换】***/
void fft(void)
{
unsigned int i,j,k,l;
complex top,bottom,xw;
reverse(); //码位倒序
for(i=0;i
{ //一级蝶形运算
l=1<
for(j=0;j
{ //一组蝶形运算
for(k=0;k
{ //一个蝶形运算
xw=mul(x[j+k+l],wn[n/(2*l)*k]); //碟间距为l
top=add(x[j+k],xw); //每组的第k个蝶形
bottom=sub(x[j+k],xw);
x[j+k]=top;
x[j+k+l]=bottom;
}
}
}
}
三、fft计算结果验证
随便输入一个64点序列，例如
x[n]={{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0}};
在keil中debug查看数组变量x的fft计算结果并与matlab计算结果进行比对，可以发现非常准确，说明程序编写正确！
我国在动力电池方面获重大突破，研发出一种NCA三元高比能量动力锂电池
苹果为什么要与整个互联网行业为敌的态度来保护用户数据？
余承东表示正在研发可折叠的5G手机 预计明年推出
岚图FREE有增程版，还有纯电版：瞄准用户需求 打造四种智能座舱模式
315曝光的骚扰电话呼出系统＂易科芯＂被查处, 龙芯科:业务已停止
快速傅里叶变换FFT的C程序代码实现
又一波半导体企业正在冲刺科创板
案例分享 | AMR的出现在很大程度上改变了传统的人工搬运方式
一加Concept One揭秘 这款手机到底有多神奇
热仿真加速新产品上市
5G时代，AI应用全面爆发，为半导体供应链注入新动能
为什么锤子新机加入iPhone“陪护”功能？老罗一语道破
FPGA数字图像显示原理与实现设计
基于Xilinx Spartan-6 FPGA加速纹理映射的实现
晶体管放大结构定义和判断
什么是系统架构 为什么要做架构设计
客车研发生产企业宇通客车发布2022第一季度报告
FL型风量仪的概述及基本技术参数
凌科连接器注意事项及产品优点
场效应管在音响电路中的合理应用
快速傅里叶变换(fast fourier transform)是离散傅里叶变换的一种快速算法，简称fft，通过fft可以将一个信号从时域变换到频域。
模拟信号经过a/d转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。
假设采样频率为fs，采样点数为n，那么fft结果就是一个n点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/n。
举例说明：用1khz的采样频率采样128点，则fft结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 hz。
这个频率点的幅值为：该点复数的模值除以n/2(n=1时是直流分量，其幅值是该点的模值除以n)。
二、傅里叶变换的c语言编程
1、对于快速傅里叶变换fft，第一个要解决的问题就是码位倒序。
假设一个n点的输入序列，那么它的序号二进制数位数就是t=log2n.
码位倒序要解决两个问题：①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。
如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢?利用c语言的移位功能!
程序如下，我不多说，看不懂者智商一定在180以下!
复数类型定义及其运算
#define n 64 //64点
#define log2n 6 //log2n=6
/*复数类型*/
typedef struct
{
float real;
float img;
}complex;
complex xdata x[n]; //输入序列
/*复数加法*/
complex add(complex a,complex b)
{
complex c;
c.real=a.real+b.real;
c.img=a.img+b.img;
return c;
}
/*复数减法*/
complex sub(complex a,complex b)
{
complex c;
c.real=a.real-b.real;
c.img=a.img-b.img;
return c;
}
/*复数乘法*/
complex mul(complex a,complex b)
{
complex c;
c.real=a.real*b.real - a.img*b.img;
c.img=a.real*b.img + a.img*b.real;
return c;
}
/***码位倒序函数***/
void reverse(void)
{
unsigned int i,j,k;
unsigned int t;
complex temp;//临时交换变量
for(i=0;i
k=i;//当前第i个序号
j=0;//存储倒序后的序号，先初始化为0
for(t=0;t
j=1;//k右移一位，次低位变为最低位
}
if(j>i)//如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换)
{
temp=x;
x=x[j];
x[j]=temp;
}
}
}
2、第二个要解决的问题就是蝶形运算
①第1级(第1列)每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，
第m级蝶形运算，每个蝶形的两节点“距离”l=2m-1。
②对于16点的fft，第1级有16组蝶形，每组有1个蝶形;第2级有4组蝶形，每组有2个蝶形;第3级有2组蝶形，每组有4个蝶形;第4级有1组蝶形，每组有8个蝶形。由此可推出，
对于n点的fft，第m级有n/2l组蝶形，每组有l=2m-1个蝶形。
③旋转因子的确定
以16点fft为例，第m级第k个旋转因子为，其中k=0～2m-1-1，即第m级共有2m-1个旋转因子，根据旋转因子的可约性，，所以第m级第k个旋转因子为，其中k=0～2m-1-1。
为提高fft的运算速度，我们可以事先建立一个旋转因子数组，然后通过查表法来实现。
complex code wn[n]=//旋转因子数组
{ //为节省cpu计算时间，旋转因子采用查表处理
//★根据实际fft的点数n，该表数据需自行修改
//以下结果通过excel自动生成
// wn[k].real=cos(2*pi/n*k);
// wn[k].img=-sin(2*pi/n*k);
{1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
{0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
{0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
{0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
{0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
{-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
{-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
{-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
{-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
{-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
{-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
{-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
{0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
{0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
{0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
{0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}
};
3、算法实现
我们已经知道，n点fft从左到右共有log2n级蝶形，每级有n/2l组，每组有l个。所以fft的c语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算(整个fft共log2n级)，中间层循环完成每一组的蝶形运算(每一级有n/2l组)，最内层循环完成单独1个蝶形运算(每一组有l个)。
/***【快速傅里叶变换】***/
void fft(void)
{
unsigned int i,j,k,l;
complex top,bottom,xw;
reverse(); //码位倒序
for(i=0;i
l=1<
for(k=0;k
xw=mul(x[j+k+l],wn[n/(2*l)*k]); //碟间距为l
top=add(x[j+k],xw); //每组的第k个蝶形
bottom=sub(x[j+k],xw);
x[j+k]=top;
x[j+k+l]=bottom;
}
}
}
}
三、fft计算结果验证
随便输入一个64点序列，例如
x[n]={{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0}};
在keil中debug查看数组变量x的fft计算结果并与matlab计算结果进行比对，可以发现非常准确，说明程序编写正确！
我国在动力电池方面获重大突破，研发出一种NCA三元高比能量动力锂电池
苹果为什么要与整个互联网行业为敌的态度来保护用户数据？
余承东表示正在研发可折叠的5G手机 预计明年推出
岚图FREE有增程版，还有纯电版：瞄准用户需求 打造四种智能座舱模式
315曝光的骚扰电话呼出系统＂易科芯＂被查处, 龙芯科:业务已停止
快速傅里叶变换FFT的C程序代码实现
又一波半导体企业正在冲刺科创板
案例分享 | AMR的出现在很大程度上改变了传统的人工搬运方式
一加Concept One揭秘 这款手机到底有多神奇
热仿真加速新产品上市
5G时代，AI应用全面爆发，为半导体供应链注入新动能
为什么锤子新机加入iPhone“陪护”功能？老罗一语道破
FPGA数字图像显示原理与实现设计
基于Xilinx Spartan-6 FPGA加速纹理映射的实现
晶体管放大结构定义和判断
什么是系统架构 为什么要做架构设计
客车研发生产企业宇通客车发布2022第一季度报告
FL型风量仪的概述及基本技术参数
凌科连接器注意事项及产品优点
场效应管在音响电路中的合理应用