使用ARM单片机进行FIR滤波
概要
在本项目中，为了降低数字信号处理过程中的主控制器成本，我们不得不使用MCU（微控制器）来替代传统的FPGA与DSP。然而对于MCU而言，实时的数字信号处理是一项非常艰巨的任务。因为每执行一次运算都要重复取指、译码、执行等操作，这些操作累加起来将会消耗非常多的时间。但是只要思想不滑坡，方法总比困难多！
举一个最简单的例子，在STM32F103C8T6上，使用定时器来触发进入中断的最小间隔为3us（HAL库）。这就意味着如果需要在这样的中断里面更新DAC的数值，它的最大采样率就仅为333Ksps，远小于手册标称1Msps。但是好在STM32为我们提供了非常多方便的外设可供调用。在上面的例子中，我们可以使用DMA+DAC的架构从数组中直接更新DAC的取值，实测采样率一直到4Msps都可以提供较好的波形。
借助DMA我们实现了传统MCU架构中无法实现的高速信号搬移（外设->内存、内存->外设、内存->内存），但是信号的实时处理仍然是一个很头疼的问题。然而这个问题已经由ARM的工程师解决了，他们在CMSIS-DSP包中提供了非常多的快速的数学运算方法，包括基本数学运算、矩阵运算、滤波、DSP、统计运算等等。直接使用这些最聪明的工程师的成果，会帮助你实现最快的运算速度与最好的兼容性（所有的ARM类型处理器都可使用）。
在本文中，我们将会介绍一个产生信号并采集滤波的例子：
- 调制波输出：使用STM32内置的DAC输出一个50Hz、100Hz与1KHz等幅度叠加的调制信号，使用DMA更新数据，采样率为100Ksps。
- 调制波采集：使用内置ADC1对上述产生的数据通过DMA进行乒乓操作采集处理，并将其通过串口发送到电脑分析，采样率为5Ksps。
- 降采样：数据采集完成后，对其进行5倍的降采样，使其采样率变为1Ksps。
- FIR滤波：再对上述降采样的信号进行一次FIR滤波，滤除200Hz的信号，仅保留50Hz的信号。
- FFT分析：在STM32上运行一个实时的FFT运算，离线分析降采样后的数据在频域表现如何。
本文的开发与测试条件如下：
- 开发板：STM32F103RCT6核心板，板载HSE晶振频率8M，HSI晶振频率32.768KHz。
- 示波器：Hantek 6022BE虚拟示波器，最大采样率48Msps，带宽20MHz。
- 示波器探头：P6200（最大带宽200M），阻抗1M$\Omega$，衰减1X（约7M带宽）。
- 调试器：ST-Link V2，购自淘宝。
- IO接口：模拟输入PA0，模拟输出PA4，LED灯PA1，串口PA9、PA10，使用SWD调试并开启HSE&HSI。
- 开发软件：编辑器Keil V5.35，编译器ARM V5（V6目前对CubeMX生成的代码不是很友好），Cube固件包STM32_F1_V1.8.4。
引脚的Pin-Map如下图所示。
调制波输出
如果我们要对ADC采集到的数据进行滤波，那么我们就应该先准备一份供测试用的模拟信号。如果手边有任意信号发生器的话，那么直接使用其生成包含若干个频率的信号即可，而在这里我们使用STM32内置的DAC输出一个包含50Hz、200Hz与1KHz的调制信号。
如前文所述，通过在定时器的中断函数中更新DAC数值的方法能达到的最高采样率约为333Ksps，这样无法产生一些频率较高的信号，浪费了STM32的DAC性能。因此在这里我们直接使用DMA把数据从内存搬移到外设，搬移的触发条件为TIM4的定时器溢出。对于STM32外设的描述，最准确的来自于STM32中文参考手册，最直观的来自于STM32CubeMX。下图是来自中文参考手册V10的DAC框图（P182），从架构图以及后文描述中就可以得到实现DMA+DAC架构的详细思路。
DAC结构描述
1、结合该框架图与文档后文描述提取关键部分，我们可以得知该DAC的触发可以有8个源：包括软件、定时器与外部中断触发；包含一个DMA请求的通道；数据只能选择对齐方式后送入DHRx，再由控制器送入DORx，用于数字->模拟转换。
2、数据输出：该DAC外设输出模拟值的流程大致为，通过DAC_DHRyyyx的数值更新DHRx，经过内部的控制逻辑更新DORx，根据DORx寄存器内部的数字值去更新模拟值。但是因为STM32是32位单片机，故其内核的寄存器位数均为32位。如果要与DAC外设的12位对应上，那么应该要选择恰当的数据对齐方式。在参考手册的184页详细描述了DAC数据格式，将其分为了：
- 8位数据右对齐：数据写入DAC_DHR8Rx[7:0]；
- 12位数据左对齐：数据写入DAC_DHR8Rx[15:4]；
- 12位数据右对齐：数据写入DAC_DHR8Rx[11:0]。
注意，下图的数据寄存器对齐仅适用于单DAC模式，对于双DAC输出其高位也需要放置数据。
3、触发源选择：在单片机内部触发DAC模拟值更新有很多种方法，分为：定时器TRGO、EXIT外部中断、SWTRIG软件触发。为了实现最快的采样率，我们这里采用定时器触发的方式，触发源的修改需要配置TSELx[2:0]（当然CubeMX已经为我们配置好了）。
3、DMA搬移：F103的外设上一共有两路DAC，均有着自己独立的DMA通道可用于快速的数据更新。在参考手册的148页的DMA2请求映像中也有对DAC通道的描述，从图中我们可以看出DAC-CH1使用到的通道为DMA2的Channel3。
CubeMX配置
前文我们花了比较大的篇幅来介绍DAC的内部架构，这些介绍的目的是为了帮助我们更好的连接这个外设，使得我们能够以一种很友好的方式使用STM32CubeMX配置该外设，而我们只需要编写上层很少的应用代码即可。
1、基础配置：对外部HSE、HSI，时钟树，SWD调试接口，GPIO外设的配置不再赘述。
2、定时器配置：对定时器进行配置的目的是通过其TRGO EVENT来触发DAC，以实现高性能与采样速率可变的DAC输出。在这里参考前文的外部触发源表格，选择TIM4作为触发源。在定时器的配置界面中，有四个关键性的参数需要选择：Prescaler（PSC）—— 预分频系数，这个系数决定着从APB总线上分频到定时器的输入时钟；Counter Period（ARR）—— 自动重装载计数值，每个时钟周期内计时器都会累加器内部的Counter，直到达到数值ARR后溢出产生中断进入下一次计数；auto-reload preload —— 决定改变ARR数值时定时器是否停止当前周期立即响应，我一般选择在Enable；Trigger Event Selection —— 用于选择定时器引起触发的事情选择，选择Update Event便可实现定时器溢出产生更新时间，该事件可用于DAC的触发。
3、DAC配置：在STM32CubeMX配置DAC的过程中，主要有如下关键的内容需要配置：OUT1 Configuration —— 输出配置使能；Output Buffer —— 输出缓冲器使能，使能该选项可以得到类似运放输出的特性，获得较低的输出阻抗方便带载；Trigger —— DAC输出的触发源选择，这里直接选择前文的Timer 4 Trigger Out Event即可。
4、DMA配置：直接在DAC外设的DMA Settings界面进行DMA的配置，通过Add选项就可以添加从内存到外设的通道。添加完成后会出现一个DMA配置信息条，表明使用DMA2 Channel3通道，方向从内存到外设。该配置界面下方有几个可供选择的参数：Mode —— Circular模式表示内存数据读取完成后会环形开始下一轮读取，Normal模式则表示单次读取；Increment Address —— 声明是否自动递增地址；Data Width —— DMA单次数据搬运的位宽，这里设置为Half Word（在STM32中为16位）满足DAC输出12位的需求。
波形生成代码
当然对于调制波形数据的生成，必须依赖于其他工具进行辅助计算，这里使用MATLAB进行计算。计算的代码如下所示，指定采样率100Ksps后设置了3个需要生成的频率：f1——50Hz、f2——200Hz、f1——1000Hz。然后对上述三个频率的波形进行等幅叠加，最后进行12位量化。
需要额外注意的一点是，上述三个频率相加后的波形包络周期会变为50Hz，为了实现最小的内存占用以及无相位跳变需要指定波形数据的长度为$\frac{fs}{f_m}=2000$。
close all;clear;clc;
n=2000;
quantify_bit=12;
fs=100e3;f1=50;f2=200;f3=1000;
t=0:1/fs:(n-1)*(1/fs);
s=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);
s=s/max(abs(s));
quantify_s=round(s*(2^(quantify_bit-1)-1));
quantify_s=quantify_s+(2^(quantify_bit-1)-1);
plot(t,quantify_s);运行上述MATLAB代码，便可以在弹出的图形窗口中看到调制波的一个单周期数据。
因为CubeMX已经为我们编写完了底层的驱动，故在Keil中我们只需要编写很少量的代码便可实现波形的DMA-DAC输出。我们需要进行的操作有：
1、导入MATLAB计算的调制波数据，可以用常数数组表示以节省内存。
//fs=100k,f1=50,f2=100,f3=1k
const uint16_t and_rom[]={2747,2800,2852,2905,2956}; //省略了绝大部分数据2、初始化TIM4定时器，使其以100K的频率触发DAC。
HAL_TIM_Base_Start(&htim4); //TIM4定时器触发DAC的DMA输出，采样率100Ksps3、开启DAC的DMA传输功能，在这里由于对调制波数组进行了一次类型转换并指定数据12位右对齐。
HAL_DAC_Start_DMA(&hdac,DAC_CHANNEL_1,(uint32_t *)and_rom,2000,DAC_ALIGN_12B_R); //通过DMA把数据从内存循环搬移到外设，为了保证周期完整性，其数据长度为2000测试
代码编译完成后下板复位，此时使用示波器可以直接测量到PA4引脚输出的模拟波形。此模拟波形包含等幅度的50Hz、200Hz与1000Hz信号，分别对应着采集信号的可辨最小->最大频率。
查看示波器采样到的波形包络（最小）频率为49.93Hz，周期20ms。
查看示波器采样到的最细波形周期为999us，对应的频率为1.001KHz。
可见上述复频率调制信号的DAC输出功能正常，当然也可以通过把数据导入MATLAB查看频谱的方式来查看频率分量。
调制波采集
为了验证ADC采集的正确性，我们考虑采集上述调制信号后通过串口输出到PC上位机，然后再通过MATLAB进行频谱分析。ADC采样同样是由定时器触发（5KHz），并通过DMA搬移到内存，最后再通过串口输出到PC上位机。
CubeMX配置
使用CubeMX对ADC、TIM、DMA与UART进行初始化，大体上的架构与DAC输出类似。
1、TIM配置：查看ADC的结构框图后选择使用Timer8 Update Event作为触发源，配置PSC为71，ARR为199即可实现一个5KHz的触发信号用于ADC采集。
2、ADC配置：在CubeMX的ADC参数设置窗口有许多可以配置的地方，在本项目中我们只改动External Trigger Conversion Source为Timer 8。其它的选项在此不做改动但略作记录：Mode —— 仅可选独立模式；Data Alignment —— 指定数据的对齐方式，因为C语言中无12bit的数据类型，必须选择如何对齐到现有的数据类型；Continuous Conversion Mode —— 连续转换模式，开启后ADC将会一直进行采集（无论是否有触发）； Number Of Conversion —— 指定转换的通道数，修改后需要使用交错采样来采集多个通道数据。
3、DMA配置：ADC的DMA配置直接在ADC的配置窗口下就可以找到，该DMA通道使用DMA1 Channel 1，模式为Circular，字长为Half Word。
4、UART配置：选择串口模式为异步后，指定波特率、数据长度、校验位与停止位即可。注意在大量数据的发送过程中尽量使用较高波特率，115200Baud对应的极限速率为14.4KB，921600Baud对应的极限速率为115.2KB。
波形采集代码
为了最大程度地提高数据采集与处理效率，我们这里引入了在FPGA中高速数据采集最常见的架构 —— 乒乓操作。即将数据的采集与处理在同一时刻进行，当然实际上处理的是前一段时间的数据，采集的是本时段的数据。在本项目中的乒乓操作使用如下思路进行：
ADC采样率5Ksps，即每秒采集5000个数据。我们将其分为阶段1：采样数1~2500；阶段2：采样数2501~5000，通过HAL库内置的
HAL_ADC_ConvHalfCpltCallback函数在采集到一般时进行跳转。
- 在时刻1，我们采集第1~2500个数据，此时不对数据进行处理。
- 在时刻2，我们采集第2500~5000个数据，此时对时刻1采样的数据进行处理。
- 在时刻3，重新采集1~2500个数据，此时对时刻2采集到的数据进行处理。
- 继续下去，直到采样结束。
可以看出来，每个时刻我们都是对前一时刻的数据进行处理，采集与处理同时进行构成一个简单的流水线架构。
为了实现上述的采集，我们需要在单片机中进行如下步骤：设置状态变量、开启TIM8定时器、开启ADC的DMA采集、设置ADC-DMA采集中断回调（注意这里的中断是由DMA触发ADC）、在主函数中对两个阶段的数据进行处理。
1、设置状态变量与采样缓冲数组，标注当前采样的阶段。
#define ADC_CAP_LENGTH 5000
uint16_t adc1_dat[ADC_CAP_LENGTH];
uint16_t adc1_dma_over_flag = 0;2、开启TIM8定时器，该定时器将会触发ADC以5Ksps的速率采集模拟信号。
HAL_TIM_Base_Start(&htim8); //TIM8定时器触发ADC的DMA采集，采样率5Ksps，注意要关闭ADC的连续采集3、开启ADC的DMA采集，数据将会被送入前文创建的缓冲数组中。
HAL_ADC_Start_DMA(&hadc1,(uint32_t *)adc1_dat,ADC_CAP_LENGTH); //开启ADC的DMA采集，采集数据存入adc1_dat4、设置半中断与全中断回调函数，以adc1_dma_over_flag数值标定采样状态。
void HAL_ADC_ConvHalfCpltCallback(ADC_HandleTypeDef *hadc){
adc1_dma_over_flag = 1;
}
void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc){
adc1_dma_over_flag = 2;
}5、在主函数的死循环中对数据进行处理，发送上个阶段采集到的数据。
while (1){
if(adc1_dma_over_flag == 1){
adc1_dma_over_flag = 0;
for(uint16_t i = 0;i < ADC_CAP_LENGTH/2;i++){
printf("%d\r\n",adc1_dat[i]);
}
}
if(adc1_dma_over_flag == 2){
adc1_dma_over_flag = 0;
for(uint16_t i = ADC_CAP_LENGTH/2;i < ADC_CAP_LENGTH;i++){
printf("%d\r\n",adc1_dat[i]);
}
}
}测试
单片机会把我们的数据打印到串口上，将其复制到MATLAB后便可进行分析。此时，我们的波形采样率为5Ksps，对应的归一化频率为$\frac{f_s}{2}=2.5Ksps$。因此频谱图中的频点0.4对应着1KHz，频点0.08对应着200Hz，频点0.02对应着50Hz。
降采样
前文提到对单片机而言实时信号处理是一个非常艰巨的任务，我们这里ADC的速率为5Ksps的时候一切看起来似乎没什么问题。但是如果采样率达到500Ksps或者更高时，对于单片机的压力会非常大，会出现数据处理延迟的问题。但是很多时候我们并不需要这么高的采样率，而且也无法降低ADC的采样率（受前端抗混叠滤波器影响），这个时候我们就需要考虑进行降采样。
当然正经的降采样一般要求先进行低通滤波避免混叠，再进行抽取。在这里我们仅进行简单的降采样，直接对连续的几个数据求和取均值即可。在这里我们考虑对上述数据进行5倍的降采样，使其采样率变为1Ksps，此时信号的频谱中应当无1KHz的信号。
代码实现
直接在前文的ADC采集代码上略作修改，添加求和取均值的操作即可。
/*宏定义、变量以及数组定义。*/
#define ADC_CAP_LENGTH 5000
#define DDC_RATE 5
uint16_t adc1_dat[ADC_CAP_LENGTH];
uint16_t adc1_ddc_dat[ADC_CAP_LENGTH/DDC_RATE];
uint16_t adc1_dma_over_flag = 0;
/*开启TIM8与ADC-DMA*/
HAL_TIM_Base_Start(&htim8); //TIM8定时器触发ADC的DMA采集，采样率5Ksps，注意要关闭ADC的连续采集
HAL_ADC_Start_DMA(&hadc1,(uint32_t *)adc1_dat,ADC_CAP_LENGTH); //开启ADC的DMA采集，采集数据存入adc1_dat
/*主循环中进行降采样后数组adc1_ddc_dat的打印*/
while (1){
if(adc1_dma_over_flag == 1){
adc1_dma_over_flag = 0;
for(uint16_t i = 0;i < (ADC_CAP_LENGTH/DDC_RATE)/2;i++){
printf("%d\r\n",adc1_ddc_dat[i]);
}
}
if(adc1_dma_over_flag == 2){
adc1_dma_over_flag = 0;
for(uint16_t i = (ADC_CAP_LENGTH/DDC_RATE)/2;i < (ADC_CAP_LENGTH/DDC_RATE);i++){
printf("%d\r\n",adc1_ddc_dat[i]);
}
}
}
/*对采样到的数据每DDC_RATE个进行求和后求平均*/
void HAL_ADC_ConvHalfCpltCallback(ADC_HandleTypeDef *hadc){
uint32_t sum; //每DDC_RATE个数据的和保存至sum
if(ADC1 == hadc->Instance){
for(int i = 0;i < (ADC_CAP_LENGTH/2)/DDC_RATE;i++){
for(int j = 0;j < DDC_RATE;j++){
sum += adc1_dat[i*DDC_RATE+j]; //每DDC_RATE个数据进行累加，数据从i*DDC_RATE~i*DDC_RATE+DDC_RATE
}
adc1_ddc_dat[i] = sum/5; //进行均值滤波
fir_in[i] = adc1_ddc_dat[i] - 2048;
sum = 0;
}
adc1_dma_over_flag = 1;
}
}
void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc){
uint32_t sum;
if(ADC1 == hadc->Instance){
for(int i = (ADC_CAP_LENGTH/2)/DDC_RATE;i < ADC_CAP_LENGTH/DDC_RATE;i++){
for(int j = 0;j < DDC_RATE;j++){
sum += adc1_dat[i*DDC_RATE+j];
}
adc1_ddc_dat[i] = sum/5;
fir_in[i] = adc1_ddc_dat[i] - 2048;
sum = 0;
}
adc1_dma_over_flag = 2;
}
}测试
编译下载后复位，此时窗口将会打印降采样后的波形数据，将其导入MATLAB分析频谱即可。可以看到此时的频谱只在0.1与0.4处有尖峰，因为此时的归一化频率已经为500Hz，故这两个频点对应的频率分别为50Hz与200Hz。
FIR滤波
前文提高我们之所以要进行降采样是因为单片机的计算能力不足，其实还有一个原因是如果在5Ksps采样率下分离50Hz与200Hz非常高品质因数的滤波器，不便于实现。
滤波器设计
因此在这里我们对降采样（1Ksps）得到的50Hz+200Hz信号进行FIR滤波，借助MATLAB使用窗函数设计滤波器的参数为：Hamming窗、采样率1000Hz、截止频率100Hz、滤波器阶数28阶。使用MATLAB设计的滤波器幅频响应如下图所示。可以看到在50Hz附近的幅值约为-0.2dB，在200Hz附近的幅值约为-73dB，可以实现非常好的滤波效果。
将其通过目标->生成C头文件（指定导出数据类型为有符号16位整型，若为F4/F7可以选择单/双度浮点）来导出FIR滤波器系数，导出的文件内容中我们主要关心长度系数BL与滤波器系数组B。
const int BL = 29;
const int16_T B[29] = {
35, 70, 104, 100, 0, -231, -541, -763, -655,
0, 1264, 2967, 4722, 6045, 6537, 6045, 4722, 2967,
1264, 0, -655, -763, -541, -231, 0, 100, 104,
70, 35
};上述的量化步骤中对于M0、M3这类无浮点单元的MCU可以选用有符号整型，量化误差可以通过滤波器设计工具的量化窗口进行比较，看起来整形量化后幅频特性还不错。
滤波实现
在前文中我们使用MATLAB导出了16位有符号量化的滤波器系数，其二进制表达为1位符号位与15位数据位。在STM32上我们可以直接用short来声明16位有符号数，当然这个short在C语言库与DSP库有非常多种等价表达，通过Keil的Go To Definition便可跳转声明。粗略地跳转了几次后，发现16位有符号数有如下几种表述，其均相互等价。
- short：短整型。
- signed short int：与上述类型声明一致。
- int16_t ：在stdint.h文件中对其进行了定义。
- q15_t：在arm_math.h文件中对齐进行了定义。
如需进行FIR滤波的计算，需要进行如下操作：导入CMSIS中的DSP库、系数及结构体的准备、调用FIR初始化函数、调用FIR滤波函数、主函数调用。
1、导入CMSIS中的DSP库，这里使用到的版本为V5.6.0。我们需要使用到的DSP库文件位于...\CMSIS\DSP\Source\，这里根据分类列举了不同功能的DSP文件夹。我们只需要导入每个文件夹下的主要源文件即可，它会自动帮我们调用其它源文件。至于头文件我们只需要包含...\CMSIS\Include这个文件夹，这个文件夹下包含了所有我们需要的头文件。
2、系数与结构体的准备，在进行FIR滤波之前我们需要对一些参数进行定义，包络FIR滤波过程参数、滤波器系数、滤波缓冲区与结构体实例化。其中结构体实例是我们进行初始化与滤波操作的关键参数，其原型为arm_fir_instance_q15。
#define FIR_DATA_LEN 500
#define BLOCK_SIZE 1
#define FIR_COE_LEN 29
#define FIR_TIMES (FIR_DATA_LEN/BLOCK_SIZE)
const q15_t COE[29] = {
35, 70, 104, 100, 0, -231, -541, -763, -655,
0, 1264, 2967, 4722, 6045, 6537, 6045, 4722, 2967,
1264, 0, -655, -763, -541, -231, 0, 100, 104,
70, 35
};
static q15_t FIR_DATA_TEMP[BLOCK_SIZE + FIR_COE_LEN - 1];
arm_fir_instance_q15 S;3、FIR初始化函数调用，直接调用CMSIS中arm_fir_init_q15函数即可。输入的参数依次为：结构体实例指针、滤波器系数长度、滤波器系数指针、缓冲区指针、单块长度。
void fir_q15_init(void){
arm_fir_init_q15(&S, FIR_COE_LEN, &COE[0], &FIR_DATA_TEMP[0], BLOCK_SIZE);
}4、FIR滤波函数调用，每BLOCK_SIZE个数据进行一次FIR滤波。调用arm_fir_q15开始滤波，输入的参数以此为：结构体实例指针、输入数据指针（根据循环递增），输出输出指针（根据循环递增）与单块长度。
void fir_q15(q15_t *in, q15_t *out){
uint16_t i;
for(i = 0; i < FIR_TIMES; i++){
arm_fir_q15(&S, in+(i*BLOCK_SIZE), out+(i*BLOCK_SIZE), BLOCK_SIZE);
}
}5、主函数调用，仍然稍微修改前文的代码，写入初始化与滤波函数。
fir_q15_init(); //初始化16位量化FIR运算
while (1)
{
if(adc1_dma_over_flag == 1){
adc1_dma_over_flag = 0;
fir_q15(&fir_in[0], &fir_out[0]);
for(uint16_t i = 0;i < (ADC_CAP_LENGTH/DDC_RATE)/2;i++){
printf("%d\r\n",fir_out[i]);
}
}
if(adc1_dma_over_flag == 2){
adc1_dma_over_flag = 0;
fir_q15(&fir_in[(ADC_CAP_LENGTH/DDC_RATE)/2], &fir_out[(ADC_CAP_LENGTH/DDC_RATE)/2]);
for(uint16_t i = (ADC_CAP_LENGTH/DDC_RATE)/2;i < (ADC_CAP_LENGTH/DDC_RATE);i++){
printf("%d\r\n",fir_out[i]);
}
}
}
}测试
编译下载上述代码，单片机会打印滤波后的数据。将打印的数据导入MATLAB进行分析，查看频谱得知200Hz信号有60dB的衰减，效果不错。
FFT处理
为了更直观地体现FFT变换的结果，我们选取降采样后的信号作为测试信号。此信号中包含50Hz与200Hz两个频点的信号，理论上在变换后会在频谱图中会出现两个尖峰。
变换实现
若要实现降采样信号的FFT变换计算频谱，需要进行如下几步操作：FFT函数初始化、FFT变换、频谱计算、主函数调用。
1、FFT函数初始化，直接调用arm_rfft_init_q15输入结构体实例即可。该初始化函数的最后两个参数为：ifftFlagR，反变换标识符；bitReverseFlag，比特反转标识符。
arm_rfft_instance_q15 S_fft;
void fft_init_q15(uint16_t len){
arm_rfft_init_q15(&S_fft, len, 0, 1);
}2、FFT变换，调用arm_rfft_q15函数即可。N点RFFT变换实际上是对做N/2点的CFFT变换，因此变换完成后的数据只包含N/2点的复数，并以实部、虚部、实部、虚部等的规则排列，最后的缓冲数组仍包含N个数据。
void fft_q15(q15_t *in, q15_t *out, uint16_t len){
arm_rfft_q15(&S_fft, in, out);
}3、幅值计算，幅值的计算实际上是在每个点计算$\sqrt{real^2+img^2}$的数值，其中输入参数为复数的个数。
void cmplx_mag_q15(q15_t *in, q15_t *out, uint16_t len){
uint32_t blkCnt = len;
q15_t real, imag; /* Temporary input variables */
q31_t acc0, acc1; /* Accumulators */
q63_t sum;
while (blkCnt > 0U)
{
real = *in++;
imag = *in++;
acc0 = ((q31_t) real * real);
acc1 = ((q31_t) imag * imag);
sum = (q63_t)acc0 + acc1;
if(sum >= 32767){
sum = 32767;
}
arm_sqrt_q15((q15_t)sum, out++);
blkCnt--;
}
}4、主函数调用，先进行初始化，再在死循环中调用fft_q15与cmplx_mag_q15进行幅度谱的计算。注意因为频谱的对称性，在这里我们只输出了一半的频谱点。
uint16_t fft_flag = 0;
q15_t fft_out[ADC_CAP_LENGTH/DDC_RATE];
q15_t fft_out_mag[ADC_CAP_LENGTH/DDC_RATE];
fft_init_q15((ADC_CAP_LENGTH/DDC_RATE)/2); //初始化16位量化FFT运算，点数(ADC_CAP_LENGTH/DDC_RATE)/2
while (1)
{
if (adc1_dma_over_flag == 1)
{
adc1_dma_over_flag = 0;
fir_q15(&fir_in[0], &fir_out[0]);
if (fft_flag == 0)
{
fft_flag = 1; //只进行一次FFT运算，点数(ADC_CAP_LENGTH/DDC_RATE)/2
fft_q15(fir_in, fft_out, (ADC_CAP_LENGTH / DDC_RATE) / 2);
cmplx_mag_q15(fft_out, fft_out_mag, (ADC_CAP_LENGTH / DDC_RATE) / 4);
for (uint16_t i = 0; i < (ADC_CAP_LENGTH / DDC_RATE) / 4; i++)
{
printf("%d\r\n", fft_out_mag[i]);
}
}
}
}测试
上述代码会将计算的幅度谱数据输出到串口，通过MATLAB绘制即可在图中看到两个尖峰。因为FFT的点数N为512，所以频率分辨率为$\frac{f_s}{N}=1.95Hz$。两个尖峰点分别位于第27与第103点，对应的频率分别为52.73Hz与200.85Hz。
资料
Ver1.0 2022.11.8 完成文档全部内容的撰写