STM32F407的DSP教程第33章 STM32F407不限制点数FFT实现
Posted Simon223
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了STM32F407的DSP教程第33章 STM32F407不限制点数FFT实现相关的知识，希望对你有一定的参考价值。
完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547
第33章 STM32F407不限制点数FFT实现
本章主要讲解不限制点数FFT的实现。
目录
33.1 初学者重要提示
- 由于ARM DSP库限制最大只能4096点，有点不够用，所以整理了个更大点数的。不限制点数，满足2^n即可，n最小值4， 即16个点的FFT，而最大值不限。
- 此FFT算法没有再使用ARM DSP库，重新实现了一个。
33.2 不限制点数FFT移植
33.2.1 移植FFT相关文件
移植下面两个文件fft.c和FFTInc.h到工程：
33.2.2 添加路径
添加路径，大家根据自己的工程来设置即可：
33.3 不限制点数FFT应用说明
33.3.1 支持的点数范围
支持最小16点FFT，最大值不限，但需满足2^n。
33.3.2 函数InitTableFFT
函数原型：
/*
*********************************************************************************************************
* 函 数 名: 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
}
#endif
函数描述：
这个函数用于FFT计算过程中需要用到的正弦和余弦表。对于8192点和16384点已经专门制作了数值表，存到内部Flash，其它点数继续使用的RAM空间，大家可以根据所使用芯片的RAM和Flash大小，选择正弦和余弦值存到RAM还是Flash。
函数参数：
- 第1个参数是FFT点数。
33.3.3 函数cfft
函数原型：
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) */
}
}
函数描述：
这个函数用于复数FFT变换。
函数参数：
- 第1个参数是复数格式。
- 第2个参数是FFT点数，最小值16，最大值不限，满足满足2^n即可。
33.3.4 FFT幅频响应举例
下面通过函数cff将正弦波做8192点FFT。
/*
*********************************************************************************************************
* 函 数 名: cfft_f32_mag
* 功能说明: 计算幅频
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
static void cfft_f32_mag(void)
{
uint32_t i;
/* 计算一批sin，cos系数 */
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
InitTableFFT(MAX_FFT_N);
#endif
for(i=0; i<MAX_FFT_N; i++)
{
/* 波形是由直流分量，500Hz正弦波组成，波形采样率MAX_FFT_N，初始相位60° */
s[i].real = 1 + cos(2*3.1415926f*500*i/MAX_FFT_N + 3.1415926f/3);
s[i].imag = 0;
}
/* MAX_FFT_N点快速FFT */
cfft(s, MAX_FFT_N);
/* 计算幅频 */
for(i=0; i<MAX_FFT_N; i++)
{
arm_sqrt_f32((float32_t)(s[i].real *s[i].real+ s[i].imag*s[i].imag ), &s[i].real);
}
/* 串口打印求解的幅频 */
for(i=0; i<MAX_FFT_N; i++)
{
printf("%f\\r\\n", s[i].real);
}
}
运行函数cfft_f32_mag可以通过串口打印FFT结果：
从上面的结果中可以出直流分量和正弦波幅值都是没有问题的。
33.3.5 FFT相频响应举例
下面通过函数cff将正弦波做8192点FFT变换：
/*
*********************************************************************************************************
* 函 数 名: PowerPhaseRadians_f32
* 功能说明: 求相位
* 形 参：_usFFTPoints 复数个数，每个复数是两个float32_t数值
* _uiCmpValue 比较值，需要求出相位的数值
* 返 回 值: 无
*********************************************************************************************************
*/
void PowerPhaseRadians_f32(uint16_t _usFFTPoints, float32_t _uiCmpValue)
{
float32_t lX, lY;
uint32_t i;
float32_t phase;
float32_t mag;
for (i=0; i <_usFFTPoints; i++)
{
lX= s[i].real; /* 实部 */
lY= s[i].imag; /* 虚部 */
phase = atan2f(lY, lX); /* atan2求解的结果范围是(-pi, pi], 弧度制 */
arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
if(_uiCmpValue > mag)
{
s[i].real = 0;
}
else
{
s[i].real= phase* 180.0f/3.1415926f; /* 将求解的结果由弧度转换为角度 */
}
}
}
/*
*********************************************************************************************************
* 函 数 名: cfft_f32_phase
* 功能说明: 计算相频
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
static void cfft_f32_phase(void)
{
uint32_t i;
/* 计算一批sin，cos系数 */
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
InitTableFFT(MAX_FFT_N);
#endif
for(i=0; i<MAX_FFT_N; i++)
{
/* 波形是由直流分量，500Hz正弦波组成，波形采样率MAX_FFT_N，初始相位60° */
s[i].real = 1 + cos(2*3.1415926f*500*i/MAX_FFT_N + 3.1415926f/3);
s[i].imag = 0;
}
/* MAX_FFT_N点快速FFT */
cfft(s, MAX_FFT_N);
/* 求相频 */
PowerPhaseRadians_f32(MAX_FFT_N, 0.5f);
/* 串口打印求解相频 */
for(i=0; i<MAX_FFT_N; i++)
{
printf("%f\\r\\n", s[i].real);
}
}
运行函数cfft_f32_phase可以通过串口打印FFT结果：
从上面的结果中可以出计算的初始相位是没有问题的。
33.4 实验例程说明（MDK）
配套例子：
V5-223_不限制点数FFT实现
实验目的：
- 学习不限制点数FFT。
实验内容：
- 启动一个自动重装软件定时器，每100ms翻转一次LED2。
- 按下按键K1，串口打印8192点FFT的幅频响应。
- 按下按键K2，串口打印8192点FFT的相频响应。
使用AC6注意事项
特别注意附件章节C的问题
上电后串口打印的信息：
波特率 115200，数据位 8，奇偶校验位无，停止位 1。
RTT方式打印信息：
程序设计：
系统栈大小分配：
硬件外设初始化
硬件外设的初始化是在 bsp.c 文件实现：
/*
*********************************************************************************************************
* 函 数 名: bsp_Init
* 功能说明: 初始化所有的硬件设备。该函数配置CPU寄存器和外设的寄存器并初始化一些全局变量。只需要调用一次
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void bsp_Init(void)
{
/*
STM32F407 HAL 库初始化，此时系统用的还是F407自带的16MHz，HSI时钟:
- 调用函数HAL_InitTick，初始化滴答时钟中断1ms。
- 设置NVIC优先级分组为4。
*/
HAL_Init();
/*
配置系统时钟到168MHz
- 切换使用HSE。
- 此函数会更新全局变量SystemCoreClock，并重新配置HAL_InitTick。
*/
SystemClock_Config();
/*
Event Recorder：
- 可用于代码执行时间测量，MDK5.25及其以上版本才支持，IAR不支持。
- 默认不开启，如果要使能此选项，务必看V5开发板用户手册第8章
*/
#if Enable_EventRecorder == 1
/* 初始化EventRecorder并开启 */
EventRecorderInitialize(EventRecordAll, 1U);
EventRecorderStart();
#endif
bsp_InitKey(); /* 按键初始化，要放在滴答定时器之前，因为按钮检测是通过滴答定时器扫描 */
bsp_InitTimer(); /* 初始化滴答定时器 */
bsp_InitUart(); /* 初始化串口 */
bsp_InitLed(); /* 初始化LED */
}
主功能：
主程序实现如下操作：
- 启动一个自动重装软件定时器，每100ms翻转一次LED2。
- 按下按键K1，串口打印8192点FFT的幅频响应。
- 按下按键K2，串口打印8192点FFT的相频响应。
/*
*********************************************************************************************************
* 函 数 名: main
* 功能说明: c程序入口
* 形 参: 无
* 返 回 值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
uint8_t ucKeyCode; /* 按键代码 */
bsp_Init(); /* 硬件初始化 */
PrintfLogo(); /* 打印例程信息到串口1 */
PrintfHelp(); /* 打印操作提示信息 */
bsp_StartAutoTimer(0, 100); /* 启动1个100ms的自动重装的定时器 */
/* 进入主程序循环体 */
while (1)
{
bsp_Idle(); /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */
if (bsp_CheckTimer(0)) /* 判断定时器超时时间 */
{
/* 每隔100ms 进来一次 */
bsp_LedToggle(4); /* 翻转LED2的状态 */
}
ucKeyCode = bsp_GetKey(); /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
if (ucKeyCode != KEY_NONE)
{
switch (ucKeyCode)
{
case KEY_DOWN_K1: /* K1键按下 */
cfft_f32_mag();
break;
case KEY_DOWN_K2: /* K2键按下 */
cfft_f32_phase();
break;
default:
/* 其它的键值不处理 */
break;
}
}
}
}
33.5 实验例程说明（IAR）
配套例子：
V5-223_不限制点数FFT实现
实验目的：
- 学习不限制点数FFT。
实验内容：
- 启动一个自动重装软件定时器，每100ms翻转一次LED2。
- 按下按键K1，串口打印8192点FFT的幅频响应。
- 按下按键K2，串口打印8192点FFT的相频响应。
上电后串口打印的信息：
波特率 115200，数据位 8，奇偶校验位无，停止位 1。
RTT方式打印信息：
程序设计：
系统栈大小分配：
硬件外设初始化
硬件外设的初始化是在 bsp.c 文件实现：
/*
*********************************************************************************************************
* 函 数 名: bsp_Init
* 功能说明: 初始化所有的硬件设备。该函数配置CPU寄存器和外设的寄存器并初始化一些全局变量。只需要调用一次
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void bsp_Init(void)
{
/*
STM32F407 HAL 库初始化，此时系统用的还是F407自带的16MHz，HSI时钟:
- 调用函数HAL_InitTick，初始化滴答时钟中断1ms。
- 设置NVIC优先级分组为4。
*/
HAL_Init();
/*
配置系统时钟到168MHz
- 切换使用HSE。
- 此函数会更新全局变量SystemCoreClock，并重新配置HAL_InitTick。
*/
SystemClock_Config();
/*
Event Recorder：
- 可用于代码执行时间测量，MDK5.25及其以上版本才支持，IAR不支持。
- 默认不开启，如果要使能此选项，务必看V5开发板用户手册第8章
*/
#if Enable_EventRecorder == 1
/* 初始化EventRecorder并开启 */
EventRecorderInitialize(EventRecordAll, 1U);
EventRecorderStart();
#endif
bsp_InitKey(); /* 按键初始化，要放在滴答定时器之前，因为按钮检测是通过滴答定时器扫描 */
bsp_InitTimer(); /* 初始化滴答定时器 */
bsp_InitUart(); /* 初始化串口 */
bsp_InitLed(); /* 初始化LED */
}
主功能：
主程序实现如下操作：
- 启动一个自动重装软件定时器，每100ms翻转一次LED2。
- 按下按键K1，串口打印8192点FFT的幅频响应。
- 按下按键K2，串口打印8192点FFT的相频响应。
/*
*********************************************************************************************************
* 函 数 名: main
* 功能说明: c程序入口
* 形 参: 无
* 返 回 值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
uint8_t ucKeyCode; /* 按键代码 */
bsp_Init(); /* 硬件初始化 */
PrintfLogo(); /* 打印例程信息到串口1 */
PrintfHelp(); /* 打印操作提示信息 */
bsp_StartAutoTimer(0, 100); /* 启动1个100ms的自动重装的定时器 */
/* 进入主程序循环体 */
while (1)
{
bsp_Idle(); /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */
if (bsp_CheckTimer(0)) /* 判断定时器超时时间 */
{
/* 每隔100ms 进来一次 */
bsp_LedToggle(4); /* 翻转LED2的状态 */
}
ucKeyCode = bsp_GetKey(); /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
if (ucKeyCode != KEY_NONE)
{
switch (ucKeyCode)
{
case KEY_DOWN_K1: /* K1键按下 */
cfft_f32_mag();
break;
case KEY_DOWN_K2: /* K2键按下 */
cfft_f32_phase();
break;
default:
/* 其它的键值不处理 */
break;
}
}
}
}
33.6 总结
本章节主要设计一个不限制点数的FFT功能，实际项目用到的地方也比较多，望初学着掌握。
以上是关于STM32F407的DSP教程第33章 STM32F407不限制点数FFT实现的主要内容，如果未能解决你的问题，请参考以下文章
STM32F407的DSP教程第30章 STM32F407复数浮点FFT（支持单精度和双精度）
STM32F407的DSP教程第29章 STM32F407移植汇编定点FFT库（64点，256点和1024点）
STM32F407的DSP教程第31章 STM32F407实数浮点FFT（支持单精度和双精度）
STM32F407的DSP教程第50章 STM32F407的样条插补实现，波形拟合丝滑顺畅