完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547

第33章 STM32F407不限制点数FFT实现

本章主要讲解不限制点数FFT的实现。

33.1 初学者重要提示

33.2 不限制点数FFT移植

33.2.1 移植FFT相关文件

33.2.2 添加路径

33.3 不限制点数FFT应用说明

33.3.1 支持的点数范围

33.3.2 函数InitTableFFT

33.3.3 函数cfft

33.3.4 FFT幅频响应举例

33.3.5 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 1for (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 */
#elseprintf("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);
#endiffor(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);
#endiffor(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();
#endifbsp_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;}}}
}