本文介紹了使用STM32F4微控制器和ARM官方數字信号處理庫（CMSIS DSP LIB）進行數字信号處理的基本操作和程式設計開發方法，并以快速傅裡葉變換（FFT）為例，在MDK環境下，對STM32F4和CMSIS DSP LIB庫的性能進行了評估。
之前用過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 }; 256點窗函數
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。