ARM DSP库CMSIS-DSP的使用——以STM32F4浮点FFT为例 [原创www.cnblogs.com/helesheng]

之前用过STM32F10x比较多，做数字信号处理用过意法半导体官方的STM32F10x_DSP_Lib_V2.0.0，总觉得这个库不太好用：

1、数字滤波器FIR和IIR的函数只能对存在缓冲中的数组滤波，且没有能够保存滤波器中间状态的数据结构，导致再次调用这些滤波函数时需要一个新的稳定期，无法实现连续实时滤波。

2、只有定点FFT函数，输出结果保存在16bits的定点寄存器中，大大降低了FFT结果的动态范围。

因此我在STM32F10x上总习惯采用自己编写的C语言数字信号处理函数，无论是时间效率还是内存效率都差强人意。

最近在项目中开始使用带有浮点处理单元（FPU）的Cortex-M4核STM32F4，发现ARM公司自己提供的CMSIS(Cortex Microcontroller Software Interface Standard)中的数字信号库（CMSIS-DSP）确实要给力得多。下面以需要浮点能力的FFT功能为例，介绍CMSIS-DSP在MDK环境中的使用。

以下原创内容欢迎网友转载，但请注明出处： https://www.cnblogs.com/helesheng

一、准备工作和准备知识

1、在意法半导体官方网站下载包含CMSIS-DSP的外设库（https://www.st.com/en/embedded-software/stsw-stm32065.html）。数字信号库包含在路径STM32F4xx_DSP_StdPeriph_Lib_V1.8.0\Libraries\CMSIS\DSP_Lib和STM32F4xx_DSP_StdPeriph_Lib_V1.8.0\Libraries\CMSIS\Lib下，将它们拷贝到目标工程文件夹下。其中Lib文件夹中包含的是经过不同编译器编译后能够运行在Cortex-M4内核上的底层数学库，DSP_Lib文件夹中包含的是调用底层函数封装而成的API函数源码。

2、Lib文件夹中包含的底层库包括：

arm_cortexM4lf_math.lib
    arm_cortexM4bf_math.lib
    arm_cortexM4l_math.lib
    arm_cortexM4b_math.lib
    arm_cortexM3l_math.lib
    arm_cortexM3b_math.lib
    arm_cortexM0l_math.lib
    arm_cortexM0b_math.lib

M0、M3、M4分别代表Cortex-M0、Cortex-M3、Cortex-M4内核上运行的库；字母l代表小端地址格式，b代表大端地址格式（所有STM32皆为小端地址格式处理器）；字母f代表使用浮点处理单元(FPU)进行计算的库。例如，我打算使用STM32F4处理器的FPU进行FFT运算，就需要将库文件arm_cortexM4lf_math.lib添加到Keil MDK工程中。

3、DSP_Lib中包含了如下所示的多个文件夹，它们中的文件功能如下：1）BasicMathFunctions

基本数学函数：提供浮点数的各种基本运算函数，如向量加减乘除等运算。

2）CommonTables

数字信号处理常用参数表。

3）ComplexMathFunctions

复数计算数学函数。

4）ControllerFunctions

控制算法函数。包括正弦余弦，PID电机控制，矢量Clarke变换，矢量Clarke逆变换等。

5）FastMathFunctions

常见快速算法的数学函数。

6）FilteringFunctions

滤波函数功能，主要为FIR和LMS（最小均方根）等滤波函数。

7）MatrixFunctions

矩阵处理函数。包括矩阵加法、矩阵初始化、矩阵反、矩阵乘法、矩阵规模、矩阵减法、矩阵转置等函数。

8）StatisticsFunctions

统计功能函数。如求平均值、最大值、最小值、计算均方根RMS、计算方差/标准差等。

9）SupportFunctions

支持功能函数，如数据拷贝，Q格式和浮点格式相互转换，Q任意格式相互转换。

10）TransformFunctions

变换功能。包括复数FFT（CFFT）/复数FFT逆运算（CIFFT）、实数FFT（RFFT）/实数FFT逆运算（RIFFT）、和DCT（离散余弦变换）和配套的初始化函数。

例如，这里我打算通过浮点FFT运算进行频谱分析，就使用了TransformFunctions文件夹下的arm_cfft_f32.c中的函数。

二、配置Keil MDK工程

1、在Options的Target选项卡中使能浮点运算单元FPU。

图1

 2、在Options的C/C++选项卡的宏定义中添加如下数学计算可能用到的宏ARM_MATH_CM4,__CC_ARM,ARM_MATH_MATRIX_CHECK,ARM_MATH_ROUNDING。

图2

 三、调用DSP_LIB中的函数实现FFT

1、从ARM官方帮助文档中查找合适的函数及其使用方法

