支持更高点数的FFT/IFFT
前言
在之前的文章中，介绍了ARM官方提供的FFT函数。尽管该函数对ARM架构处理器做了针对性加速，但最高4096点数的FFT运算稍显羸弱，在更复杂的信号处理领域，这些函数就不能满足需求了。
曾经本人就使用过65536点数的FFT运算暴力提高测量精度。
本文提供8192点数的FFT/IFFT运算算法（理论可至更高）及使用说明。
算法
FFT
unsigned char FFT(Complex *d,int m)
{
#ifndef __EXTERN_W__
static Complex *w;
static int mw = 0;
float arg, w_real, w_imag, wr_real, wr_imag, wtemp;
#endif
static int n = 1;
Complex temp, tmp1, tmp2;
Complex *di, *dip, *dj, *wp;
int i, j, k, l, le, wi;
#ifndef __EXTERN_W__
if(m != mw)
{
if(mw != 0)
free(w);
mw = m;
if(m == 0)
return 0;
n = 1 << m;
le = n >> 1;
w =q;
if(!w)
return 0;
arg = 4.0 * atan(1.0) / le;
wr_real = w_real = cos(arg);
wr_imag = w_imag = -sin(arg);
dj = w;
for(j = 1; j < le; j++)
{
dj->real = (float)wr_real;
dj->imag = (float)wr_imag;
dj++;
wtemp = wr_real * w_real - wr_imag * w_imag;
wr_imag = wr_real * w_imag + wr_imag * w_real;
wr_real = wtemp;
}
}
#else
n = 1 << m;
#endif
le = n;
wi = 1;
for(l = 0; l < m; l++)
{
le >>= 1;
for(i = 0; i < n; i += (le << 1))
{
di = d + i;
dip = di + le;
temp.real = (di->real + dip->real);
temp.imag = (di->imag + dip->imag);
dip->real = (di->real - dip->real);
dip->imag = (di->imag - dip->imag);
*di = temp;
}
wp = (Complex*)w + wi - 1;
for(j = 1; j < le; j++)
{
tmp1 = *wp;
for(i = j; i < n; i += (le << 1))
{
di = d + i;
dip = di + le;
temp.real = (di->real + dip->real);
temp.imag = (di->imag + dip->imag);
tmp2.real = (di->real - dip->real);
tmp2.imag = (di->imag - dip->imag);
dip->real = (tmp2.real * tmp1.real - tmp2.imag * tmp1.imag);
dip->imag = (tmp2.real * tmp1.imag + tmp2.imag * tmp1.real);
*di = temp;
}
wp += wi;
}
wi <<= 1;
}
for(i = 0; i < n; i++)
{
j = 0;
for(k = 0; k < m; k++)
j = (j << 1) | ((i >> k) & 1);
if(i < j)
{
di = d + i;
dj = d + j;
temp = *dj;
*dj = *di;
*di = temp;
}
}
return 1;
}
IFFT
unsigned char IFFT(Complex *d, int m) {
int n = 1 << m; // 计算点数N=2^m
// 1. 对输入数据取共轭
for (int i = 0; i < n; i++) {
d[i].imag = -d[i].imag;
}
// 2. 执行FFT（复用原FFT函数）
if (!FFT(d, m)) {
return 0;
}
// 3. 再次取共轭并除以N
float inv_n = 1.0f / n;
for (int i = 0; i < n; i++) {
d[i].real = d[i].real * inv_n;
d[i].imag = -d[i].imag * inv_n;
}
return 1;
}
使用说明
定义复数数组，实部赋值为原始信号数据，虚部全部赋值为0，调用FFT函数后得到结果。 根据数字信号处理的基本知识，所得结果数组实部与虚部的平方和开方结果即信号幅值信息。 下面提供FFT使用示例
void FFT_Init(float* fundamental_magnitude,float* fundamental_phase)
{
for(int i=0; i<FFT_LENGTH; i++)
{
q1[i].real = adc1_data[i] *3.3/65536;
q1[i].imag = 0;
}
FFT(q1, 13);
/* 计算幅度和相位 */
for(uint32_t i=0; i<FFT_LENGTH/2; i++)
{
// 计算模值
float real = q1[i].real;
float imag = q1[i].imag;
// 缓存中间结果，减少重复计算
float real_square = real * real;
float imag_square = imag * imag;
arm_sqrt_f32(real_square + imag_square, &magnitude11[i]);
phase11[i] = atan2f(imag, real);
}
/* 找出基频成分（幅度最大的非零频率） */
uint32_t fundamental_idx = 0;
float max_mag = 0;
for(uint32_t i=2; i<FFT_LENGTH/2; i++) // 跳过直流分量
{
if(magnitude11[i] > max_mag)
{
max_mag = magnitude11[i];
fundamental_idx = i;
}
}
/* 获取基频的幅度和相位 */
*fundamental_magnitude = magnitude11[fundamental_idx];
*fundamental_phase = phase11[fundamental_idx]; // 弧度
//printf("%d\r\n",fundamental_idx);
//printf("%f\r\n",240000000.0*fundamental_idx/(FFT_LENGTH*146));
}IFFT使用原理相同 代码工程同步上传至github
Thanks for reading!