IIR濾波器這名字初聽起來像教科書里的概念但我敢說幾乎所有接觸過數(shù)字信號處理的人最終都得跟它打交道。你手上那塊STM32讀進(jìn)來的電壓、麥克風(fēng)收到的人聲、傳感器采集的振動波形想要去掉噪聲或者提取特征頻率簡單粗暴的均值濾波不夠用、FFT又不實(shí)時的時候IIR濾波器就是最直接的那把刀。這篇內(nèi)容我從工程落地的角度把IIR濾波器從設(shè)計(jì)原理、系數(shù)提取、代碼實(shí)現(xiàn)到SOS矩陣整理成一條完整鏈路重點(diǎn)覆蓋STM32這類MCU上用“直接I型”和“SOS矩陣”實(shí)現(xiàn)濾波器時會遇到的坑和心得適合剛接觸數(shù)字濾波的嵌入式工程師、音頻開發(fā)者以及做信號采集相關(guān)項(xiàng)目但還沒系統(tǒng)梳理過濾波器設(shè)計(jì)的朋友。1. 為什么選擇IIR與FIR的本質(zhì)差異1.1 從傳遞函數(shù)角度看兩者的不同IIR濾波器的全稱是無限脈沖響應(yīng)濾波器它的核心特征可以用z域傳遞函數(shù)來表示。一個N階IIR系統(tǒng)的傳遞函數(shù)長這樣H(z) (b0 b1·z?1 b2·z?2 ... bM·z??) / (1 a1·z?1 a2·z?2 ... aN·z??)分母里那串系數(shù)a就構(gòu)成了反饋回路這是IIR和FIR最根本的區(qū)別。FIR濾波器的傳遞函數(shù)只有分子沒有分母也就沒有反饋它的脈沖響應(yīng)在有限個采樣點(diǎn)后就歸零了所以叫“有限脈沖響應(yīng)”。而IIR因?yàn)橛蟹帜?、有反饋理論上一個脈沖輸入會在輸出端留下一串無限延長的響應(yīng)尾巴這才是“無限脈沖響應(yīng)”這個名字的由來。別看這只是數(shù)學(xué)形式上的差別它帶來的工程影響非常大。分母多了一項(xiàng)意味著IIR可以用低得多的階數(shù)實(shí)現(xiàn)同樣陡峭的過渡帶。舉個我實(shí)際對比過的例子一個采樣率48kHz、截止頻率10kHz的低通濾波器如果用FIR要達(dá)到60dB阻帶衰減和較窄過渡帶至少需要一百多階每輸入一個采樣點(diǎn)就要做一百多次乘加運(yùn)算而同樣規(guī)格的IIR我用一個4階的橢圓濾波器就搞定了單拍計(jì)算量只有FIR的零頭。這個差距在音頻、實(shí)時控制這類對延時和算力敏感的場景里基本就是能不能跑得動的區(qū)別。1.2 IIR真實(shí)性能與代價當(dāng)然天下沒有免費(fèi)的午餐。IIR用低階換來的代價是“相位非線性”和“穩(wěn)定性風(fēng)險(xiǎn)”。相位非線性這一點(diǎn)很多人一開始不在意但真正用起來就會發(fā)現(xiàn)問題。你在示波器上看波形可能覺得沒差別可如果處理的是心電信號、振動信號或者任何需要保持波形原始形態(tài)的場景IIR會把不同頻率成分的延時變得不一樣導(dǎo)致輸出波形“變形”——這不是幅度上的變形而是時間軸上的錯位。FIR因?yàn)橄禂?shù)對稱設(shè)計(jì)天然能做到線性相位所有頻率延時完全一致波形不會走樣。所以選型的時候我一般先問自己一個問題這個應(yīng)用對相位敏感嗎如果只是去除電源紋波、平滑傳感器噪聲、做音頻音調(diào)調(diào)整IIR完全合適如果要做帶通濾波后的時域特征分析、波達(dá)時間估計(jì)那我建議要么用FIR要么接受IIR非線性相位帶來的誤差要么用零相位雙向?yàn)V波離線場景彌補(bǔ)。穩(wěn)定性風(fēng)險(xiǎn)則是另一個大坑。既然IIR有反饋分母那就天然存在“會不會發(fā)散”的問題。從z域看系統(tǒng)穩(wěn)定的條件是所有極點(diǎn)都落在z平面單位圓內(nèi)。問題在于當(dāng)我們把設(shè)計(jì)好的濾波器搬進(jìn)16位或32位定點(diǎn)芯片、用有限精度的浮點(diǎn)運(yùn)算實(shí)現(xiàn)時系數(shù)的量化誤差可能導(dǎo)致原本在單位圓內(nèi)、離邊界很遠(yuǎn)的極點(diǎn)漂移到圓外系統(tǒng)就變成了一個“振蕩器”。這也是后面要重點(diǎn)聊SOS矩陣的直接原因——把高階系統(tǒng)拆成多個二階子系統(tǒng)級聯(lián)能顯著降低系數(shù)量化誤差帶來的極點(diǎn)漂移風(fēng)險(xiǎn)讓濾波器在真實(shí)硬件上更穩(wěn)定。1.3 選型判斷什么時候該用IIR說了那么多我可以把這些年做工程選型的經(jīng)驗(yàn)整理成一張表幫助快速決策對比項(xiàng)IIRFIR相同指標(biāo)所需階數(shù)低通常4~10階高動輒幾十上百階單點(diǎn)計(jì)算量小大相位特征非線性可做到線性相位穩(wěn)定性有極點(diǎn)需關(guān)注全零點(diǎn)無條件穩(wěn)定適合場景實(shí)時控制、音頻均衡、噪聲抑制相位敏感測量、多速率處理、心理聲學(xué)典型工具巴特沃斯、切比雪夫、橢圓窗函數(shù)法、等效波法這張表不是說要大家死記硬背我更想強(qiáng)調(diào)的是“實(shí)時性”這個維度。在STM32這類MCU上一個中斷里可能同時要做ADC采集、濾波、PID控制、通信處理留給濾波器的預(yù)算往往只有幾個微秒到幾十個微秒。FIR要是上千階光這部分的乘加運(yùn)算就能把CPU吃滿。IIR用個小幾十次的乘加就完成了同樣指標(biāo)這種差距是實(shí)際測出來的不是紙面數(shù)據(jù)能體現(xiàn)的。結(jié)論很簡單高性能計(jì)算平臺、相位敏感場景選FIR嵌入式實(shí)時場景、噪聲抑制與平滑場景IIR是性價比極高的選擇。而且IIR配合SOS結(jié)構(gòu)穩(wěn)定性和可維護(hù)性都能得到保障這也是它至今在工業(yè)控制、音頻、傳感器處理領(lǐng)域占據(jù)重要位置的原因。2. 設(shè)計(jì)IIR濾波器的完整流程與工具實(shí)操2.1 設(shè)計(jì)規(guī)格的確定很多人拿到一個需求就急著打開工具箱敲代碼我建議先花十分鐘把設(shè)計(jì)規(guī)格寫清楚這一步省了后面返工十倍的功夫。所謂設(shè)計(jì)規(guī)格核心就是五個參數(shù)采樣率fs、通帶邊緣頻率、阻帶邊緣頻率、通帶最大紋波、阻帶最小衰減。采樣率是這一切的基礎(chǔ)奈奎斯特頻率fs/2決定了你能處理的最高頻率分量。在數(shù)字濾波器設(shè)計(jì)里所有頻率參數(shù)最終都要?dú)w一化到這個奈奎斯特頻率上。比如采樣率1000Hz、想保留100Hz以內(nèi)的信號那歸一化截止頻率就是100 / (1000/2) 0.2。這個歸一化過程用工具的時候也會遇到只是大多數(shù)時候工具幫你做了而已。通帶紋波和阻帶衰減這兩個參數(shù)決定了濾波器的“質(zhì)量”。通帶紋波1dB意味著信號通過后在幅度上會有大約±5.6%的波動這在音頻應(yīng)用里通常感知不到但在測量儀器上就可能影響精度。阻帶衰減40dB意味著阻帶信號能削弱到原來的1%這對大多數(shù)工業(yè)應(yīng)用已經(jīng)夠了想做到60dB、80dB也不是不行只是濾波器階數(shù)或者設(shè)計(jì)算法的復(fù)雜度會上升。設(shè)計(jì)規(guī)格不是越嚴(yán)越好跟機(jī)械加工的公差一樣越嚴(yán)成本越高這里的“成本”就是階數(shù)和運(yùn)算量。2.2 用Python工具箱完成系數(shù)設(shè)計(jì)我自己最常用的設(shè)計(jì)工具是Python的SciPy庫因?yàn)樗饶茉赑C上快速驗(yàn)證又能把設(shè)計(jì)出的系數(shù)直接移植到C代碼里。下面是一個實(shí)際可運(yùn)行的低通濾波器設(shè)計(jì)示例import numpy as np from scipy.signal import ellip, butter, cheby1, sosfreqz fs 1000.0 # 采樣率 1kHz fc 100.0 # 通帶截止頻率 100Hz order 4 # 濾波器階數(shù) rp 1.0 # 通帶紋波 1dB rs 40.0 # 阻帶衰減 40dB # 用橢圓濾波器設(shè)計(jì)輸出SOS矩陣格式 sos ellip(order, rp, rs, fc/(fs/2), outputsos) print(sos) # 打印頻率響應(yīng)做驗(yàn)證 w, h sosfreqz(sos, worN2048, fsfs)這里用橢圓濾波器是因?yàn)樗谕瑯拥碾A數(shù)下過渡帶最窄適合在低階條件下實(shí)現(xiàn)盡可能陡的衰減。如果你對通帶內(nèi)的紋波非常敏感、不想要橢圓濾波器的等紋波紋波特征可以把ellip換成butter巴特沃斯它的通帶響應(yīng)是最平坦的但過渡帶會寬一些。切比雪夫I型在通帶內(nèi)有等紋波、阻帶單調(diào)切比雪夫II型則相反。選型邏輯很簡單追求最陡過渡帶選橢圓追求通帶平坦選巴特沃斯兩種極端之間選切比雪夫。2.3 設(shè)計(jì)結(jié)果校驗(yàn)?zāi)玫较禂?shù)之后別急著抄進(jìn)代碼。我習(xí)慣先打印頻率響應(yīng)曲線做人眼驗(yàn)證——確認(rèn)通帶內(nèi)增益接近0dB、截止頻率位置正確、阻帶衰減達(dá)到設(shè)計(jì)目標(biāo)。有一回我設(shè)計(jì)一個帶通濾波器系數(shù)算出來一切正常結(jié)果畫幅頻圖發(fā)現(xiàn)中心頻率偏了5Hz仔細(xì)排查發(fā)現(xiàn)是設(shè)計(jì)時歸一化頻率的參考搞錯了。這個錯誤要是直接燒進(jìn)板子可能要調(diào)大半天才能發(fā)現(xiàn)。SciPy的sosfreqz函數(shù)可以直接對SOS格式的系數(shù)計(jì)算頻率響應(yīng)也可以先用tf2sos把普通的分子分母系數(shù)轉(zhuǎn)成SOS。如果做的是離線處理還可以直接用sosfilt函數(shù)一次性跑完整段數(shù)據(jù)來驗(yàn)證濾波效果把輸入信號和輸出信號放在一起對比看看時域波形是不是符合預(yù)期。這個“PC端設(shè)計(jì)→PC端驗(yàn)證→移植到MCU”的流程我到現(xiàn)在還會用因?yàn)镸CU上調(diào)試濾波的代價比PC上高太多。3. 直接I型濾波器在MCU/STM32上的代碼實(shí)現(xiàn)3.1 直接I型結(jié)構(gòu)的數(shù)學(xué)推導(dǎo)拿到系數(shù)之后怎么把它變成能在MCU上跑的代碼呢有一個概念要先講清楚差分方程。IIR濾波器的傳遞函數(shù)可以等價轉(zhuǎn)換成時域里的差分方程代碼就是逐樣本地執(zhí)行這個差分方程。假設(shè)我們有這樣一個2階傳遞函數(shù)H(z) (b0 b1·z?1 b2·z?2) / (1 a1·z?1 a2·z?2)對應(yīng)的差分方程為y[n] b0·x[n] b1·x[n-1] b2·x[n-2] - a1·y[n-1] - a2·y[n-2]注意分母里的a1、a2在方程里變成了減法而且之前提到的歸一化要求a0 1如果設(shè)計(jì)工具給出的系數(shù)里a0不是1需要先把所有系數(shù)除以a0。x[n-1]、x[n-2]是前兩個輸入樣本y[n-1]、y[n-2]是前兩個輸出樣本。實(shí)現(xiàn)這個方程的時候只需要維護(hù)四個歷史變量就夠了。這種結(jié)構(gòu)之所以叫“直接I型”是因?yàn)樗苯诱罩罘址匠讨痦?xiàng)實(shí)現(xiàn)輸入歷史先走分子輸出歷史再走分母。另外還有“直接II型”和“轉(zhuǎn)置直接II型”它們的數(shù)學(xué)等價但數(shù)值特性和內(nèi)存占用略有差別。在MCU上我通常用直接I型或轉(zhuǎn)置直接II型因?yàn)樗鼈兊臓顟B(tài)變量天然是信號延遲線上的數(shù)值便于初始化和調(diào)試。3.2 第一個能用單二階節(jié)的C語言實(shí)現(xiàn)在MCU上用C寫IIR濾波器我建議先做成一個簡單的結(jié)構(gòu)體加處理函數(shù)不要一開始就寫很復(fù)雜的級聯(lián)框架。下面這個例子就是最典型的單二階節(jié)實(shí)現(xiàn)也就是俗稱的biquad雙二階節(jié)typedef struct { float b0, b1, b2; // 分子系數(shù) float a1, a2; // 分母系數(shù) float x1, x2; // 輸入歷史 float y1, y2; // 輸出歷史 } Biquad; float biquad_process(Biquad *f, float in) { float out f-b0 * in f-b1 * f-x1 f-b2 * f-x2 - f-a1 * f-y1 - f-a2 * f-y2; // 狀態(tài)更新 f-x2 f-x1; f-x1 in; f-y2 f-y1; f-y1 out; return out; }這段代碼的核心就是差分方程的逐樣本執(zhí)行。我特意把系數(shù)a0省略了因?yàn)樵O(shè)計(jì)工具輸出的a0通常已經(jīng)是1如果遇到a0不等于1的情況一定要先做歸一化再填入這個結(jié)構(gòu)體。我在給板子移植代碼的時候踩過一次直接從MATLAB導(dǎo)出的系數(shù)里a0 1.05沒歸一化就填進(jìn)去了結(jié)果濾波器增益整體偏移響了好久的困惑才排查到問題。同時要注意結(jié)構(gòu)體里的狀態(tài)值x1/x2/y1/y2在濾波器運(yùn)行前應(yīng)該清零。這個“清狀態(tài)”的步驟有時候比算法本身還關(guān)鍵——如果不清零上電后可能會有一段未知的瞬態(tài)輸出尤其在控制回路里這種瞬態(tài)可能導(dǎo)致執(zhí)行機(jī)構(gòu)亂動一下風(fēng)險(xiǎn)不小。3.3 與STM32硬件集成的注意事項(xiàng)在STM32上跑這個過濾器不同的人有不同的集成路徑我自己的習(xí)慣是放在ADC的中斷回調(diào)或者DMA傳輸完成回調(diào)里每采集到一個樣本就調(diào)用一次biquad_process。這樣濾波是逐樣本實(shí)時完成的輸出可以直接給到后續(xù)的PID、FFT或者顯示刷新。實(shí)時性和代碼效率是兩個需要同時考慮的問題。在Cortex-M4及以上內(nèi)核中CMSIS-DSP庫提供了arm_biquad_cascade_df1_f32函數(shù)它對多階IIR濾波做了優(yōu)化利用SIMD指令可以讓多個濾波器并行執(zhí)行。如果項(xiàng)目里只是簡單的一兩個二階節(jié)手寫的biquad代碼就夠了沒必要引入整個DSP庫但如果要處理多通道音頻或者高階數(shù)濾波強(qiáng)烈建議換到CMSIS-DSP的級聯(lián)接口省下的不只是代碼量更是實(shí)打?qū)嵉腃PU占用率。還有一點(diǎn)關(guān)系到ADC數(shù)據(jù)處理的細(xì)節(jié)讀取ADC值之后通常需要先減去直流偏置再做濾波。要是對原始ADC碼值直接濾波那個恒定的直流分量會在濾波結(jié)果里保留導(dǎo)致輸出一直偏在一個非零基線上。對于需要判斷閾值或者計(jì)算有效值的應(yīng)用這個直流偏置會讓所有后續(xù)判斷都出錯。我一般在濾波前做一次簡單的直流扣除或者輸入本來就是交流耦合信號就不用擔(dān)心。3.4 定點(diǎn)化MCU上無法回避的精度問題很多STM32型號沒有FPU或者FPU頻率低、不適合大量浮點(diǎn)運(yùn)算。這時候就要考慮把浮點(diǎn)濾波器轉(zhuǎn)成定點(diǎn)實(shí)現(xiàn)最常用的格式是Q15和Q31對應(yīng)16位和32位定點(diǎn)。定點(diǎn)化的核心思想是把系數(shù)和信號都乘以一個縮放因子放大成整數(shù)來做乘法累加最后再縮回去。比如Q15格式就是把浮點(diǎn)數(shù)乘以32768再取整。但這里有兩個必須注意的點(diǎn)系數(shù)放大后會引入量化誤差。原本0.9999的系數(shù)可能量化成0.9997這在反饋回路里可能讓極點(diǎn)位置發(fā)生微小偏移。階數(shù)高的時候這種微小偏移會累積甚至導(dǎo)致不穩(wěn)定。這也是我為什么強(qiáng)烈推薦SOS結(jié)構(gòu)的原因——每一級只有2階極點(diǎn)離單位圓的敏感度低得多。中間乘加的溢出問題。biquad內(nèi)部有5次乘法和4次加法如果用Q15做乘法結(jié)果需要32位來保存中間值。如果你每個數(shù)都是Q15那乘出來的積是Q30再加上另一個Q15就需要注意數(shù)據(jù)寬度。很多人在定點(diǎn)化的時候栽在這里要么丟精度要么溢出。如果MCU沒有FPU但又不想寫定點(diǎn)我有一條折中建議用float32試試。Cortex-M4以上的內(nèi)核雖然有FPU但即使沒有硬件浮點(diǎn)用軟件浮點(diǎn)實(shí)現(xiàn)2階IIR的耗時通常在幾十微秒級別對采樣率不高比如1kHz的應(yīng)用完全夠用。真正必須定點(diǎn)化的場景是那種采樣率幾十kHz、又要在中斷里干很多活的場合。這時候再耐心做定點(diǎn)不要一上來就把自己繞進(jìn)位寬地獄。4. SOS矩陣讓高階濾波器穩(wěn)定落地的關(guān)鍵4.1 高階直接型的數(shù)值穩(wěn)定性隱患先做個思想實(shí)驗(yàn)。想象一個10階IIR濾波器直接用傳遞函數(shù)分子分母那一堆系數(shù)去實(shí)現(xiàn)相當(dāng)于一個10階的反饋系統(tǒng)。任何微小的系數(shù)量化誤差都像在一根長竹竿的頂端加重量——竹竿越高頂端輕輕一晃底部就產(chǎn)生巨大的偏差。極點(diǎn)分布對系數(shù)誤差的敏感度跟階數(shù)成正比階數(shù)越高單位圓附近的極點(diǎn)越容易被推出圓外。這就是為什么高階IIR濾波器直接用“直接型”結(jié)構(gòu)幾乎必出問題的根本原因。你可以在MATLAB里看理論頻率響應(yīng)畫得完美無缺但把同樣的系數(shù)寫進(jìn)32位定點(diǎn)MCU跑出來的可能就是自激振蕩的噪聲。我用MATLAB做過一個實(shí)驗(yàn)一個10階巴特沃斯低通系數(shù)保留6位小數(shù)后直接實(shí)現(xiàn)極點(diǎn)位置跟原始設(shè)計(jì)差得不算大但一對共軛極點(diǎn)已經(jīng)落到了單位圓外系統(tǒng)響應(yīng)變成增長振蕩。這個實(shí)驗(yàn)特別直觀地說明不是設(shè)計(jì)的問題是實(shí)現(xiàn)結(jié)構(gòu)的問題。4.2 什么是SOS矩陣二階節(jié)的級聯(lián)SOS是Second-Order Sections的縮寫中文通常叫二階節(jié)級聯(lián)。它的核心思想是把一個N階的傳遞函數(shù)分解成N/2個2階系統(tǒng)的級聯(lián)每個2階系統(tǒng)用biquad結(jié)構(gòu)實(shí)現(xiàn)每個biquad的系數(shù)單獨(dú)歸一化。SOS矩陣的每一行長這樣[ b0, b1, b2, a0, a1, a2 ]其中a0通常等于1工具會在輸出時自動歸一化。比如一個4階濾波器輸出的SOS矩陣是2行6列就代表兩個級聯(lián)的biquad。第一個biquad的輸出喂給第二個biquad的輸入串聯(lián)完成整個濾波。從數(shù)值穩(wěn)定性角度看這個分解的意義在于每一級只承受2階的極點(diǎn)敏感度即使系數(shù)有量化誤差影響也限制在本級不會跨級級聯(lián)放大。極點(diǎn)從“一根高竹竿的一端”變成了“幾截短竹竿的連接”穩(wěn)定裕度和抗量化能力都大幅提升。此外SOS結(jié)構(gòu)還有一個工程上的好處單位“一節(jié)”就是最基本的biquad代碼結(jié)構(gòu)統(tǒng)一、調(diào)試方便每一個節(jié)的輸入輸出都可以單獨(dú)斷點(diǎn)查看過濾器到底哪一級出了問題一目了然。4.3 從普通系數(shù)轉(zhuǎn)SOS以及實(shí)際代碼結(jié)構(gòu)設(shè)計(jì)工具通常都支持直接輸出SOS。用SciPy的時候在ellip、butter這些函數(shù)里指定outputsos就可以了如果手里只有普通的分子分母系數(shù)b和a數(shù)組也可以用scipy.signal.tf2sos做轉(zhuǎn)換from scipy.signal import butter, tf2sos, sosfilt b, a butter(4, 0.2, outputba) sos tf2sos(b, a) print(sos)在C代碼里SOS級聯(lián)實(shí)現(xiàn)其實(shí)非常直觀——就是把上一節(jié)的輸出作為下一節(jié)的輸入逐節(jié)調(diào)用同一個biquad_process函數(shù)。我一般用一個數(shù)組保存所有SOS行用循環(huán)依次處理typedef struct { uint8_t num_sections; float coeffs[MAX_SECTIONS][6]; float state[MAX_SECTIONS][4]; // x1, x2, y1, y2 } SosFilter; float sos_process(SosFilter *f, float in) { float out in; for (int i 0; i f-num_sections; i) { out biquad_process((Biquad){...}, out); // 使用第i組系數(shù) } return out; }需要注意的一個細(xì)節(jié)是SOS的級聯(lián)順序。級聯(lián)順序不是隨便排的不同排序會影響數(shù)值精度。通常原則是把Q因子較高帶寬窄、諧振峰尖銳的節(jié)放在前面把增益較大的節(jié)放在前面也有助于信噪比。SciPy的sosfilt內(nèi)部會自動排序但如果你手動處理系數(shù)建議參考這個原則。另外SOS各節(jié)之間可以穿插分配濾波器總增益避免某一節(jié)固定增益特別大導(dǎo)致中間信號飽和——這個問題在定點(diǎn)實(shí)現(xiàn)里格外突出。5. 容易踩的坑與排查技巧實(shí)錄5.1 現(xiàn)象輸出爆炸先查極點(diǎn)而不是查算法我第一次把IIR濾波器燒進(jìn)板子的時候內(nèi)心是有點(diǎn)慌的——上電后輸出直接變成接近滿幅度的振蕩用示波器一看整個波形都在瘋狂抖動。當(dāng)時第一反應(yīng)是代碼寫錯了來回檢查了好幾遍biquad函數(shù)的乘加順序結(jié)果都沒問題。后來仔細(xì)一查是設(shè)計(jì)的時候用了10階直接型系數(shù)經(jīng)過去浮點(diǎn)轉(zhuǎn)換之后精密值流失極點(diǎn)跑出了單位圓。這種發(fā)散問題的排查套路現(xiàn)在已經(jīng)固定了先用PC端MATLAB或者Python把固定系數(shù)的極點(diǎn)畫出來看是不是都在單位圓內(nèi)系數(shù)量化之后再看一次極點(diǎn)位置。如果極點(diǎn)沒問題再查代碼里是否忘記先歸一化a0、是否初始化了狀態(tài)變量、是否加了錯誤的負(fù)號。這個順序基本能定位90%以上的發(fā)散問題?？偨Y(jié)起來就是IIR濾波器只要輸出不對勁永遠(yuǎn)先懷疑數(shù)值和極點(diǎn)不要懷疑硬件和主循環(huán)這是這個領(lǐng)域最值錢的一條經(jīng)驗(yàn)。5.2 現(xiàn)象截止頻率偏移記得雙線性變換的預(yù)畸變另一個高頻問題是按設(shè)計(jì)算出來的截止頻率是100Hz實(shí)測出來變成95Hz或者105Hz而且偏差可能隨截止頻率升高而變大。在我遇到過的大多數(shù)情況下原因出在模擬原型變換到數(shù)字域時沒有做“頻率預(yù)畸變”。IIR濾波器通常是從模擬濾波器原型巴特沃斯、切比雪夫、橢圓出發(fā)通過雙線性變換映射到數(shù)字域的。雙線性變換會把模擬頻率軸的無限范圍壓縮到數(shù)字頻率軸上的0到π之間這種頻率壓縮是非線性的如果不做補(bǔ)償最終數(shù)字濾波器的截止頻率就會偏移。很多時候設(shè)計(jì)工具已經(jīng)自動做了預(yù)畸變但如果你手寫轉(zhuǎn)換過程或者從模擬設(shè)計(jì)表直接查系數(shù)就容易踩雷。工具選擇錯了也會導(dǎo)致看起來“偏移”。比如設(shè)計(jì)采樣率是48kHz結(jié)果你按44.1kHz的歸一化頻率設(shè)計(jì)了那實(shí)際效果當(dāng)然對不上。如果排除代碼問題后頻率還是有偏差先檢查你的設(shè)計(jì)工具用的采樣率、歸一化方式再檢查預(yù)畸變處理大概率能解決。5.3 現(xiàn)象起始瞬態(tài)振蕩還有一個看起來嚇人但其實(shí)很常見的情況濾波器一跑起來輸出的前幾十個樣本會出現(xiàn)很大的過沖感覺像“爆炸”但過了這段之后就恢復(fù)正常了。這不是濾波器不穩(wěn)定而是起始瞬態(tài)響應(yīng)——濾波器狀態(tài)變量從零開始相當(dāng)于給濾波器輸入了一個階躍信號系統(tǒng)自然會產(chǎn)生瞬態(tài)響應(yīng)。在實(shí)時控制應(yīng)用里這種瞬態(tài)可能讓執(zhí)行機(jī)構(gòu)產(chǎn)生一次意外的運(yùn)動必須處理。我的做法是在系統(tǒng)啟動階段先運(yùn)行一段時間的濾波但丟棄這段時間的輸出或者把狀態(tài)初始化到輸入信號的平均值附近減少階躍感。在離線處理場景可以用scipy.signal.filtfilt做零相位濾波但它內(nèi)部也會處理邊緣效應(yīng)原理是反向再濾波一次相位為零但計(jì)算量翻倍。5.4 現(xiàn)象信號經(jīng)過IIR后幅值不對有時候?yàn)V波器跑得“很穩(wěn)”但仔細(xì)對比輸入輸出某個頻段的信號不是大了就是小了。這個現(xiàn)象往往是設(shè)計(jì)規(guī)格和實(shí)際需求不匹配造成的而不是代碼問題。比如你做的是Butterworth 1dB帶寬的低通但在通帶邊緣1kHz-3dB你卻拿設(shè)計(jì)規(guī)格里那條“-1dB頻點(diǎn)”跟“-3dB頻點(diǎn)”混淆了自然會覺得不對。要規(guī)避這類問題最好是拿到濾波器后先把頻率響應(yīng)打印一下在0Hz到奈奎斯特頻率范圍內(nèi)畫一條完整的響應(yīng)曲線確認(rèn)不同頻點(diǎn)的增益值符合預(yù)期。高通的低頻、帶通的兩側(cè)滾降尤其容易讓人產(chǎn)生誤解因?yàn)椴煌瑸V波器定義“通帶頻率”的方式不同有的按-3dB有的按阻帶邊緣不做校驗(yàn)就是給自己挖坑。我把這幾年在IIR濾波器上遇到過的典型問題整理成一張速查表方便直接對照現(xiàn)象可能原因快速處置輸出發(fā)散/自激振蕩極點(diǎn)跑出單位圓、系數(shù)未歸一化、定點(diǎn)溢出畫極點(diǎn)圖檢查極點(diǎn)位置先改用雙精度浮點(diǎn)驗(yàn)證截止頻率偏移未做頻率預(yù)畸變、設(shè)計(jì)工具采樣率設(shè)置錯誤檢查設(shè)計(jì)參數(shù)改用工具自動預(yù)畸變起始幅度過沖狀態(tài)變量初始化為零產(chǎn)生的瞬態(tài)啟動時丟棄前N個樣本或初始化狀態(tài)到穩(wěn)態(tài)值輸出幅度整體偏低/偏高增益分配錯誤、a0未歸一化、通帶定義混淆打印頻率響應(yīng)曲線驗(yàn)證0Hz通帶增益是否為1定點(diǎn)實(shí)現(xiàn)噪聲大系數(shù)量化誤差大、中間累加溢出改用更高的定點(diǎn)位寬使用SOS級聯(lián)結(jié)構(gòu)5.5 排查技巧從“先畫圖”到“分段定點(diǎn)”關(guān)于IIR的調(diào)試我有一個堅(jiān)持了很多年的習(xí)慣任何濾波器改動第一步一定是先畫頻率響應(yīng)圖第二步是給一個已知信號正弦波、方波、階躍看時域輸出第三步才連接到真實(shí)數(shù)據(jù)源上。三步走下來80%的問題能在接入真實(shí)系統(tǒng)前暴露剩下的20%再配合逐級斷點(diǎn)、打印狀態(tài)變量等手段也能快速定位。分段定點(diǎn)這個技巧也分享給大家。定點(diǎn)化之前先用浮點(diǎn)在PC端完整跑一遍系統(tǒng)把每一級SOS節(jié)點(diǎn)的輸出最大值記錄下來。這些值就是定點(diǎn)化的縮放依據(jù)——該級要用多大Q值不會溢出、不會浪費(fèi)精度都有依據(jù)了。分段定點(diǎn)比整體定點(diǎn)效果好得多因?yàn)椴煌?jié)之間信號動態(tài)范圍差異可能很大用統(tǒng)一的Q值意味著低幅度節(jié)浪費(fèi)精度、高幅度節(jié)面臨溢出風(fēng)險(xiǎn)。寫在最后回頭看IIR濾波器這條路最讓我感慨的一點(diǎn)是很多項(xiàng)目失敗不是設(shè)計(jì)思路錯了而是被“差一點(diǎn)點(diǎn)”的系數(shù)精度、“差一點(diǎn)點(diǎn)”的實(shí)現(xiàn)結(jié)構(gòu)、或者“差一點(diǎn)點(diǎn)”的調(diào)試方法拖垮的。IIR本身是個非常成熟的技術(shù)從巴特沃斯到SOS每一步都有明確的數(shù)學(xué)基礎(chǔ)和實(shí)踐路徑但真正把這條路走通、走穩(wěn)還是需要一點(diǎn)一點(diǎn)踩坑換來的手感。我個人現(xiàn)在做濾波器項(xiàng)目的習(xí)慣是需求拿過來先寫規(guī)格然后PC端設(shè)計(jì)校驗(yàn)確認(rèn)無誤再落到MCU代碼全部用SOS結(jié)構(gòu)浮點(diǎn)優(yōu)先、定點(diǎn)按需最后一定要做一次長時間運(yùn)行測試。這個流程看起來繁瑣但它幫我擋掉了太多半夜調(diào)板子的尷尬。希望這篇從設(shè)計(jì)原理到工程落地的經(jīng)驗(yàn)梳理能讓你在IIR濾波器上少走幾條我當(dāng)年走過的彎路。如果你在STM32或者音頻項(xiàng)目里遇到過跟IIR相關(guān)的怪問題歡迎沿著文中的思路排查一遍——很多時候答案就藏在那張速查表里。