// 快速傅里 叶变换FFT的C语言算法彻底研究
// LED音乐频谱显示的核心算法就是快速傅里叶变换，FFT的理解和编程还是比较难的，特地撰写此文分享一下研究成果。
// 一、彻底理解傅里叶变换
// 快速傅里叶变换（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位二进制数倒序；②将倒序后的两个存储单元进行交换。

// ①将t=3位二进制数倒序

// ②将倒序后的两个存储单元进行交换
// 如果输入序列的自然顺序号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.realb.real - a.imgb.img;
c.img=a.realb.img + a.imgb.real;
return c;
}
/码位倒序函数/
void Reverse(void)
{
unsigned int i,j,k;
unsigned int t;
complex temp;//临时交换变量
for(i=0;i<N;i++)//从第0个序号到第N-1个序号
{
k=i;//当前第i个序号
j=0;//存储倒序后的序号，先初始化为0
for(t=0;t<log2N;t++)//共移位t次，其中log2N是事先宏定义算好的
{
j<<=1;
j|=(k&1);//j左移一位然后加上k的最低位
k>>=1;//k右移一位，次低位变为最低位
}
if(j>i)//如果倒序后大于原序数，就将两个存储单元进行交换（判断j>i是为了防止重复交换）
{
temp=x[ i];
x[ i]=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(2PI/Nk);
// WN[k].img=-sin(2PI/Nk);
{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<log2N;i++) /共log2N级/
{ //一级蝶形运算
l=1<<i;//l等于2的i次方
for(j=0;j<N;j+=2l) /每L个蝶形是一组，每级有N/2L组/
{ //一组蝶形运算
for(k=0;k<l;k++) /每组有L个/
{ //一个蝶形运算
xW=mul(x[j+k+l],WN[N/(2l)*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计算结果进行比对，可以发现非常准确，说明程序编写正确！

假设采样频率为 Fs，信号频率 F，采样点数为 N。那么 FFT 之后结果就是一个为 N 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关
系呢？假设原始信号的峰值为 A，那么 FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是 A的 N/2 倍。 而第一个点就是直流分量，它的模值就是直流分量的 N 倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即 0Hz），而最后一个点 N 的再下一个点（实际上这个点是不存在的，这里是假设的第 N+1 个点，可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率 Fs，这中间被 N-1 个点平均分成 N 等份，每个点的频率依次增加。例如某点 n 所表示的频率为：F n =(n − 1)xF s
根据 FFT 结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。总的来说，这个过程就是这样：假设采样频率为 Fs，采样点数为 N，做 FFT 之后，某一点 n（n 从 1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以 N/2 就是对应该频率下的信号的幅度（对于直流信号是除以 N）；该点的相位即是对应该频率下的信号的相位。相位的计算可用函数 atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi 到 pi。要精确到 xHz，则需要采样长度为 1/x 秒的信号，并做 FFT。要提高频率分辨率，就需要增加采样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的 0，使其长度达到需要的点数，再做 FFT，这在一定程度上能够提高频率分辨力。具体的频率细分法大家可参考相关文献。
频谱图像显示的就是信号的组成正弦波的峰峰值（即模值）和各自组成信号的频率。
以下图片是方波的粗略频谱。

ThinkPadW530的STM32傅里叶函数.c文件路径D:\keil5\ARM\Pack\ARM\CMSIS\4.3.0\CMSIS\DSP_Lib\Examples
已经调试成功的傅里叶变换程序代码如下：
*uint32_t input[1024], output[1024], Mag[1024]; / 输入，输出和幅值 /
/

---

• 函 数 名: PowerMag
• 功能说明: 求模值
• 形 参：_usFFTPoints FFT点数
• 返 回 值: 无

---

/
void PowerMag(uint16_t _usFFTPoints)
{
int16_t lX,lY;
uint16_t i;
float32_t mag;
/ 计算幅值 /
for (i=0; i < _usFFTPoints/2; i++)
{
lX= (output[i]<<16)>>16; / 实部*/
lY= (output[i]>> 16); /* 虚部 /
{
float X = _usFFTPoints((float)lX)/32768;
float Y = _usFFTPoints*((float)lY)/32768;
mag = sqrt(XX+ YY)/_usFFTPoints; /* 求模 /
Mag[i]=(u32)(mag65536); /* 求模后乘以2才是实际模值，直流分量不需要乘2 /
}
}
/ 由于上面多乘了2，所以这里直流分量要除以2 /
Mag[0] = Mag[0]/2;
}
/

---

• 函 数 名: DSP_FFT1024
• 功能说明: 1024点FFT实现
• 形 参：无
• 返 回 值: 无

---

/
//void PowerMag(uint16_t _usFFTPoints);
void DSP_FFT1024(void)
{
uint16_t i;
/ 获得1024个采样点 /
for (i = 0; i < 1024; i++)
{
// input[i] = 0;
/ 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率1KHz /
input[i] = 1024 + 1024sin(23.1415926f50i/1000) + 512sin(23.1415926f20i/1000)+ 2048sin(23.1415926f100i/1000) ;
}
/ 计算1024点FFT
output：输出结果，高16位是虚部，低16位是实部。
input ：输入数据，高16位是虚部，低16位是实部。
第三个参数必须是1024。
/
cr4_fft_1024_stm32(output, input, 1024);
/ 求幅值 /
PowerMag(1024);
/ 打印输出结果 /
for (i = 0; i < 1024; i++)
{
printf("%4d, %.2f %10d\n",i,(((float)i1024/1024)-1),Mag[i]);
}
//(((float)i1024/1024)-1)该项为频率计算，用i采样频率（1024）/（采样点数）
}**
找到以上代码的网址：https://wenku.baidu.com/view/08ccee0984868762cbaed532.html