：采樣、變換到特征提取)
簡介基于STM32的FFT頻譜分析與波形識別工程面向嵌入式開發(fā)者與信號處理學習者可完成50Hz~200Hz范圍內(nèi)基頻及3、5、7次諧波峰值測量并自動區(qū)分正弦、方波、鋸齒波、三角波。資源包含基礎版與可觸屏調(diào)采樣頻率的優(yōu)化版后者能針對不同信號實時調(diào)整采樣率以改善頻率分辨率適合需要兼顧算法實現(xiàn)與界面交互的實戰(zhàn)項目。壓縮包共395個文件以c/h源碼、o/d中間文件、crf編譯產(chǎn)物、uvprojx工程配置為主另含hex、map、axf、sct等燒錄與鏈接腳本整體約9.78MB。代碼基于正點原子精英版3.5寸TFTLCD燒寫即可運行已有3.3萬余人學習下載。通過閱讀工程可掌握STM32的ADC采樣、FFT變換、LCD波形繪制及觸屏交互設計同時附帶Keil工程備份與批處理腳本便于不同環(huán)境下的編譯與清理是學習數(shù)字信號處理與嵌入式顯示的完整范例。 STM32做FFT頻譜分析這個話題我前后調(diào)了一個多月才算真正跑順。你要問這個項目到底有什么用簡單說就一句話給單片機一雙“眼睛”讓它不僅能看波形長什么樣還能拆出這個波形背后由哪些頻率成分組成最后自動判斷出這是什么類型的信號。往實際應用上靠音頻頻譜顯示、信號發(fā)生器輸出自檢、電機振動監(jiān)測、電力諧波分析底層都是這套“采樣 FFT 特征識別”的組合。所以這不是一個只能用來交差畢業(yè)設計的玩具項目它更像一個隨時能塞進各種儀器設備里的核心模塊。這篇不打算從頭推導傅里葉變換的數(shù)學公式那玩意教材里講得夠多了我重點講從零把“STM32 FFT頻譜分析 波形識別”跑起來的過程中真正影響成敗的選型思路、參數(shù)計算、代碼實現(xiàn)和那些文檔里不會寫的坑。1. 項目拆解這套系統(tǒng)到底由哪些環(huán)節(jié)構(gòu)成1.1 三件事采樣、變換、識別看標題好像就兩個功能實際拆開是三件事采樣要保證數(shù)據(jù)均勻且準確FFT要把時域數(shù)據(jù)變成頻域數(shù)據(jù)波形識別則要綜合時域和頻域特征做判決。這三件事層層遞進前一步不扎實后面全是錯的。我見過不少人一上來就急著調(diào)FFT庫結(jié)果ADC采樣間隔亂跳FFT出來的頻譜跟撒了一把芝麻一樣根本沒法看。反過來也有人把FFT結(jié)果打印出來就以為完事了完全沒想怎么從一堆幅值數(shù)據(jù)里判定波形類型。其實整個項目的硬骨頭不在FFT本身而在“采樣鏈路的穩(wěn)定性”和“特征怎么提取”這兩塊。1.2 我最終采用的系統(tǒng)框圖整體架構(gòu)不復雜模擬信號經(jīng)過簡單的偏置和限幅電路送入STM32內(nèi)部ADC由定時器觸發(fā)轉(zhuǎn)換用DMA把結(jié)果搬到內(nèi)存攢夠一個FFT幀后做浮點變換再把幅值譜送去顯示同時跑一個識別函數(shù)判斷波形類型。軟件方面我用的是STM32CubeMX HAL庫 CMSIS-DSP這套組合是目前最主流的路徑網(wǎng)上資料多出問題也好查。標準庫新建工程那套流程我也走過但平心而論CubeMX生成初始化代碼的效率要高得多適合把精力集中在核心算法上。開發(fā)環(huán)境用Keil MDK或者VSCode EIDE都可以關(guān)鍵是把CMSIS-DSP的庫加進去這一步后面會細說。2. 選型對比芯片、FFT庫和采樣鏈路怎么選2.1 F103還是F407這不是糾結(jié)的問題如果你手里有F103的開發(fā)板說實話不推薦用在這個項目上。做FFT需要大量的乘加運算尤其是浮點運算。F103主頻72MHz沒有FPU軟件浮點做512點FFT勉強能跑1024點就要幾十毫秒加上顯示和識別刷新率會掉到每秒幾幀卡頓感非常明顯。F407就不一樣了168MHz主頻帶單精度硬件FPU和DSP指令跑1024點實數(shù)FFT在CMSIS-DSP庫的加持下實測大概是百微秒級別完全不影響實時刷新。還有一個容易被忽略的點F407的SRAM有192KBF103只有20KB。做4096點FFT時浮點數(shù)組緩沖區(qū)就要占用兩路N個float4096點大約是32KBF103不但性能跟不上內(nèi)存空間也直接卡死。如果預算有限退一步選帶FPU的G4系列或者F3系列也可以型號不同但思路完全一致。核心結(jié)論就一句這個項目盡量別用不帶FPU的老型號硬扛。2.2 FFT庫的選擇與配置FFT實現(xiàn)方式有三個層次自己寫基2算法、移植KissFFT等第三方庫、直接用CMSIS-DSP。我的建議是直接用CMSIS-DSPST官方在Cortex-M4上做了大量指令級優(yōu)化自己寫的代碼很難超過它。CMSIS-DSP里最常用的是arm_rfft_fast_f32專門處理實數(shù)序列因為ADC采出來的本來就是實數(shù)用這個函數(shù)比先構(gòu)造復數(shù)再調(diào)arm_cfft_f32要省一半運算。配置時要在工程里添加CMSIS-DSP源碼或預編譯庫用Keil的話在Manage Run-Time Environment里勾選DSP庫的Transform和Statistics功能就能自動引入。注意一個版本差異CMSIS-DSP有舊版arm_math.h和新版arm_rfft_fast_instance_f32結(jié)構(gòu)體定義位置有調(diào)整之分用新版時有些函數(shù)名一樣但頭文件路徑和宏定義不同編譯報錯多為“找不到arm_math.h”或者“identifier undefined”優(yōu)先檢查是否把DSP庫的Include路徑加全了。2.3 采樣鏈路為什么要用定時器觸發(fā)ADCFFT對采樣點的間隔均勻性非常敏感。如果靠main循環(huán)里的delay去采樣循環(huán)里每條指令的執(zhí)行時間不同會產(chǎn)生抖動這些抖動在頻域里表現(xiàn)為噪聲底座抬高小信號直接被淹沒。所以正確的做法是定時器產(chǎn)生更新事件把這個事件接到ADC的外部觸發(fā)引腳每次觸發(fā)啟動一次采樣轉(zhuǎn)換轉(zhuǎn)換完成由DMA自動把結(jié)果搬到內(nèi)存數(shù)組。整個過程中CPU不參與數(shù)據(jù)搬運只等DMA傳輸完成中斷來通知“一幀數(shù)據(jù)齊了”。這種硬件自動化的方式采樣間隔的抖動可以做到納秒級FFT效果會干凈很多。3. 參數(shù)計算與核心代碼實現(xiàn)3.1 采樣率、FFT點數(shù)和頻率分辨率怎么定玩FFT必須理解三個參數(shù)的關(guān)系采樣率Fs決定能分析的頻率范圍根據(jù)奈奎斯特定理最高分析頻率是Fs/2FFT點數(shù)N決定頻率分辨率相鄰兩根譜線的間隔是Fs/N采集一幀數(shù)據(jù)需要的時間就是N/Fs。舉個例子我用來做音頻頻段的實驗信號最高頻率按20kHz算Fs取48kHz就有富余。N取1024時頻率分辨率是48000/102446.875Hz這意味著50Hz和60Hz這種頻率根本分不開。如果你要觀察工頻或低頻信號就得把采樣率降下來比如Fs2048HzN1024分辨率就是2Hz能很清楚地分辨50Hz和60Hz。內(nèi)存占用也要提前算1024點FFT需要一個1024點的浮點輸入數(shù)組和一個1024點的輸出數(shù)組每個float占4字節(jié)總共8KB。再加上幅值數(shù)組和其他緩沖最好不要超過芯片SRAM的一半免得堆棧溢出。我用的是512點起步調(diào)試通后再升到1024點或者2048點循序漸進比較穩(wěn)妥。3.2 定時器觸發(fā)ADC DMA的配置要點CubeMX里的配置順序有講究。先把定時器時鐘源設為內(nèi)部時鐘然后算預分頻PSC和自動重裝ARR。以F407主頻168MHz、目標采樣率48kHz為例168000000 / 48000 3500那么可以取PSC1ARR1749因為168MHz / (11) / (17491) 48000Hz剛好精確命中。之后把定時器的Trigger Output設為Update Event這個TRGO信號就是給ADC用的。ADC側(cè)要關(guān)閉連續(xù)轉(zhuǎn)換模式外部觸發(fā)源選擇對應的Timer Trigger Out事件采樣周期我一般選較短的檔位比如15個ADC時鐘周期左右既能降低源阻抗要求又不會拖慢最大采樣率。DMA配置為Circular循環(huán)模式數(shù)據(jù)寬度設為Half Word因為12位ADC結(jié)果存在16位半字里。這里要強調(diào)一個檢查點ADC的觸發(fā)方式、DMA的外設地址和內(nèi)存地址不要搞反。外設地址要填ADC的數(shù)據(jù)寄存器地址內(nèi)存地址填你自己定義的數(shù)組首地址。很多工程用DMA采出來全是0多半是地址配錯或者寬度不一致。3.3 浮點FFT計算與幅值校準代碼做完采樣配置核心算法代碼其實很簡短。下面是我實際在用的FFT處理函數(shù)關(guān)鍵行都有注釋。#include arm_math.h #define FFT_SIZE 1024 static float32_t fft_input[FFT_SIZE]; static float32_t fft_output[FFT_SIZE]; static float32_t fft_mag[FFT_SIZE / 2]; static arm_rfft_fast_instance_f32 fft_inst; void FFT_Init(void) { arm_rfft_fast_init_f32(fft_inst, FFT_SIZE); } void FFT_Process(uint16_t *adc_buf) { uint16_t i; float32_t avg 0.0f; // ADC值轉(zhuǎn)電壓并去直流 for (i 0; i FFT_SIZE; i) { fft_input[i] (float32_t)adc_buf[i] * 3.3f / 4095.0f; avg fft_input[i]; } avg / FFT_SIZE; for (i 0; i FFT_SIZE; i) { fft_input[i] - avg; } // 0表示正變換 arm_rfft_fast_f32(fft_inst, fft_input, fft_output, 0); // 輸出排列: out[0]直流實部, out[1]Nyquist實部, // out[2]Re1, out[3]Im1, out[4]Re2, out[5]Im2 ... fft_mag[0] fft_output[0] / FFT_SIZE; // 直流 fft_mag[FFT_SIZE / 2 - 1] fft_output[1] / FFT_SIZE; // Nyquist for (i 1; i FFT_SIZE / 2; i) { float32_t re fft_output[2 * i]; float32_t im fft_output[2 * i 1]; // 單邊譜除直流和Nyquist外要乘2 fft_mag[i] sqrtf(re * re im * im) * 2.0f / FFT_SIZE; } }這里最容易踩的坑就是輸出數(shù)據(jù)排列格式。很多人把arm_rfft_fast_f32的輸出直接丟給arm_cmplx_mag_f32去求幅值結(jié)果完全不對。因為這個函數(shù)的輸出并不是普通復數(shù)交錯排列而是把直流和Nyquist單獨放在前兩個位置從第三個元素開始才是Re1, Im1, Re2, Im2這種交錯格式。要手動按上面的方式取數(shù)。去直流那一步也不能省。硬件上為了測交流信號經(jīng)常把信號偏置到1.65V這個直流分量如果不減掉會在0Hz處形成一個大尖峰不僅占掉顯示量程還可能掩蓋低頻信號。3.4 加窗處理頻譜泄漏的救星在采樣頻率和信號頻率不是整數(shù)倍關(guān)系時信號截斷會產(chǎn)生頻譜泄漏表現(xiàn)為真實譜線附近拖出一大片“裙邊”。最常用的處理是加漢寧窗。static float32_t window[FFT_SIZE]; void Window_Init(void) { for (uint16_t i 0; i FFT_SIZE; i) { window[i] 0.5f - 0.5f * cosf(2.0f * PI * i / (FFT_SIZE - 1)); } } // 在FFT_Process里、去直流之后執(zhí)行 for (i 0; i FFT_SIZE; i) { fft_input[i] * window[i]; }加窗是有代價的窗函數(shù)會攤寬主瓣降低頻率分辨率而且幅度會乘一個系數(shù)。漢寧窗的幅度恢復因子是2所以加窗后如果要做精確的幅值測量要在前面的乘2基礎上再乘2。不過如果只是看相對頻譜形態(tài)或者做波形識別乘不乘恢復因子影響不大因為所有譜線被同等縮放。如果信號頻率正好能對齊到bin中心比如Fs/N能整除信號頻率那不加窗也問題不大。但實際信號總有漂移所以我默認都加窗省心。4. 波形識別算法從特征到判決4.1 先做時域特征峰值系數(shù)區(qū)分大類FFT結(jié)果搞定了接下來是重頭戲——怎么識別波形類型。我的方法分兩步時域粗判頻域細判。時域最有效的特征之一是峰值系數(shù)也就是峰值除以有效值Crest Factor。理想情況下正弦波峰值系數(shù)是√2約1.414方波是1三角波是√3約1.732鋸齒波也是√3左右。只靠這個指標就能把方波和正弦波分得比較干凈。計算有效值需要在時域做不能直接用FFT幅值。我一般在ADC數(shù)據(jù)上算float CalcRMS(uint16_t *adc_buf, float dc_offset) { float sum 0.0f; for (uint16_t i 0; i FFT_SIZE; i) { float x (float32_t)adc_buf[i] * 3.3f / 4095.0f - dc_offset; sum x * x; } return sqrtf(sum / FFT_SIZE); } // 峰值就是整個數(shù)組里偏離直流最大的那個值 float peak GetPeak(adc_buf, dc_offset); float crest peak / CalcRMS(adc_buf, dc_offset);還要利用過零檢測測出基波頻率和FFT最大譜峰對應的頻率互相對照。兩者一致則說明FFT基頻找對了不一致說明可能采到了噪聲或者多頻信號識別置信度要打折。4.2 頻域諧波指紋奇次偶次和衰減規(guī)律時域粗判之后再用FFT結(jié)果看諧波結(jié)構(gòu)。不同波形的諧波指紋差異很明顯我把判斷依據(jù)整理成一張表波形諧波特征幅度衰減規(guī)律峰值系數(shù)正弦波基本只有基波諧波很少無約1.414方波只有奇次諧波1/k約1三角波只有奇次諧波1/k2約1.732鋸齒波奇次偶次都有1/k約1.732具體做法是找到最大譜峰作為基波記錄幅值h1然后在2倍、3倍、4倍……頻率處搜索對應譜峰記錄h2、h3、h4等。由此計算總諧波失真THD和奇偶次諧波能量比。float fundamental fft_mag[max_index]; float thd 0.0f; float even_energy 0.0f, odd_energy 0.0f; for (int k 2; k 8; k) { int idx (int)(max_index * k); if (idx FFT_SIZE / 2) break; float h fft_mag[idx]; thd h * h; if (k % 2 0) even_energy h; else odd_energy h; } thd sqrtf(thd) / fundamental;4.3 判決流程與容錯處理把時域和頻域信息合并我用一個簡單的決策樹判斷波形類型順序是先看峰值系數(shù)。如果接近1優(yōu)先判為方波即使頻譜上奇次諧波不全也要信時域的大方向。如果峰值系數(shù)接近1.414且THD小判為正弦波。如果峰值系數(shù)接近1.732再看頻域結(jié)構(gòu)存在明顯偶次諧波時判為鋸齒波只有奇次諧波時比較相鄰奇次諧波衰減速度接近1/k判方波接近1/k2判三角波。這里有個很現(xiàn)實的坑真實方波經(jīng)過低通濾波或帶寬限制后高階諧波會被削掉看起來會越來越像三角波。所以不要只依賴頻域衰減規(guī)律一定要結(jié)合時域峰值系數(shù)和波形包絡的斜率判斷。我實測中遇到過一次方波和三角波來回跳的情況后來增加了“連續(xù)3幀結(jié)果一致才更新識別結(jié)果”的濾波機制界面才穩(wěn)定下來。5. 調(diào)試實錄那些文檔里沒有的坑5.1 下載器連不上Error: No STM32 target found這個報錯幾乎是每個人都會被教育一次的經(jīng)典問題。完整報錯是“Error: No STM32 target found! If your product embeds Debug Authentication, please ...”。第一次遇到別慌按下面對照檢查SWD兩根線SWDIO和SWCLK有沒有接反GND是否共地。目標板電壓是否正常很多下載失敗其實是板子沒供電。按住板子復位鍵的同時點擊下載能連上就說明程序里把SWD引腳復用掉了。實在不行就把BOOT0拉高再上電用燒錄工具執(zhí)行全片擦除然后BOOT0拉回低電平。檢查一下Keil或CubeIDE里的Debug設置是不是選了Serial Wire有些新建工程默認沒勾選管腳全變普通GPIO自然連不上。另外如果電腦上“STM32 Virtual COM Port”設備顯示黃色感嘆號那是ST-Link的VCP驅(qū)動沒裝好重新裝一下STSW-LINK009驅(qū)動就解決了不會影響下載但串口打印數(shù)據(jù)時缺了它還真不行。5.2 頻譜圖全是噪點或大尖峰頻譜一團糟的排查順序很重要。先確認采樣率是否和預期一致。用示波器看定時器觸發(fā)引腳或者直接從一個已知頻率的信號發(fā)生器灌1kHz正弦波看FFT最大譜峰是不是正好落在1000Hz附近偏了就是時鐘樹或者預分頻算錯。再看FFT結(jié)果是不是只顯示前半部分了。很多人直接把全部N/2個點畫出來結(jié)果看到“鏡像頻譜”。記住實數(shù)FFT的結(jié)果關(guān)于Fs/2對稱只需顯示0到Fs/2這一段。噪聲底座過高的問題重點檢查有沒有用定時器硬件觸發(fā)采樣以及信號源阻抗是否太高導致ADC采樣保持電容充放電不充分。后者可以加一級運放電壓跟隨器解決。5.3 波形識別誤判的常見原因識別誤判這一塊我見過最多的情況是把帶限方波識別成三角波。原因是方波經(jīng)過前端抗混疊濾波器或者隔直電路后高次諧波衰減嚴重譜峰衰減規(guī)律接近1/k2。我的解決辦法是給判定邏輯加一個先驗條件先看時域的上升沿陡峭程度如果邊沿變化率明顯大于三角波的線性斜率即使諧波衰減快也優(yōu)先判定為方波。還有一個容易被忽略的點ADC采樣到的信號有直流偏置時波形識別會失敗因為正負半周不再對稱峰值系數(shù)完全失真。所以識別前一定要把直流偏置減掉而且要在偏置電壓穩(wěn)定的前提下進行。5.4 其他隱患延時卡死和系統(tǒng)穩(wěn)定性HAL_Delay卡死這個問題在FFT項目里也遇到過。情況通常是在DMA傳輸完成中斷里調(diào)用了HAL_Delay而HAL_Delay依賴SysTick中斷SysTick優(yōu)先級如果比DMA中斷低程序就死在等待標志位上。解決辦法很簡單中斷里不要調(diào)用HAL_Delay用狀態(tài)機或者變量打時間戳代替。另外如果刷新頻譜時發(fā)現(xiàn)畫面閃爍明顯多半是直接在顯示函數(shù)里做了大量浮點運算導致主循環(huán)變慢。合理做法是讓FFT和識別跑在主循環(huán)顯示環(huán)節(jié)用DMA送顯或者用雙緩沖交替渲染。做到這一步整個系統(tǒng)才算真正能拿出去用。這個項目做完后我最大的體會是FFT本身不難難的是把采樣鏈路做扎實、把特征閾值調(diào)到合適。調(diào)試波形識別的那個星期我把四種波形來回灌了上百遍慢慢才摸清不同容差對結(jié)果的影響。如果你也想做類似的東西建議先用電腦上的Python或Matlab把FFT和識別算法模擬一遍確認閾值靠譜后再往STM32上移植會省掉很多燒錄調(diào)試的時間。照著這套方案走你差不多兩三天就能看到自己的單片機上顯示出一條干凈的頻譜再花一兩天把波形識別調(diào)穩(wěn)這個項目就算真正落地了。本文還有配套的精品資源點擊獲取 STM32 FFT頻譜分析與波形識別實戰(zhàn)：采樣、變換到特征提取
簡介基于STM32的FFT頻譜分析與波形識別工程面向嵌入式開發(fā)者與信號處理學習者可完成50Hz~200Hz范圍內(nèi)基頻及3、5、7次諧波峰值測量并自動區(qū)分正弦、方波、鋸齒波、三角波。資源包含基礎版與可觸屏調(diào)采樣頻率的優(yōu)化版后者能針對不同信號實時調(diào)整采樣率以改善頻率分辨率適合需要兼顧算法實現(xiàn)與界面交互的實戰(zhàn)項目。壓縮包共395個文件以c/h源碼、o/d中間文件、crf編譯產(chǎn)物、uvprojx工程配置為主另含hex、map、axf、sct等燒錄與鏈接腳本整體約9.78MB。代碼基于正點原子精英版3.5寸TFTLCD燒寫即可運行已有3.3萬余人學習下載。通過閱讀工程可掌握STM32的ADC采樣、FFT變換、LCD波形繪制及觸屏交互設計同時附帶Keil工程備份與批處理腳本便于不同環(huán)境下的編譯與清理是學習數(shù)字信號處理與嵌入式顯示的完整范例。 STM32做FFT頻譜分析這個話題我前后調(diào)了一個多月才算真正跑順。你要問這個項目到底有什么用簡單說就一句話給單片機一雙“眼睛”讓它不僅能看波形長什么樣還能拆出這個波形背后由哪些頻率成分組成最后自動判斷出這是什么類型的信號。往實際應用上靠音頻頻譜顯示、信號發(fā)生器輸出自檢、電機振動監(jiān)測、電力諧波分析底層都是這套“采樣 FFT 特征識別”的組合。所以這不是一個只能用來交差畢業(yè)設計的玩具項目它更像一個隨時能塞進各種儀器設備里的核心模塊。這篇不打算從頭推導傅里葉變換的數(shù)學公式那玩意教材里講得夠多了我重點講從零把“STM32 FFT頻譜分析 波形識別”跑起來的過程中真正影響成敗的選型思路、參數(shù)計算、代碼實現(xiàn)和那些文檔里不會寫的坑。1. 項目拆解這套系統(tǒng)到底由哪些環(huán)節(jié)構(gòu)成1.1 三件事采樣、變換、識別看標題好像就兩個功能實際拆開是三件事采樣要保證數(shù)據(jù)均勻且準確FFT要把時域數(shù)據(jù)變成頻域數(shù)據(jù)波形識別則要綜合時域和頻域特征做判決。這三件事層層遞進前一步不扎實后面全是錯的。我見過不少人一上來就急著調(diào)FFT庫結(jié)果ADC采樣間隔亂跳FFT出來的頻譜跟撒了一把芝麻一樣根本沒法看。反過來也有人把FFT結(jié)果打印出來就以為完事了完全沒想怎么從一堆幅值數(shù)據(jù)里判定波形類型。其實整個項目的硬骨頭不在FFT本身而在“采樣鏈路的穩(wěn)定性”和“特征怎么提取”這兩塊。1.2 我最終采用的系統(tǒng)框圖整體架構(gòu)不復雜模擬信號經(jīng)過簡單的偏置和限幅電路送入STM32內(nèi)部ADC由定時器觸發(fā)轉(zhuǎn)換用DMA把結(jié)果搬到內(nèi)存攢夠一個FFT幀后做浮點變換再把幅值譜送去顯示同時跑一個識別函數(shù)判斷波形類型。軟件方面我用的是STM32CubeMX HAL庫 CMSIS-DSP這套組合是目前最主流的路徑網(wǎng)上資料多出問題也好查。標準庫新建工程那套流程我也走過但平心而論CubeMX生成初始化代碼的效率要高得多適合把精力集中在核心算法上。開發(fā)環(huán)境用Keil MDK或者VSCode EIDE都可以關(guān)鍵是把CMSIS-DSP的庫加進去這一步后面會細說。2. 選型對比芯片、FFT庫和采樣鏈路怎么選2.1 F103還是F407這不是糾結(jié)的問題如果你手里有F103的開發(fā)板說實話不推薦用在這個項目上。做FFT需要大量的乘加運算尤其是浮點運算。F103主頻72MHz沒有FPU軟件浮點做512點FFT勉強能跑1024點就要幾十毫秒加上顯示和識別刷新率會掉到每秒幾幀卡頓感非常明顯。F407就不一樣了168MHz主頻帶單精度硬件FPU和DSP指令跑1024點實數(shù)FFT在CMSIS-DSP庫的加持下實測大概是百微秒級別完全不影響實時刷新。還有一個容易被忽略的點F407的SRAM有192KBF103只有20KB。做4096點FFT時浮點數(shù)組緩沖區(qū)就要占用兩路N個float4096點大約是32KBF103不但性能跟不上內(nèi)存空間也直接卡死。如果預算有限退一步選帶FPU的G4系列或者F3系列也可以型號不同但思路完全一致。核心結(jié)論就一句這個項目盡量別用不帶FPU的老型號硬扛。2.2 FFT庫的選擇與配置FFT實現(xiàn)方式有三個層次自己寫基2算法、移植KissFFT等第三方庫、直接用CMSIS-DSP。我的建議是直接用CMSIS-DSPST官方在Cortex-M4上做了大量指令級優(yōu)化自己寫的代碼很難超過它。CMSIS-DSP里最常用的是arm_rfft_fast_f32專門處理實數(shù)序列因為ADC采出來的本來就是實數(shù)用這個函數(shù)比先構(gòu)造復數(shù)再調(diào)arm_cfft_f32要省一半運算。配置時要在工程里添加CMSIS-DSP源碼或預編譯庫用Keil的話在Manage Run-Time Environment里勾選DSP庫的Transform和Statistics功能就能自動引入。注意一個版本差異CMSIS-DSP有舊版arm_math.h和新版arm_rfft_fast_instance_f32結(jié)構(gòu)體定義位置有調(diào)整之分用新版時有些函數(shù)名一樣但頭文件路徑和宏定義不同編譯報錯多為“找不到arm_math.h”或者“identifier undefined”優(yōu)先檢查是否把DSP庫的Include路徑加全了。2.3 采樣鏈路為什么要用定時器觸發(fā)ADCFFT對采樣點的間隔均勻性非常敏感。如果靠main循環(huán)里的delay去采樣循環(huán)里每條指令的執(zhí)行時間不同會產(chǎn)生抖動這些抖動在頻域里表現(xiàn)為噪聲底座抬高小信號直接被淹沒。所以正確的做法是定時器產(chǎn)生更新事件把這個事件接到ADC的外部觸發(fā)引腳每次觸發(fā)啟動一次采樣轉(zhuǎn)換轉(zhuǎn)換完成由DMA自動把結(jié)果搬到內(nèi)存數(shù)組。整個過程中CPU不參與數(shù)據(jù)搬運只等DMA傳輸完成中斷來通知“一幀數(shù)據(jù)齊了”。這種硬件自動化的方式采樣間隔的抖動可以做到納秒級FFT效果會干凈很多。3. 參數(shù)計算與核心代碼實現(xiàn)3.1 采樣率、FFT點數(shù)和頻率分辨率怎么定玩FFT必須理解三個參數(shù)的關(guān)系采樣率Fs決定能分析的頻率范圍根據(jù)奈奎斯特定理最高分析頻率是Fs/2FFT點數(shù)N決定頻率分辨率相鄰兩根譜線的間隔是Fs/N采集一幀數(shù)據(jù)需要的時間就是N/Fs。舉個例子我用來做音頻頻段的實驗信號最高頻率按20kHz算Fs取48kHz就有富余。N取1024時頻率分辨率是48000/102446.875Hz這意味著50Hz和60Hz這種頻率根本分不開。如果你要觀察工頻或低頻信號就得把采樣率降下來比如Fs2048HzN1024分辨率就是2Hz能很清楚地分辨50Hz和60Hz。內(nèi)存占用也要提前算1024點FFT需要一個1024點的浮點輸入數(shù)組和一個1024點的輸出數(shù)組每個float占4字節(jié)總共8KB。再加上幅值數(shù)組和其他緩沖最好不要超過芯片SRAM的一半免得堆棧溢出。我用的是512點起步調(diào)試通后再升到1024點或者2048點循序漸進比較穩(wěn)妥。3.2 定時器觸發(fā)ADC DMA的配置要點CubeMX里的配置順序有講究。先把定時器時鐘源設為內(nèi)部時鐘然后算預分頻PSC和自動重裝ARR。以F407主頻168MHz、目標采樣率48kHz為例168000000 / 48000 3500那么可以取PSC1ARR1749因為168MHz / (11) / (17491) 48000Hz剛好精確命中。之后把定時器的Trigger Output設為Update Event這個TRGO信號就是給ADC用的。ADC側(cè)要關(guān)閉連續(xù)轉(zhuǎn)換模式外部觸發(fā)源選擇對應的Timer Trigger Out事件采樣周期我一般選較短的檔位比如15個ADC時鐘周期左右既能降低源阻抗要求又不會拖慢最大采樣率。DMA配置為Circular循環(huán)模式數(shù)據(jù)寬度設為Half Word因為12位ADC結(jié)果存在16位半字里。這里要強調(diào)一個檢查點ADC的觸發(fā)方式、DMA的外設地址和內(nèi)存地址不要搞反。外設地址要填ADC的數(shù)據(jù)寄存器地址內(nèi)存地址填你自己定義的數(shù)組首地址。很多工程用DMA采出來全是0多半是地址配錯或者寬度不一致。3.3 浮點FFT計算與幅值校準代碼做完采樣配置核心算法代碼其實很簡短。下面是我實際在用的FFT處理函數(shù)關(guān)鍵行都有注釋。#include arm_math.h #define FFT_SIZE 1024 static float32_t fft_input[FFT_SIZE]; static float32_t fft_output[FFT_SIZE]; static float32_t fft_mag[FFT_SIZE / 2]; static arm_rfft_fast_instance_f32 fft_inst; void FFT_Init(void) { arm_rfft_fast_init_f32(fft_inst, FFT_SIZE); } void FFT_Process(uint16_t *adc_buf) { uint16_t i; float32_t avg 0.0f; // ADC值轉(zhuǎn)電壓并去直流 for (i 0; i FFT_SIZE; i) { fft_input[i] (float32_t)adc_buf[i] * 3.3f / 4095.0f; avg fft_input[i]; } avg / FFT_SIZE; for (i 0; i FFT_SIZE; i) { fft_input[i] - avg; } // 0表示正變換 arm_rfft_fast_f32(fft_inst, fft_input, fft_output, 0); // 輸出排列: out[0]直流實部, out[1]Nyquist實部, // out[2]Re1, out[3]Im1, out[4]Re2, out[5]Im2 ... fft_mag[0] fft_output[0] / FFT_SIZE; // 直流 fft_mag[FFT_SIZE / 2 - 1] fft_output[1] / FFT_SIZE; // Nyquist for (i 1; i FFT_SIZE / 2; i) { float32_t re fft_output[2 * i]; float32_t im fft_output[2 * i 1]; // 單邊譜除直流和Nyquist外要乘2 fft_mag[i] sqrtf(re * re im * im) * 2.0f / FFT_SIZE; } }這里最容易踩的坑就是輸出數(shù)據(jù)排列格式。很多人把arm_rfft_fast_f32的輸出直接丟給arm_cmplx_mag_f32去求幅值結(jié)果完全不對。因為這個函數(shù)的輸出并不是普通復數(shù)交錯排列而是把直流和Nyquist單獨放在前兩個位置從第三個元素開始才是Re1, Im1, Re2, Im2這種交錯格式。要手動按上面的方式取數(shù)。去直流那一步也不能省。硬件上為了測交流信號經(jīng)常把信號偏置到1.65V這個直流分量如果不減掉會在0Hz處形成一個大尖峰不僅占掉顯示量程還可能掩蓋低頻信號。3.4 加窗處理頻譜泄漏的救星在采樣頻率和信號頻率不是整數(shù)倍關(guān)系時信號截斷會產(chǎn)生頻譜泄漏表現(xiàn)為真實譜線附近拖出一大片“裙邊”。最常用的處理是加漢寧窗。static float32_t window[FFT_SIZE]; void Window_Init(void) { for (uint16_t i 0; i FFT_SIZE; i) { window[i] 0.5f - 0.5f * cosf(2.0f * PI * i / (FFT_SIZE - 1)); } } // 在FFT_Process里、去直流之后執(zhí)行 for (i 0; i FFT_SIZE; i) { fft_input[i] * window[i]; }加窗是有代價的窗函數(shù)會攤寬主瓣降低頻率分辨率而且幅度會乘一個系數(shù)。漢寧窗的幅度恢復因子是2所以加窗后如果要做精確的幅值測量要在前面的乘2基礎上再乘2。不過如果只是看相對頻譜形態(tài)或者做波形識別乘不乘恢復因子影響不大因為所有譜線被同等縮放。如果信號頻率正好能對齊到bin中心比如Fs/N能整除信號頻率那不加窗也問題不大。但實際信號總有漂移所以我默認都加窗省心。4. 波形識別算法從特征到判決4.1 先做時域特征峰值系數(shù)區(qū)分大類FFT結(jié)果搞定了接下來是重頭戲——怎么識別波形類型。我的方法分兩步時域粗判頻域細判。時域最有效的特征之一是峰值系數(shù)也就是峰值除以有效值Crest Factor。理想情況下正弦波峰值系數(shù)是√2約1.414方波是1三角波是√3約1.732鋸齒波也是√3左右。只靠這個指標就能把方波和正弦波分得比較干凈。計算有效值需要在時域做不能直接用FFT幅值。我一般在ADC數(shù)據(jù)上算float CalcRMS(uint16_t *adc_buf, float dc_offset) { float sum 0.0f; for (uint16_t i 0; i FFT_SIZE; i) { float x (float32_t)adc_buf[i] * 3.3f / 4095.0f - dc_offset; sum x * x; } return sqrtf(sum / FFT_SIZE); } // 峰值就是整個數(shù)組里偏離直流最大的那個值 float peak GetPeak(adc_buf, dc_offset); float crest peak / CalcRMS(adc_buf, dc_offset);還要利用過零檢測測出基波頻率和FFT最大譜峰對應的頻率互相對照。兩者一致則說明FFT基頻找對了不一致說明可能采到了噪聲或者多頻信號識別置信度要打折。4.2 頻域諧波指紋奇次偶次和衰減規(guī)律時域粗判之后再用FFT結(jié)果看諧波結(jié)構(gòu)。不同波形的諧波指紋差異很明顯我把判斷依據(jù)整理成一張表波形諧波特征幅度衰減規(guī)律峰值系數(shù)正弦波基本只有基波諧波很少無約1.414方波只有奇次諧波1/k約1三角波只有奇次諧波1/k2約1.732鋸齒波奇次偶次都有1/k約1.732具體做法是找到最大譜峰作為基波記錄幅值h1然后在2倍、3倍、4倍……頻率處搜索對應譜峰記錄h2、h3、h4等。由此計算總諧波失真THD和奇偶次諧波能量比。float fundamental fft_mag[max_index]; float thd 0.0f; float even_energy 0.0f, odd_energy 0.0f; for (int k 2; k 8; k) { int idx (int)(max_index * k); if (idx FFT_SIZE / 2) break; float h fft_mag[idx]; thd h * h; if (k % 2 0) even_energy h; else odd_energy h; } thd sqrtf(thd) / fundamental;4.3 判決流程與容錯處理把時域和頻域信息合并我用一個簡單的決策樹判斷波形類型順序是先看峰值系數(shù)。如果接近1優(yōu)先判為方波即使頻譜上奇次諧波不全也要信時域的大方向。如果峰值系數(shù)接近1.414且THD小判為正弦波。如果峰值系數(shù)接近1.732再看頻域結(jié)構(gòu)存在明顯偶次諧波時判為鋸齒波只有奇次諧波時比較相鄰奇次諧波衰減速度接近1/k判方波接近1/k2判三角波。這里有個很現(xiàn)實的坑真實方波經(jīng)過低通濾波或帶寬限制后高階諧波會被削掉看起來會越來越像三角波。所以不要只依賴頻域衰減規(guī)律一定要結(jié)合時域峰值系數(shù)和波形包絡的斜率判斷。我實測中遇到過一次方波和三角波來回跳的情況后來增加了“連續(xù)3幀結(jié)果一致才更新識別結(jié)果”的濾波機制界面才穩(wěn)定下來。5. 調(diào)試實錄那些文檔里沒有的坑5.1 下載器連不上Error: No STM32 target found這個報錯幾乎是每個人都會被教育一次的經(jīng)典問題。完整報錯是“Error: No STM32 target found! If your product embeds Debug Authentication, please ...”。第一次遇到別慌按下面對照檢查SWD兩根線SWDIO和SWCLK有沒有接反GND是否共地。目標板電壓是否正常很多下載失敗其實是板子沒供電。按住板子復位鍵的同時點擊下載能連上就說明程序里把SWD引腳復用掉了。實在不行就把BOOT0拉高再上電用燒錄工具執(zhí)行全片擦除然后BOOT0拉回低電平。檢查一下Keil或CubeIDE里的Debug設置是不是選了Serial Wire有些新建工程默認沒勾選管腳全變普通GPIO自然連不上。另外如果電腦上“STM32 Virtual COM Port”設備顯示黃色感嘆號那是ST-Link的VCP驅(qū)動沒裝好重新裝一下STSW-LINK009驅(qū)動就解決了不會影響下載但串口打印數(shù)據(jù)時缺了它還真不行。5.2 頻譜圖全是噪點或大尖峰頻譜一團糟的排查順序很重要。先確認采樣率是否和預期一致。用示波器看定時器觸發(fā)引腳或者直接從一個已知頻率的信號發(fā)生器灌1kHz正弦波看FFT最大譜峰是不是正好落在1000Hz附近偏了就是時鐘樹或者預分頻算錯。再看FFT結(jié)果是不是只顯示前半部分了。很多人直接把全部N/2個點畫出來結(jié)果看到“鏡像頻譜”。記住實數(shù)FFT的結(jié)果關(guān)于Fs/2對稱只需顯示0到Fs/2這一段。噪聲底座過高的問題重點檢查有沒有用定時器硬件觸發(fā)采樣以及信號源阻抗是否太高導致ADC采樣保持電容充放電不充分。后者可以加一級運放電壓跟隨器解決。5.3 波形識別誤判的常見原因識別誤判這一塊我見過最多的情況是把帶限方波識別成三角波。原因是方波經(jīng)過前端抗混疊濾波器或者隔直電路后高次諧波衰減嚴重譜峰衰減規(guī)律接近1/k2。我的解決辦法是給判定邏輯加一個先驗條件先看時域的上升沿陡峭程度如果邊沿變化率明顯大于三角波的線性斜率即使諧波衰減快也優(yōu)先判定為方波。還有一個容易被忽略的點ADC采樣到的信號有直流偏置時波形識別會失敗因為正負半周不再對稱峰值系數(shù)完全失真。所以識別前一定要把直流偏置減掉而且要在偏置電壓穩(wěn)定的前提下進行。5.4 其他隱患延時卡死和系統(tǒng)穩(wěn)定性HAL_Delay卡死這個問題在FFT項目里也遇到過。情況通常是在DMA傳輸完成中斷里調(diào)用了HAL_Delay而HAL_Delay依賴SysTick中斷SysTick優(yōu)先級如果比DMA中斷低程序就死在等待標志位上。解決辦法很簡單中斷里不要調(diào)用HAL_Delay用狀態(tài)機或者變量打時間戳代替。另外如果刷新頻譜時發(fā)現(xiàn)畫面閃爍明顯多半是直接在顯示函數(shù)里做了大量浮點運算導致主循環(huán)變慢。合理做法是讓FFT和識別跑在主循環(huán)顯示環(huán)節(jié)用DMA送顯或者用雙緩沖交替渲染。做到這一步整個系統(tǒng)才算真正能拿出去用。這個項目做完后我最大的體會是FFT本身不難難的是把采樣鏈路做扎實、把特征閾值調(diào)到合適。調(diào)試波形識別的那個星期我把四種波形來回灌了上百遍慢慢才摸清不同容差對結(jié)果的影響。如果你也想做類似的東西建議先用電腦上的Python或Matlab把FFT和識別算法模擬一遍確認閾值靠譜后再往STM32上移植會省掉很多燒錄調(diào)試的時間。照著這套方案走你差不多兩三天就能看到自己的單片機上顯示出一條干凈的頻譜再花一兩天把波形識別調(diào)穩(wěn)這個項目就算真正落地了。本文還有配套的精品資源點擊獲取