ARM官方的CMSIS-DSP库的帮助文档是HTML格式的网页，保存在.. \STM32F4xx_DSP_StdPeriph_Lib_V1.8.0\Libraries\CMSIS路径下，打开后如下图所示。

图3

我选择32位浮点（float）数据类型的函数arm_cfft_f32();来实现FFT，其原型如网页右侧视图所示。

从上面的网页中可以查得： 

1 void arm_cfft_f32    (
2 const arm_cfft_instance_f32 *     S,
3 float32_t *     p1,
4 uint8_t     ifftFlag,
5 uint8_t     bitReverseFlag 
6 );

函数名称中的字母“c”表明，这个函数的输入和输出都是复数（complex number）。输入复数向量存放在第二个参数p1指向的浮点数组中，函数执行完后输出复数向量亦存放在相同地址。输入输出复数数组的存储结构相同，偶数地址用于存储复数实部，随后的奇数地址用于存储同一个复数的虚部。

第一个参数则是一个指向复数浮点FFT常数结构体的指针，该结构体已经定义在头文件arm_const_structs.h中，可以直接通过&arm_cfft_sR_f32_len256、&arm_cfft_sR_f32_len512、&arm_cfft_sR_f32_len1024等形式应用，FFT变换的点数由改参数决定。

第三个参数ifftFlag用于指定函数执行的FFT逆变换（ifftFlag=1），还是FFT正变换（ifftFlag=0）。第四个参数bitReverseFlag用于指定是否对结果进行位翻转变换，根据FFT理论知识可知，对自然顺序输入进行FFT蝶形运算的结果的顺序，是自然顺序输出的地址高低位翻转的结果。当bitReverseFlag=1时，该函数能够对输出的结果执行地址高低位翻转和顺序置换；当bitReverseFlag=0时，该功能被禁止。大多数情况下，为得到频率依次递增的频域数据，需要使能该函数的位翻转功能。

2、调用函数arm_cfft_f32();实现FFT

未验证STM32F4xx_DSP_StdPeriph_Lib中DSP功能调用的正确性，我对STM32F401RC片上ADC采集得到的信号进行了FFT变换和显示。具体实现步骤如下：

1）通过DMA控制ADC实现高速采样，代码如下（其中DMA配置函数ADC_Config();不是本文重点，具体代码此处未列出）：

1 ADC_Config();//ADC和DMA配置
2 ADC_SoftwareStartConv(ADCx);    //启动adc和DMA传输
3 while(DMA_GetFlagStatus(DMA_STREAMx, DMA_FLAG_TCIF0) == RESET) ;
4 //等待采样和DMA传输完成            
5 ADC_DMACmd(ADCx, DISABLE);
6 DMA_ClearFlag(DMA_STREAMx, DMA_FLAG_TCIF0);

2）将采样结果放置在复数结构的缓冲区中。其中数组uhADCxConvertedValue[]中存放的A/D采样的结果是DMA控制器在此之前存入的，它们被放置到即将进行FFT的浮点数组fft_buf_float[ ]的偶数地址中。而奇数地址中存放的是FFT输入的虚部，这里被全部置0。

1 //FFT函数的输入和输出都是复数，因此还有虚部，将输入填入实部，虚部为0
2 for(i = 0; i<256; i++){
3     fft_buf_float[2*i] = uhADCxConvertedValue[i];
4     fft_buf_float[2*i + 1] = 0;
5 }    

3）由于STM32的ADC是单极性的（只有正数结果），其输入一定含有直流分量，其能量将远大于交流分量，影响FFT结果显示，因此在进行FFT之前先减去结果中的直流分量。（不去除直流分量的FFT结果，读者可以参见本博文后续实验结果部分所示）

1 //减去平均值，去除输入信号的直流成分（输入信号在浮点数组fft_buf_float[ ]的偶数地址中）
2         for(i = 0; i<256; i++)
3             sum = sum + fft_buf_float[2*i];
4         average = sum/256;
5         for(i = 0; i<256; i++)
6             fft_buf_float[2*i] = fft_buf_float[2*i] - average;

4）对需要进行FFT计算的时域数据进行加窗

加窗的目的，是为了防止时域截断造成的能量泄露现象。（对不加窗数据进行FFT的结果，读者可以参见本博文后续实验结果部分所示）我使用了hanning窗，代码如下：

1 //加窗防止能量泄露
2 for(i = 0; i<256; i++)
3     fft_buf_float[2*i] = fft_buf_float[2*i] * hanning_win_table[i] ;    

我没有在CMSIS-DSP库中找到窗函数，好在用Matlab计算浮点汉宁窗数值并不复杂，我将数组hanning_win_table[i]的数值罗列于此。

