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