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