1 //256点汉宁窗
2 const float hanning_win_table[256] = 
3     {0.000149421078453171, 0.000597595007177931, 0.00134425391964726, 0.00238895154954133,        0.00373106349747415, 0.00536978760418699, 0.00730414442998667, 0.00953297784014107,         0.0120549556958828, 0.0148685706506078, 0.0179721410507924, 0.0213638119410917,         0.0250415561730168, 0.0290031756165303, 0.0332463024738333, 0.0377684006945618,         0.0425667674915437, 0.0476385349562126, 0.0529806717727114, 0.0585899850296627,         0.0644631221285217, 0.0705965727873714, 0.0769866711389634, 0.0836295979217494,         0.0905213827625935, 0.0976579065498020, 0.105034903895052, 0.112647965682748,         0.120492541705279, 0.128563943382607, 0.136857346564560, 0.145367794414147, 0.154090200370186,         0.163019351187458, 0.172149910052584, 0.181476419773753, 0.190993306042403, 0.200694880764894,         0.210575345462196, 0.220628794735546, 0.230849219796014, 0.241230512055860, 0.251766466779543,         0.262450786792195, 0.273277086243339, 0.284238894423618, 0.295329659632231, 0.306542753092784,         0.317871472915207, 0.329309048101366, 0.340848642591985, 0.352483359352449, 0.364206244495054,         0.376010291435238, 0.387888445079305, 0.399833606041146, 0.411838634885427, 0.423896356394721,         0.435999563858022, 0.448141023378083, 0.460313478195001, 0.472509653023471, 0.484722258401111,         0.496943995045254, 0.509167558215622, 0.521385642080249, 0.533590944082063, 0.545776169303514,         0.557934034826625, 0.570057274085885, 0.582138641211355, 0.594170915359415, 0.606146905028547,         0.618059452357584, 0.629901437403849, 0.641665782398637, 0.653345455977481, 0.664933477382692,         0.676422920635650, 0.687806918676346, 0.699078667467724, 0.710231430062342, 0.721258540628942,         0.732153408436510, 0.742909521793458, 0.753520451939554, 0.763979856888295, 0.774281485217411,         0.784419179805244, 0.794386881510759, 0.804178632795004, 0.813788581281829, 0.823210983255770,         0.832440207094966, 0.841470736637101, 0.850297174476322, 0.858914245189185, 0.867316798487695,         0.875499812297549, 0.883458395759753, 0.891187792153812, 0.898683381740745, 0.905940684524234,         0.912955362928245, 0.919723224389528, 0.926240223863451, 0.932502466241655, 0.938506208680101,         0.944247862836110, 0.949723997013056, 0.954931338211443, 0.959866774085119, 0.964527354801480,         0.968910294804540, 0.973012974479810, 0.976832941720003, 0.980367913390621, 0.983615776694547,         0.986574590434830, 0.989242586174910, 0.991618169295584, 0.993699919948085, 0.995486593902703,         0.996977123292440, 0.998170617251262, 0.999066362446550, 0.999663823505453, 0.999962643334866,         0.999962643334866, 0.999663823505453, 0.999066362446550, 0.998170617251262, 0.996977123292440,         0.995486593902703, 0.993699919948085, 0.991618169295584, 0.989242586174910, 0.986574590434830,         0.983615776694547, 0.980367913390621, 0.976832941720003, 0.973012974479810, 0.968910294804540,         0.964527354801480, 0.959866774085119, 0.954931338211443, 0.949723997013056, 0.944247862836110,         0.938506208680101, 0.932502466241655, 0.926240223863451, 0.919723224389528, 0.912955362928245,         0.905940684524234, 0.898683381740745, 0.891187792153812, 0.883458395759753, 0.875499812297549,         0.867316798487695, 0.858914245189185, 0.850297174476322, 0.841470736637101, 0.832440207094966,         0.823210983255770, 0.813788581281829, 0.804178632795004, 0.794386881510759, 0.784419179805244,         0.774281485217411, 0.763979856888295, 0.753520451939554, 0.742909521793458, 0.732153408436510,         0.721258540628942, 0.710231430062342, 0.699078667467724, 0.687806918676346, 0.676422920635650,         0.664933477382692, 0.653345455977481, 0.641665782398637, 0.629901437403849, 0.618059452357584,         0.606146905028547, 0.594170915359415, 0.582138641211355, 0.570057274085885, 0.557934034826625,         0.545776169303514, 0.533590944082063, 0.521385642080249, 0.509167558215622, 0.496943995045254,         0.484722258401111, 0.472509653023471, 0.460313478195001, 0.448141023378083, 0.435999563858022,         0.423896356394721, 0.411838634885427, 0.399833606041146, 0.387888445079305, 0.376010291435238,         0.364206244495054, 0.352483359352449, 0.340848642591985, 0.329309048101366, 0.317871472915207,         0.306542753092784, 0.295329659632231, 0.284238894423618, 0.273277086243339, 0.262450786792195,         0.251766466779543, 0.241230512055860, 0.230849219796014, 0.220628794735546, 0.210575345462196,         0.200694880764894, 0.190993306042403, 0.181476419773753, 0.172149910052584, 0.163019351187458,         0.154090200370186, 0.145367794414147, 0.136857346564560, 0.128563943382607, 0.120492541705279,         0.112647965682748, 0.105034903895052, 0.0976579065498020, 0.0905213827625935,         0.0836295979217494, 0.0769866711389634, 0.0705965727873714, 0.0644631221285217,     0.0585899850296627, 0.0529806717727114, 0.0476385349562126, 0.0425667674915437,     0.0377684006945618, 0.0332463024738333, 0.0290031756165303, 0.0250415561730168,     0.0213638119410917, 0.0179721410507924, 0.0148685706506078, 0.0120549556958828,     0.00953297784014107, 0.00730414442998667, 0.00536978760418699, 0.00373106349747415,     0.00238895154954133, 0.00134425391964726, 0.000597595007177931, 0.000149421078453171    };

