【STM32F407的DSP教程】第33章 STM32F407不限制点数FFT实现
完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547
第33章 STM32F407不限制点数FFT实现
本章主要讲解不限制点数FFT的实现。
33.1 初学者重要提示
33.2 不限制点数FFT移植
33.3 不限制点数FFT应用说明
33.4 实验例程说明(MDK)
33.5 实验例程说明(IAR)
33.6 总结
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功能，实际项目用到的地方也比较多，望初学着掌握。
微信公众号：armfly_com
安富莱论坛：www.armbbs.cn
安富莱淘宝：https://armfly.taobao.com
