:SOS矩陣與直接I型實戰(zhàn)指南)
開場別被老技術(shù)這三個字騙了IIR濾波器無限脈沖響應(yīng)濾波器大學(xué)數(shù)字信號處理課上最勸退的那一章也是我工作這些年里用得最多的濾波器。你可能覺得FIR才是萬金油什么場合都能套一個窗函數(shù)上去但我告訴你在很多資源受限的嵌入式場景里IIR才是真正的保命方案。為什么一句話同樣的濾波效果IIR用的階數(shù)更低算得快省內(nèi)存。一個3階的巴特沃斯低通性能大致能抵得上十幾階甚至幾十階的FIR。對跑在STM32這種Cortex-M內(nèi)核上的程序來說這差別不是一點點是實實在在的算力開銷和RAM占用。這篇文章我不會從Z變換的嚴格定義開始講那套東西教材里多得是。我按自己實際用下來的思路來先搞明白IIR到底是個什么東西它和FIR的核心差異在哪然后講清楚SOS矩陣和直接I型這兩種最常用的實現(xiàn)方式再落到STM32上系數(shù)怎么算、代碼怎么寫、怎么避開那些害死人的坑。讀完你能直接上手至少不會再用錯結(jié)構(gòu)、算錯系數(shù)。1. IIR濾波器的本質(zhì)輸出不僅取決于輸入還取決于過去的輸出1.1 從差分方程看IIR的反饋本質(zhì)FIR濾波器的差分方程長這樣y[n] b0·x[n] b1·x[n-1] ... bN·x[n-N]注意輸出只和當(dāng)前及過去的輸入有關(guān)沒有輸出反饋。所以一個單位脈沖進去經(jīng)過N個采樣點之后輸出就歸零了脈沖響應(yīng)是有限長的這就是Finite Impulse Response名字的由來。IIR濾波器的差分方程多了一項y[n] b0·x[n] b1·x[n-1] ... bM·x[n-M] - a1·y[n-1] - a2·y[n-2] - ... - aN·y[n-N]后面這一串帶a系數(shù)的項是把過去的輸出又加權(quán)加回來這就是反饋。因為輸出會不斷反饋到輸入端所以一個單位脈沖進去理論上輸出永遠不會完全歸零雖然實際上因為有限精度會逐漸衰減到零脈沖響應(yīng)是無限長的這就是Infinite Impulse Response。這個反饋是IIR的靈魂也是它所有優(yōu)缺點的根源。反饋讓同樣的濾波效果只需要很少的系數(shù)——階數(shù)低、運算量小、內(nèi)存占用小。但反饋也帶來了穩(wěn)定性問題如果反饋系數(shù)不合適輸出可能發(fā)散直接飄到天上去。FIR是絕對穩(wěn)定的只要系數(shù)是有限值IIR則必須檢查極點位置。1.2 頻域視角IIR的優(yōu)勢從哪來IIR濾波器在頻域上可以做到極陡的過渡帶。比如你有一個50Hz的工頻干擾旁邊就是你要的100Hz信號你用FIR想把這個干擾壓下去可能需要60階甚至80階但用IIR2階到4階就差不多能做到。原因是IIR的系統(tǒng)函數(shù)可以寫成H(z) B(z) / A(z)分子多項式B(z)決定零點分母多項式A(z)決定極點。FIR只有分子相當(dāng)于只有零點IIR既有多項式又有分母可以實現(xiàn)更復(fù)雜的頻率形狀。極點可以看成是讓濾波器在某些頻率上諧振從而實現(xiàn)高增益的窄帶特性或者急劇變化的相位特征。用個生活化的類比FIR是純粹靠記住更多歷史數(shù)據(jù)來平滑結(jié)果的算法像一個記性好但反應(yīng)慢的人IIR是記得住過去的結(jié)果并據(jù)此調(diào)整當(dāng)前判斷的算法像一個經(jīng)驗豐富但偶爾會一根筋不穩(wěn)定的老手。1.3 一個看得見的例子同參數(shù)下直接對比我用一個具體的例子來說明。假設(shè)采樣率Fs 1000Hz截止頻率Fc 50Hz的低通濾波器MATLAB/Octave里用同樣的設(shè)計規(guī)格來對比。如果用切比雪夫I型IIR3階就能做到通帶波紋0.5dB、阻帶衰減40dB。但如果用FIR要達到差不多的過渡帶寬度和阻帶衰減用窗函數(shù)法設(shè)計估算需要的階數(shù)大約是N ≈ (Attenuation_dB - 8) / (2.285 × Δω)其中Δω是歸一化過渡帶寬度。算下來至少需要20多階。在STM32F103這種72MHz的MCU上每秒鐘采樣1000次每個采樣點要算20次乘加運算IIR只需要算6次乘加。高下立判。2. 從傳遞函數(shù)到實際代碼IIR的三種常見實現(xiàn)結(jié)構(gòu)這一節(jié)是實操的基石。很多新手直接拿高階IIR系數(shù)往代碼里一塞發(fā)現(xiàn)輸出全是NaN根本不知道自己錯在哪。問題往往出在實現(xiàn)結(jié)構(gòu)上。2.1 直接I型Direct Form I最直觀但也最浪費直接I型是最容易理解的實現(xiàn)方式就是把差分方程照抄成代碼y[n] b0·x[n] b1·x[n-1] ... bM·x[n-M] - a1·y[n-1] - ... - aN·y[n-N]在代碼里你需要維護兩個數(shù)組一個是輸入的歷史緩沖x_buffer一個是輸出的歷史緩沖y_buffer。處理每個新樣本的偽代碼如下// Direct Form I 偽代碼 y b0*x b1*x_buf[0] b2*x_buf[1] - a1*y_buf[0] - a2*y_buf[1]; // 更新緩沖 x_buf[1] x_buf[0]; x_buf[0] x; y_buf[1] y_buf[0]; y_buf[0] y;這樣有兩個歷史緩沖各占M和N個位置。看起來邏輯簡單但要注意兩個問題第一存儲浪費。雖然IIR階數(shù)低但每個濾波器都要維護兩個緩沖。如果是多通道的音頻處理或者多路傳感器信號內(nèi)存消耗會成倍增加。第二數(shù)值靈敏度。直接I型在系數(shù)取值范圍很大時中間結(jié)果可能非常大在定點DSP或單片機上很容易溢出。這是它最大的問題。但直接I型也有好處它對系數(shù)量化誤差的敏感度相對較低相比直接II型來說而且移植簡單、容易查錯。STM32這類帶FPU的MCU上用浮點運算實現(xiàn)直接I型性能基本不是問題出不了大亂子。2.2 直接II型Direct Form II省內(nèi)存但更矯情直接II型也叫Canonical Form是一種更節(jié)省內(nèi)存的結(jié)構(gòu)。它首先計算一個中間變量w[n]w[n] x[n] - a1·w[n-1] - ... - aN·w[n-N]然后再計算輸出y[n] b0·w[n] b1·w[n-1] ... bM·w[n-M]這樣做的好處是你只需要維護一組狀態(tài)變量w_buffer而不是兩組內(nèi)存開銷減少了一半。在很多老的DSP課程里它被大大推崇但實際工程中用得反而少原因在于它對系數(shù)誤差更敏感。特別是在把高階濾波器系數(shù)直接量化成16位定點數(shù)時直接II型的極點位置偏移可能比直接I型更大稍不留神濾波器就不符合設(shè)計規(guī)格了。2.3 級聯(lián)型SOS高階濾波器的正確打開方式如果你搜索過iir濾波器 sos 矩陣你看到的SOSSecond-Order Sections就是級聯(lián)型的標(biāo)準格式。它的核心思想是不把高階IIR濾波器作為一個整體去實現(xiàn)而是把它拆成多個二階濾波器串行級聯(lián)起來。為什么這么做因為高階多項式在數(shù)值計算上非常脆弱。假設(shè)你設(shè)計了一個8階的巴特沃斯濾波器分母多項式A(z)有8個系數(shù)這些系數(shù)的動態(tài)范圍可能非常大小到10的負幾次方大到幾百。在浮點運算中可能還好但在定點處理器上量化誤差會迅速放大極點的實際位置和理論位置偏差很大濾波器可能變成振蕩器——輸出持續(xù)抖動甚至發(fā)散。如果把8階濾波器拆成4個二階節(jié)每個二階節(jié)的系數(shù)范圍小得多、數(shù)值穩(wěn)定性好得多依次計算每一級的輸出是下一級的輸入整體效果不變但數(shù)值表現(xiàn)優(yōu)秀得多。SOS矩陣的格式通常是這樣的每一行表示一個二階節(jié)[b0, b1, b2, 1, a1, a2]舉個例子一個4階巴特沃斯低通濾波器Fs1000HzFc50Hz在MATLAB/Octave中設(shè)計后用sos函數(shù)輸出的就是一個Nx6的矩陣sos 0.0201, 0.0402, 0.0201, 1.0000, -1.5606, 0.6414 1.0000, 2.0000, 1.0000, 1.0000, -1.3643, 0.5098每一行代表一個二階節(jié)。第一行增益較低b系數(shù)小第二行增益較高b系數(shù)接近1。兩級串聯(lián)在一起整體的傳遞函數(shù)就是這兩個二階節(jié)的乘積。級聯(lián)SOS是實際工程中最推薦的選擇。TI的DSP庫、CMSIS-DSP里的arm_biquad_cascade_df1_f32函數(shù)以及各種CMSIS濾波器庫都是基于SOS級聯(lián)結(jié)構(gòu)實現(xiàn)的。3. 用SOS矩陣串聯(lián)還是直接算代碼實現(xiàn)的關(guān)鍵區(qū)別現(xiàn)在問題來了已知SOS矩陣怎么在代碼里實現(xiàn)3.1 逐級處理的實現(xiàn)方法最直接的方法是按順序處理每一級。注意每一級的輸出就是下一級的輸入處理完第一級得到中間信號再傳給第二級。偽代碼如下// 假設(shè)有numSections個二階節(jié)sos是numSections x 6的系數(shù)矩陣 // x_in是當(dāng)前采樣值state是numSections x 4的狀態(tài)緩沖 float filter_process(float x_in) { float y x_in; for (int i 0; i numSections; i) { y biquad_process(y, sos[i], state[i]); } return y; }其中biquad_process輸入一個值輸出這個二階節(jié)的輸出。它的內(nèi)部實現(xiàn)直接I型二階節(jié)float biquad_process(float x) { float y b0*x b1*state[0] b2*state[1] - a1*state[2] - a2*state[3]; state[1] state[0]; state[0] x; state[3] state[2]; state[2] y; return y; }這就是CMSIS-DSP中arm_biquad_cascade_df1_f32函數(shù)的思路。這種逐級處理的方法也有額外的調(diào)試好處你可以在任意一級輸出處加打印或者斷點確認是哪一級出了問題排錯方便很多。3.2 直接實現(xiàn)一個塞滿系數(shù)的高階濾波器為什么危險有人會想既然我有8階濾波器的全部系數(shù)為什么不直接把差分方程里所有a和b系數(shù)寫進代碼理論上可以但實際工程中我強烈不建議。原因還是之前提到的數(shù)值穩(wěn)定性問題。用MATLAB/Octave設(shè)計一個8階切比雪夫II型濾波器它的極點和零點可能非常接近單位圓任何微小的量化誤差都可能把極點推出單位圓外導(dǎo)致濾波器不穩(wěn)定。就算不推到單位圓外極點位置偏移也會讓實際頻響曲線和設(shè)計規(guī)格差異很大可能原來要求-40dB處實際上只有-32dB。還有一個隱患如果輸入的x信號太大中間各級的輸出可能很大在級聯(lián)結(jié)構(gòu)里你可以在每一級之間重新歸一化很多庫會自動做但直接實現(xiàn)的高階結(jié)構(gòu)沒有這個靈活度。如果中間變量超出浮點數(shù)范圍直接就NaN了。3.3 從MATLAB/Octave導(dǎo)出SOS矩陣的操作流程這段操作值得仔細看因為它是搜索iir濾波器 sos 矩陣的人最想找到的答案。在MATLAB/Octave中的標(biāo)準流程是設(shè)計原型濾波器比如一個5階的巴特沃斯低通采樣率Fs1000Hz截止頻率50Hz[z, p, k] butter(5, 50/(1000/2), low);注意50/(1000/2)是歸一化截止頻率數(shù)字信號處理中頻率必須用奈奎斯特頻率Fs/2歸一化這是最容易搞錯的地方。把零極點增益模型轉(zhuǎn)換成SOS矩陣[sos, g] zp2sos(z, p, k);sos就是二級節(jié)矩陣g是全局增益。這里有一個重要技巧zp2sos可以帶參數(shù)up或down指定級聯(lián)順序通常用up把最靠近單位圓的極點放在最后數(shù)值穩(wěn)定性更好[sos, g] zp2sos(z, p, k, up);把全局增益分配進第一級和第二級。許多庫的biquad結(jié)構(gòu)自帶增益系數(shù)每個節(jié)都有b0/b1/b2所以你可以把增益g乘到第一級的b系數(shù)上也可以平均分配到所有級。為了縮小每級的動態(tài)范圍通常建議g乘到第一級sos(1, 1:3) sos(1, 1:3) * g;這樣后面每一節(jié)都不需要額外乘增益。但如果g特別大或特別小建議在級間觀察信號幅度必要時手動調(diào)整增益分配。如果要在STM32上用定點實現(xiàn)還需要把系數(shù)轉(zhuǎn)換成Q格式。比如用16位Q15表示系數(shù)取值范圍-1到1之間但要注意IIR系數(shù)的負數(shù)范圍可能超過-1這需要額外處理。通常建議直接用浮點STM32F4以上帶FPU或者用32位定點而不是用16位。后面我會詳細講這個問題。4. FIR和IIR到底怎么選從幾個實際工程場景看4.1 一張表站在需求角度做對比先給一個實用對比表這是我做選型時必看的對比維度IIRFIR相同濾波效果的階數(shù)低通常2~8階高通常幾十階計算量小大內(nèi)存占用小大線性相位不支持相位非線性天然支持關(guān)于中心對稱穩(wěn)定性需要檢查極點始終穩(wěn)定數(shù)值敏感性高需小心實現(xiàn)低適合場景資源受限、實時性要求高的嵌入式音頻處理、需要無相位失真的數(shù)據(jù)采集表格里最關(guān)鍵的指標(biāo)是線性相位。4.2 相位敏感場景FIR勝出如果你的應(yīng)用是對信號做高精度測量比如振動分析、心電信號處理、音頻濾波波形的時域形狀很重要那么IIR的非線性相位可能是個問題。IIR濾波器在某些頻率點會有較大的群延遲波動信號經(jīng)過濾波器后不同頻率成分的延遲不同時域波形會發(fā)生畸變。舉個例子你用IIR濾波器處理心電信號它的P波、QRS波群、T波含有不同的頻率成分IIR濾波器會導(dǎo)致這些波的相對時間關(guān)系發(fā)生偏移醫(yī)生一看波形就覺得不對。你以為濾波后信號變干凈了實際上它已經(jīng)被扭曲了。這種情況下FIR的線性相位特性就是剛需。但如果你只是從傳感器數(shù)據(jù)里濾掉一些噪聲不關(guān)心波形的具體相位比如做溫控的PID控制里濾掉高頻抖動、電池管理系統(tǒng)里讀取電流電壓的平均值那IIR完全沒有問題。4.3 資源有限場景IIR勝出這又要說回STM32了。假設(shè)你用STM32F103主頻72MHz沒有FPU只有單精度浮點仿真庫做一個電機電流環(huán)的采樣濾波。電流環(huán)的采樣頻率可能要到10kHz甚至20kHz。如果你用FIR濾波器階數(shù)50每次采樣要做50次乘加運算在無FPU的芯片上這已經(jīng)是很重的開銷了。但如果用IIR 2階濾波器每次采樣只需要大約6次乘加運算少了一個數(shù)量級。更關(guān)鍵的是狀態(tài)緩沖只需要4個浮點數(shù)而FIR需要50個。所以在硬實時系統(tǒng)中IIR往往是唯一的合理選擇。你可能覺得反正計算量也不大但你要考慮整個系統(tǒng)的預(yù)算——中斷里除了濾波還有PID計算、通信協(xié)議處理、狀態(tài)機判斷。每多1微秒的開銷在高速控制回路里都是要命的。5. STM32上實現(xiàn)IIR從系數(shù)獲取到實測這是搜索stm32 iir濾波器直接i型 系數(shù)的人最關(guān)心的部分我直接按完整流程來。5.1 在STM32上用浮點直接I型實現(xiàn)2階IIR先提供一個可以直接用的函數(shù)這是最經(jīng)典的2階直接I型實現(xiàn)系數(shù)從MATLAB/Octave導(dǎo)出后手動填入typedef struct { float b0, b1, b2; float a1, a2; float x1, x2; // 輸入歷史 float y1, y2; // 輸出歷史 } iir_biquad_t; float iir_biquad_process(iir_biquad_t *f, float x) { float y f-b0 * x f-b1 * f-x1 f-b2 * f-x2 - f-a1 * f-y1 - f-a2 * f-y2; // 更新歷史 f-x2 f-x1; f-x1 x; f-y2 f-y1; f-y1 y; return y; }注意這里的符號約定在MATLAB中濾波器系數(shù)形式一般是y[n] b0*x[n] ... - a1*y[n-1] - ...所以代碼里面減法要對應(yīng)好。很多人就是在這里搞錯了符號導(dǎo)致濾波器完全不對。5.2 從MATLAB/Octave獲取系數(shù)并檢查以STC/STM32項目中常用到的一個50Hz工頻陷波器為例采樣率Fs500Hz希望濾除50Hz干擾??梢杂胕irnotch函數(shù)Fs 500; fo 50; bw 5; % 3dB帶寬 [b, a] iirnotch(fo/(Fs/2), bw/(Fs/2));得到類似這樣的系數(shù)b [0.97551, -1.10060, 0.97551] a [1.00000, -1.10060, 0.95102]注意這里a的第一個元素是1代碼里不需要用它但需要確認。把b0、b1、b2、a1、a2分別填到結(jié)構(gòu)體的對應(yīng)字段里就行。在把系數(shù)燒到單片機之前強烈建議先在電腦上做一次快速仿真驗證??梢栽贛ATLAB/Octave里用freqz(b, a, 1024, Fs)看頻響曲線確認陷波點位置和帶寬正確再用stepz或impz看時域響應(yīng)是否收斂。這一步花5分鐘能省去后面在板子上調(diào)試幾小時。5.3 定點還是浮點STM32上怎么選STM32F0、F1這類不帶FPU的芯片用浮點數(shù)做實時的IIR是有代價的。雖然編譯器有軟件浮點庫但乘法運算會擴展到幾十條匯編指令速度慢。這種情況下有兩種方案方案一用Q15或Q31格式實現(xiàn)在STM32上不太推薦自己搞。STM32官方庫里的arm_biquad_cascade_df1_f32是單精度浮點版本需要用帶FPU的Cortex-M4/M7/M33。CMSIS-DSP里也有arm_biquad_cascade_df1_q15和arm_biquad_cascade_df1_q31定點版本可以直接用省心很多。方案二如果你必須用STM32F1系列說實話我建議先算算實際負載再決定。一個2階IIR每秒鐘跑1000次軟件浮點也就需要幾百微秒對于慢速采樣場景其實無所謂。只有當(dāng)采樣率很高或多路濾波時才需要考慮軟件浮點性能瓶頸。方案三用定點手動實現(xiàn)系數(shù)非常不建議初學(xué)者做因為會遇到飽和、舍入誤差、極限環(huán)振蕩等問題。我個人的偏好是如果成本允許直接上STM32G4或STM32F4帶FPU用單精度浮點實現(xiàn)IIR開發(fā)速度和調(diào)試體驗遠超定點方案價格差距也就幾塊錢。5.4 一個簡易測試方法板子上怎么確認濾波器沒寫錯代碼寫完了板子跑起來了怎么確認濾波器正常工作我的做法是在濾波函數(shù)入口處注入一個已知信號觀察輸出是否和MATLAB對同一信號的仿真結(jié)果一致。具體操作在單片機里生成一個固定頻率的正弦波數(shù)組比如放在Flash里的const數(shù)組或者直接在代碼里用簡單的sin函數(shù)生成。給定一個包含50Hz低頻和300Hz高頻混合的信號經(jīng)過濾波器后觀察輸出波形是否只留下低頻部分。把UF問題在輸出端加一個調(diào)試變量通過串口或者JTAG調(diào)試器把數(shù)據(jù)點抓出來畫波形。如果輸出波形形狀和MATLAB仿真基本一致幅值允許有微小誤差濾波器就通了。這一步別跳曾經(jīng)我身邊不止一個人跳過這一步直接接到實際傳感器上結(jié)果濾波器系數(shù)里符號弄錯折騰了一天最后發(fā)現(xiàn)低頻全被濾掉了。6. 實際踩坑記錄IIR濾波器工程化的五個教訓(xùn)最后這部分是我自己踩坑踩出來的每一件都付出過時間成本。6.1 初始瞬態(tài)問題濾波器的啟動沖擊IIR濾波器有反饋所以它有一個建立時間。假設(shè)初始化時狀態(tài)變量全部清零突然輸入一個大的階躍信號輸出會有一個很大的過沖可能持續(xù)幾毫秒甚至更長。在很多控制系統(tǒng)中這個過沖會導(dǎo)致執(zhí)行機構(gòu)猛地動一下這是不可接受的。我的處理辦法是系統(tǒng)啟動時先不要馬上讓濾波器輸出控制量而是先讓傳感器信號穩(wěn)定一段時間或者在初始化時把濾波器的輸入歷史狀態(tài)賦值為當(dāng)前采樣值的合理估計讓濾波器從工作點開始。比如傳感器剛上電時輸出接近0就把x1、x2、y1、y2都初始化成0如果預(yù)期傳感器穩(wěn)定在某個偏置電壓你可以先采N個樣本求平均然后設(shè)置初始狀態(tài)。6.2 系數(shù)量化誤差16位定點的災(zāi)難我第一次在STM32F103上做IIR時天真地把浮點系數(shù)直接轉(zhuǎn)成Q15格式結(jié)果濾波器特性完全亂套。原因很簡單IIR系數(shù)的動態(tài)范圍經(jīng)常超過Q15能表示的范圍比如系數(shù)可能在0.998到1.002之間Q15量化后直接變成一個小整數(shù)極點的精度損失太嚴重。后來我改用浮點加上CMSIS-DSP的庫問題就消失了。如果確實只能用定點至少要選Q31并且考慮用級聯(lián)SOS結(jié)構(gòu)。16位定點做IIR尤其是高階IIR基本等于自找麻煩。6.3 采樣率不穩(wěn)定的影響IIR系數(shù)是基于固定采樣率設(shè)計的。如果一個系統(tǒng)用定時器中斷觸發(fā)采樣但是定時器優(yōu)先級設(shè)置不當(dāng)導(dǎo)致采樣周期抖動濾波器實際表現(xiàn)會和設(shè)計差很多。最典型的例子是你用HAL_Delay或者簡單的while循環(huán)來做ADC采樣觸發(fā)時鐘不準確采樣率漂移你以為是5kHz采樣并計算截止頻率實際可能是4.7kHz濾波器特性偏移還不算大但如果觸發(fā)被其他中斷打斷采樣間隔忽長忽短IIR的反饋會累積誤差輸出出現(xiàn)抖動。解決辦法是用硬件定時器觸發(fā)ADC采樣或者用DMA保證采樣間隔精確。這是嵌入式工程師做任何數(shù)字濾波都必須刻進DNA的規(guī)范。6.4 狀態(tài)變量精度float還是doubleCortex-M4的FPU有單精度浮點速度極快但精度只有大約7位有效數(shù)字。對于2階IIR來說這精度足夠。但對于8階以上的高階IIR單浮點可能不夠用狀態(tài)變量累積誤差會導(dǎo)致輸出噪聲增加。我的建議8階以下放心用float更高階建議用double或者保證SOS級聯(lián)順序合理。還有如果你需要極窄帶濾波器比如Q值很高的帶通或陷波器float可能不夠因為極點離單位圓太近對精度極其敏感。我的經(jīng)驗閾值是如果濾波器的極點模值超過0.999用double。6.5 濾波器穩(wěn)定性必須在真實運行范圍內(nèi)驗證很多濾波器的極點位置是在設(shè)計階段算好的但實際運行過程中信號頻率和幅度的變化可能導(dǎo)致濾波器中間變量過大出現(xiàn)飽和等非線性現(xiàn)象使得系統(tǒng)不穩(wěn)定。特別是定點實現(xiàn)中飽和溢出后狀態(tài)變量變成極大值反饋恢復(fù)不過來濾波器就卡在抱死狀態(tài)。所以穩(wěn)定性驗證不能只看理論要加到最惡劣的輸入信號進行測試比如信號滿幅值、接近截止頻率的階躍信號。如果出現(xiàn)異常首先檢查各節(jié)輸出是否觸及數(shù)據(jù)范圍上界必要時加入飽和保護邏輯。結(jié)尾一些關(guān)于IIR的個人心得IIR濾波器不是銀彈但它絕對是嵌入式工程師和算法工程師工具箱里不可缺少的一把錘子。我用它做過電機電流環(huán)的陷波、心率傳感器的噪聲抑制、電源紋波數(shù)據(jù)的平滑、音頻回聲消除的前置濾波每個場景都踩過不同的坑但核心方法論是一樣的先算清需求指標(biāo)選對階數(shù)和類型用MATLAB/Octave設(shè)計并用SOS級聯(lián)實現(xiàn)最后在板子上用已知信號驗證。如果你想從IIR開始入手我建議先從2階巴特沃斯低通開始把差分方程、直接I型實現(xiàn)、頻響驗證這套鏈路走通再嘗試高階的切比雪夫、橢圓濾波器。在這個過程中要特別留意那個符號約定的問題——MATLAB里a系數(shù)帶著負號出現(xiàn)在差分方程里很多初學(xué)者在這里翻車不是少數(shù)。最后分享一個小技巧如果你的系統(tǒng)里信號頻率范圍比較寬又擔(dān)心IIR的相位失真影響波形可以在信號鏈路里先做一次正向濾波再做一次反向濾波所謂零相位濾波但這只適合離線數(shù)據(jù)處理實時系統(tǒng)里就別想了。實時系統(tǒng)的話要么接受IIR的相位特性要么踏踏實實上FIR。搞清楚自己的需求邊界比糾結(jié)哪個技術(shù)更高級重要得多。