5）调用函数arm_cfft_f32();计算FFT结果

 1 arm_cfft_f32(&arm_cfft_sR_f32_len256 , fft_buf_float , 0 , 1); 2 //FFT正变换，输出结果需要进行位翻转 

其中第三个参数被置0，表示FFT算法是从时域向频域的正变换；第四个参数为1，表示需要对输出结果进行位翻转，以直接得到频率一次递增的频域数值。

6）将FFT结果转换为幅度谱

调用CMSIS-DSP库提供的函数arm_cmplx_mag_f32();计算信号的幅度谱（这将导致FFT结果中不同频率频谱的相位信息）。

 1 arm_cmplx_mag_f32(fft_buf_float , fft_out_float , 256); 2 //FFT函数输出也是复数，这里只取了幅度 

该函数的第一个参数指定了需要计算复数模的数组指针，第二个参数指定了计算结果存放的数组指针，第三个参数是需要计算复数模的数据个数。

7）将结果发送到PC进行显示

为更方便的观察和比较FFT结果，我使用以太网通过简单的UDP协议将数据上传到PC中显示。以下代码将结果存入缓冲区，其中float_un_buf是我事先定义的一个联合体，它由一个float数据和四个uchar数据构成。上述数据对象作为联合体成员，由于地址交叠，可以通过通过访问uchar数组将float数据转换为发送缓冲区所需的uchar类型。

 1 //把浮点结果放到无符号字节型的发送缓冲当中
 2 for(i = 0; i<256; i++)
 3 {
 4     float_un_buf.data_float = fft_out_float[i];
 5     Tx_Buffer[4*i] = float_un_buf.char_data[0];
 6     Tx_Buffer[4*i + 1] = float_un_buf.char_data[1];
 7     Tx_Buffer[4*i + 2] = float_un_buf.char_data[2];
 8     Tx_Buffer[4*i + 3] = float_un_buf.char_data[3];
 9 }
10 Write_SOCK_Data_Buffer(0, Tx_Buffer, 256*4);

四、实验结果

我使用STM32F4片上ADC（3.3V基准电压，采样率为1.4MSPS，DMA控制）对信号发生器产生的均值为1.65V，峰峰值为2V，频率202KHz的正弦信号进行了采样和FFT变换。然后将结果转换而成的幅度谱，通过UDP协议发送到PC上（接口芯片为W5500）。PC端程序开发平台为Visual Studio C# + NI Measurement studio。

1、观察频谱图（注意以下显示结果横轴数值不正确，由于是测试没有详细实现频域对应频率代码）

加汉宁窗并去除直流分量后显示结果如下图所示。

图4

 不去除直流分量后显示结果如下图5所示。可以明显观察到加到的直流分量。

图5

 不加汉宁窗结果如图6所示，可以观察到明显的能量泄露。

图6

 若将输入改为锯齿波，如图7所示。可以明显的观察到几次谐波分量。

图7

2、STM32F40x浮点FFT转换时间性能

通过J-Link仿真器和Keil-MDK环境，对我使用的STM32F401RC@84MHz打开FPU和关闭FPU条件下256点的浮点FFT转换时间进行了测试，结果为——打开FPU：251.4±5us；关闭FPU：4.9±100us

可见使用FPU后浮点FFT性能提高了约20倍。另外，读者如果想验证不实用FPU的算法时间，应将MDK工程中的库文件替代为arm_cortexM4l_math.lib，且options for target中的Floating Point Hardware选项（如图1所示）改为Not used。