目录
一、ST 官方汇编 FFT 库（64点， 256 点和 1024 点）
1、cr4_fft_xxx_stm32
2、计算幅频响应
3、计算相频响应
二、复数浮点 FFT、IFFT（支持单精度和双精度）
1、基础支持
2、单精度函数 arm_cfft_f32
3、双精度函数 arm_cfft_f64
三、实数浮点 FFT（支持单精度和双精度）
1、基础支持
2、单精度函数 arm_rfft_fast_f32
3、双精度函数 arm_rfft_fast_f64
四、不限制点数 FFT 实现
1、函数 InitTableFFT
2、函数 cfft
一、ST 官方汇编 FFT 库（64点， 256 点和 1024 点）
这个汇编的FFT库是来自STM32F10x DSP library，由于是汇编实现的，而且是基4算法，所以实现FFT在速度上比较快。
1、cr4_fft_xxx_stm32
其实现32位定点数的运算，其中高16位为虚部，低16位为实部。
| cr4_fft_1024_stm32 | 函数用于实现1024点数据的FFT计算 |
| cr4_fft_256_stm32 | 函数用于实现256点数据的FFT计算 |
| cr4_fft_64_stm32 | 函数用于实现64点数据的FFT计算 |
2、计算幅频响应
uint32_t input[1024], output[1024], Mag[1024];/* 输入，输出和幅值 */
float32_t Phase[1024]; /* 相位*/
/*
*********************************************************************************************************
* 函 数 名: PowerMag
* 功能说明: 求模值
* 形 参： _usFFTPoints FFT点数
* 返 回 值: 无
*********************************************************************************************************
*/
void PowerMag(uint16_t _usFFTPoints)
{
int16_t lX,lY;
uint16_t i;
float32_t mag;
/* 计算幅值 */
for (i=0; i < _usFFTPoints; i++)
{
lX= (output[i]<<16)>>16; /* 实部*/
lY= (output[i]>> 16); /* 虚部 */
arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
Mag[i]= mag*2; /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
}
/* 由于上面多乘了2，所以这里直流分量要除以2 */
Mag[0]= Mag[0]>>1;
}3、计算相频响应
/*
*********************************************************************************************************
* 函 数 名: Power_Phase_Radians
* 功能说明: 求相位
* 形 参： _usFFTPoints FFT点数， uiCmpValue 阀值
* 返 回 值: 无
*********************************************************************************************************
*/
void Power_Phase_Radians(uint16_t _usFFTPoints, uint32_t _uiCmpValue)
{
int16_t lX, lY;
uint16_t i;
float32_t phase;
float32_t mag;
for (i=0; i <_usFFTPoints; i++)
{
lX= (output[i]<<16)>>16; /* 实部 */
lY= (output[i] >> 16); /* 虚部 */
phase = atan2(lY, lX); /* atan2求解的结果范围是(-pi, pi], 弧度制 */
arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
if(_uiCmpValue > mag)
{
Phase[i] = 0;
}
else
{
Phase[i] = phase* 180.0f/3.1415926f; /* 将求解的结果由弧度转换为角度 */
}
}
二、复数浮点 FFT、IFFT（支持单精度和双精度）
新版 DSP 库浮点 FFT 推荐使用混合基函数 arm_cfft_f32，而基 2 函数 arm_cfft_radix2_f32 和基 4函数 arm_cfft_radix4_f32 将废弃。 ARM 说明如下：
DSP 库的早期发行版提供了单独的 radix-2 和 radix-4 对浮点数据进行运算的算法。 这些功能仍然提供，但已弃用。 相比新版函数，老版的功能较慢且通用性较低
1、基础支持
当前复数 FFT 函数支持三种数据类型，分别是浮点，定点 Q31 和 Q15。这些 FFT 函数有一个共同的特点，就是用于输入信号的缓冲，在转化结束后用来存储输出结果。这样做的好处是节省了 RAM 空间，不需要为输入和输出结果分别设置缓存。由于是复数 FFT，所以输入和输出缓存要存储实部和虚部。存储顺序如下： {real[0], imag[0], real[1], imag[1],………………} ，在使用中切记不要搞错。
2、单精度函数 arm_cfft_f32
支持16、32、64、128、256、512、1024、2048、4096点单精度复数FFT、IFFT。
3、双精度函数 arm_cfft_f64
支持16、32、64、128、256、512、1024、2048、4096点双精度复数FFT、IFFT
三、实数浮点 FFT（支持单精度和双精度）
CMSIS DSP 库里面包含一个专门用于计算实数序列的 FFT 库，很多情况下，用户只需要计算实数序列即可。计算同样点数 FFT 的实数序列要比计算同样点数的虚数序列有速度上的优势。
快速的 rfft 算法是基于混合基 cfft 算法实现的。
1、基础支持
一个 N 点的实数序列 FFT 正变换采用下面的步骤实现：
由上面的框图可以看出，实数序列的 FFT 是先计算 N/2 个实数的 CFFT，然后再重塑数据进行处理从而获得半个 FFT 频谱即可（利用了 FFT 变换后频谱的对称性）。
一个 N 点的实数序列 FFT 逆变换采用下面的步骤实现：
实数 FFT 支持浮点， Q31 和 Q15 三种数据类型。
2、单精度函数 arm_rfft_fast_f32
支持32、64、128、256、512、1024、2048、4096点单精度实数FFT、IFFT。
3、双精度函数 arm_rfft_fast_f64
支持32、64、128、256、512、1024、2048、4096点双精度实数FFT、IFFT
四、不限制点数 FFT 实现
可以看到前面的由于 ARM DSP 库限制最大只能 4096 点，有点不够用，所以整理了个更大点数的。不限制点数，满足 2^n 即可， n 最小值 4， 即 16 个点的 FFT，而最大值不限。
此 FFT 算法没有再使用 ARM DSP 库，重新实现了一个。
1、函数 InitTableFFT
这个函数用于 FFT 计算过程中需要用到的正弦和余弦表。 对于 8192 点和 16384 点已经专门制作了数值表，存到内部 Flash，其它点数继续使用的 RAM 空间。
/*
*********************************************************************************************************
* 函 数 名: Int_FFT_TAB
* 功能说明: 正弦和余弦表
* 形 参: 点数
* 返 回 值: 无
*********************************************************************************************************
*/
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
float32_t costab[MAX_FFT_N/2];
float32_t sintab[MAX_FFT_N/2];
void InitTableFFT(uint32_t n)
{
uint32_t i;
/* 正常使用下面获取 cos 和 sin 值 */
#if 1
for (i = 0; i < n; i ++ )
{
sintab[ i ]= sin( 2 * PI * i / MAX_FFT_N );
costab[ i ]= cos( 2 * PI * i / MAX_FFT_N );
}
/* 打印出来是方便整理 cos 值和 sin 值数组，将其放到内部 Flash，从而节省 RAM */
#else
printf("const float32_t sintab[%d] = {\r\n", n);
for (i = 0; i < n; i ++ )
{
sintab[ i ]= sin( 2 * PI * i / MAX_FFT_N );
printf("%.11ff,\r\n", sintab[ i ]);
}
printf("};\r\n");
printf("const float32_t costab[%d] = {\r\n", n);
for (i = 0; i < n; i ++ )
{
sintab[ i ]= cos( 2 * PI * i / MAX_FFT_N );
printf("%.11ff,\r\n", sintab[ i ]);
}
printf("};\r\n");
#endif
}
#endif2、函数 cfft
这个函数用于复数 FFT 变换，如果需要实数则直接把虚部设置为0即可。
void cfft(struct compx *_ptr, uint32_t FFT_N )
/*
*********************************************************************************************************
* 函 数 名: cfft
* 功能说明: 对输入的复数组进行快速傅里叶变换（FFT）
* 形 参: *_ptr 复数结构体组的首地址指针 struct 型
* FFT_N 表示点数
* 返 回 值: 无
*********************************************************************************************************
*/
void cfft(struct compx *_ptr, uint32_t FFT_N )
{
float32_t TempReal1, TempImag1, TempReal2, TempImag2;
uint32_t k,i,j,z;
uint32_t Butterfly_NoPerColumn; /* 每级蝶形的蝶形组数 */
uint32_t Butterfly_NoOfGroup; /* 蝶形组的第几组 */
uint32_t Butterfly_NoPerGroup; /* 蝶形组的第几个蝶形 */
uint32_t ButterflyIndex1,ButterflyIndex2,P,J;
uint32_t L;
uint32_t M;
z=FFT_N/2; /* 变址运算，即把自然顺序变成倒位序，采用雷德算法 */
for(i=0,j=0;i<FFT_N-1;i++)
{
/*
如果 i<j,即进行变址 i=j 说明是它本身， i>j 说明前面已经变换过了，不许再变化，注意这里一般是实数有虚数部分 设置成结合体
*/
if(i<j)
{
TempReal1 = _ptr[j].real;
_ptr[j].real= _ptr[i].real;
_ptr[i].real= TempReal1;
}
k=z; /*求 j 的下一个倒位序 */
while(k<=j) /* 如果 k<=j,表示 j 的最高位为 1 */
{
j=j-k; /* 把最高位变成 0 */
k=k/2; /* k/2，比较次高位，依次类推，逐个比较，直到某个位为 0，通过下面那句 j=j+k 使其变为 1 */
}
j=j+k; /* 求下一个反序号，如果是 0，则把 0 改为 1 */
}
/* 第 L 级蝶形(M)第 Butterfly_NoOfGroup 组(Butterfly_NoPerColumn)第 J 个蝶形(Butterfly_NoPerGroup)****** */
/* 蝶形的组数以 2 的倍数递减 Butterfly_NoPerColumn，每组中蝶形的个数以 2 的倍数递增 Butterfly_NoPerGroup */
/* 在计算蝶形时，每 L 列的蝶形组数,一共有 M 列，每组蝶形中蝶形的个数,蝶形的阶数(0,1,2.....M-1) */
Butterfly_NoPerColumn = FFT_N;
Butterfly_NoPerGroup = 1;
M = log2(FFT_N);
for ( L = 0;L < M; L++ )
{
Butterfly_NoPerColumn >>= 1; /* 蝶形组数 假如 N=8，则(4,2,1) */
//第 L 级蝶形 第 Butterfly_NoOfGroup 组（0,1， ....Butterfly_NoOfGroup-1）
for ( Butterfly_NoOfGroup = 0;Butterfly_NoOfGroup < Butterfly_NoPerColumn;Butterfly_NoOfGroup++ )
{
/* 第 Butterfly_NoOfGroup 组中的第 J 个 */
for ( J = 0;J < Butterfly_NoPerGroup;J ++ )
{ /* 第 ButterflyIndex1 和第 ButterflyIndex2 个元素作蝶形运算,WNC */
/* (0,2,4,6)(0,1,4,5)(0,1,2,3) */
/* 两个要做蝶形运算的数相距 Butterfly_NoPerGroup， ge=1,2,4 */
/* 这里相当于 P=J*2^(M-L),做了一个换算下标都是 N (0,0,0,0)(0,2,0,2)(0,1,2,3) */
ButterflyIndex1 = ( ( Butterfly_NoOfGroup * Butterfly_NoPerGroup ) << 1 ) + J;
ButterflyIndex2 = ButterflyIndex1 + Butterfly_NoPerGroup;
P = J * Butterfly_NoPerColumn;
/* 计算和转换因子乘积 */
TempReal2 = _ptr[ButterflyIndex2].real * costab[ P ] + _ptr[ButterflyIndex2].imag *
sintab[ P ];
TempImag2 = _ptr[ButterflyIndex2].imag * costab[ P ] - _ptr[ButterflyIndex2].real *
sintab[ P ] ;
TempReal1 = _ptr[ButterflyIndex1].real;
TempImag1 = _ptr[ButterflyIndex1].imag;
/* 蝶形运算 */
_ptr[ButterflyIndex1].real = TempReal1 + TempReal2;
_ptr[ButterflyIndex1].imag = TempImag1 + TempImag2;
_ptr[ButterflyIndex2].real = TempReal1 - TempReal2;
_ptr[ButterflyIndex2].imag = TempImag1 - TempImag2;
}
}
Butterfly_NoPerGroup<<=1; /* 一组中蝶形的个数(1,2,4) */
}
