1. 項目緣起當“固定”的濾波器遇上“變化”的世界在數(shù)字信號處理的世界里我們常常面臨一個經(jīng)典的矛盾設(shè)計的靈活性與實現(xiàn)的效率。就拿一個最常見的需求來說——你需要一個截止頻率可變的低通濾波器。最直接的想法是什么沒錯每次頻率參數(shù)一變我就重新計算一遍濾波器的系數(shù)然后更新到我的濾波器結(jié)構(gòu)中。這在離線處理或者對實時性要求不高的場合或許可行但一旦放到FPGA、DSP或者高性能嵌入式系統(tǒng)中頻繁地重新計算并加載一組可能多達幾十甚至上百個的濾波器系數(shù)比如一個高階FIR濾波器帶來的計算開銷和延遲往往是不可接受的。尤其是在軟件無線電、雷達信號處理、音頻效果器這類對實時性要求極高的領(lǐng)域參數(shù)需要連續(xù)、平滑地調(diào)整這種“硬切換”系數(shù)的方式會引入信號的不連續(xù)產(chǎn)生可聞的咔嗒聲或影響系統(tǒng)性能。那么有沒有一種方法能讓濾波器的某個關(guān)鍵參數(shù)比如截止頻率、分數(shù)延遲像擰旋鈕一樣連續(xù)可調(diào)而無需在每次調(diào)整時都大動干戈地重新計算所有系數(shù)呢這就是Farrow結(jié)構(gòu)濾波器要解決的核心問題。它本質(zhì)上是一種實現(xiàn)可變分數(shù)延遲濾波器或可變參數(shù)濾波器的高效結(jié)構(gòu)。我第一次在項目中接觸它是在為一個多速率采樣系統(tǒng)設(shè)計一個采樣率轉(zhuǎn)換模塊時傳統(tǒng)的多項式插值方法在精度和復(fù)雜度上難以平衡而Farrow結(jié)構(gòu)提供了一種優(yōu)雅的折中方案。今天我們就來徹底拆解這個在工程上極具魅力的結(jié)構(gòu)從它要解決什么問題到它的核心原理再到如何一步步設(shè)計并實現(xiàn)它最后分享一些在硬件實現(xiàn)中容易踩的坑。2. Farrow結(jié)構(gòu)的核心思想將參數(shù)變化“固化”到結(jié)構(gòu)里要理解Farrow結(jié)構(gòu)我們得先忘掉那些復(fù)雜的公式從一個更直觀的視角來看。想象一下一個濾波器的輸出是其系數(shù)與輸入信號卷積的結(jié)果。如果這個濾波器的某個特性如群延遲需要連續(xù)變化傳統(tǒng)做法是讓系數(shù)本身成為這個變化參數(shù)的函數(shù)。也就是說每一個濾波器系數(shù)h[k]不再是固定值而是一個關(guān)于可變參數(shù)d通常代表分數(shù)延遲范圍在0到1之間的函數(shù)h[k] f_k(d)。Farrow結(jié)構(gòu)的巧妙之處在于它對這個函數(shù)f_k(d)做了一個關(guān)鍵的假設(shè)和簡化每個系數(shù)關(guān)于可變參數(shù)d的變化可以用一個低階多項式來近似。也就是說h[k] ≈ c_{k,0} c_{k,1} * d c_{k,2} * d^2 ... c_{k,L} * d^L這里L是多項式的階數(shù)c_{k,l}就是我們需要預(yù)先計算并固定下來的“子濾波器”系數(shù)?？吹搅藛嶙兓牟糠謉及其冪次被抽離出來了而需要存儲和參與實時卷積運算的變成了固定不變的系數(shù)c_{k,l}。這個思想帶來了巨大的優(yōu)勢實時性當需要改變延遲d時我們不再需要重新計算或加載任何濾波器系數(shù)c_{k,l}。我們只需要根據(jù)新的d值實時計算出一組“權(quán)重”(1, d, d^2, ..., d^L)然后用這組權(quán)重去組合那些固定的子濾波器的輸出。結(jié)構(gòu)化與模塊化整個濾波器可以被分解為(L1)個并行的、系數(shù)固定的子濾波器每個對應(yīng)多項式的一項以及一個多項式求值即權(quán)重組合模塊。這種結(jié)構(gòu)非常規(guī)整特別適合用硬件如FPGA進行并行流水線實現(xiàn)。設(shè)計靈活性我們可以通過設(shè)計多項式階數(shù)L和子濾波器系數(shù)c_{k,l}來權(quán)衡逼近精度、計算復(fù)雜度和濾波器性能。那么這個結(jié)構(gòu)具體長什么樣呢一個典型的、用于實現(xiàn)分數(shù)延遲的Farrow結(jié)構(gòu)框圖如下所示我們以三次多項式即L3為例輸入 x[n] | |---- 固定子濾波器 H0(z) (系數(shù)為 c_{k,0}) ---- 乘 1 (即 d^0) ---\ |---- 固定子濾波器 H1(z) (系數(shù)為 c_{k,1}) ---- 乘 d ----- 加法器 ---- 輸出 y[n] |---- 固定子濾波器 H2(z) (系數(shù)為 c_{k,2}) ---- 乘 d^2 -----/ ---- 固定子濾波器 H3(z) (系數(shù)為 c_{k,3}) ---- 乘 d^3 -----/關(guān)鍵點H0(z),H1(z), ...,H3(z)都是普通的、系數(shù)固定的FIR濾波器。它們的輸出分別乘以d^0,d^1,d^2,d^3然后求和就得到了最終具有分數(shù)延遲d的輸出y[n]。d可以在每個采樣時刻動態(tài)改變而所有子濾波器的系數(shù)是焊死在硬件或代碼里的。3. 從零開始Farrow濾波器系數(shù)設(shè)計方法論理解了思想下一步就是如何得到那些固定的子濾波器系數(shù)c_{k,l}。這是Farrow濾波器設(shè)計的核心。設(shè)計目標通常是讓整個濾波器在感興趣的頻帶內(nèi)其頻率響應(yīng)盡可能地逼近一個理想的、延遲為d的分數(shù)延遲器。理想分數(shù)延遲器的頻率響應(yīng)是H_ideal(e^{jω}) e^{-jωd}它具有線性相位群延遲恒為d。設(shè)計方法主要有兩大類基于最小二乘誤差準則和基于拉格朗日插值。這里我重點講工程上最常用、也相對直觀的最小二乘設(shè)計法。3.1 設(shè)計問題建模假設(shè)我們設(shè)計一個長度為N階數(shù)為N-1的FIR濾波器用于近似分數(shù)延遲d。其傳遞函數(shù)為H(z, d) Σ_{k0}^{N-1} h[k] * z^{-k}而根據(jù)Farrow結(jié)構(gòu)我們假設(shè)h[k]是d的L階多項式h[k] Σ_{l0}^{L} c_{k,l} * d^l我們的目標是在d ∈ [0, 1]和ω ∈ [0, απ]α是帶寬因子比如0.8表示利用80%的奈奎斯特帶寬的范圍內(nèi)讓實際響應(yīng)H(e^{jω}, d)逼近理想響應(yīng)e^{-jωd}。定義一個誤差函數(shù)E(ω, d) H(e^{jω}, d) - e^{-jωd}設(shè)計目標就是找到一組系數(shù)c_{k,l}使得在定義的(ω, d)區(qū)域上誤差E的某種范數(shù)如平方誤差的積分最小。這是一個標準的優(yōu)化問題。3.2 具體設(shè)計步驟與MATLAB示例雖然推導(dǎo)過程涉及積分和矩陣運算但幸運的是我們可以借助MATLAB等工具來高效完成。下面是一個手把手的步驟步驟一確定設(shè)計參數(shù)N: 主FIR濾波器長度抽頭數(shù)。越大逼近精度越高計算量也越大。L: 多項式階數(shù)。越高對參數(shù)d變化的擬合能力越強但結(jié)構(gòu)也更復(fù)雜。通常L3三次是一個很好的折中。alpha: 歸一化帶寬0到1之間。我們只關(guān)心這個帶寬內(nèi)的性能。設(shè)計網(wǎng)格在d和ω的二維平面上需要劃分密集的網(wǎng)格點來進行離散化優(yōu)化。例如d從0到1步進0.01ω從0到alpha*pi步進pi/500。步驟二構(gòu)建最小二乘問題對于每一個網(wǎng)格點(ω_i, d_j)我們可以寫出H(e^{jω_i}, d_j) Σ_{k0}^{N-1} (Σ_{l0}^{L} c_{k,l} * d_j^l) * e^{-jω_i k}令C是一個將所有系數(shù)c_{k,l}按特定順序例如先按k再按l排列成的列向量。那么對于所有網(wǎng)格點我們可以建立一個巨大的線性方程組A * C ≈ B其中A矩陣的每一行對應(yīng)一個(ω_i, d_j)點其元素由d_j^l * e^{-jω_i k}構(gòu)成B向量對應(yīng)每個點的理想響應(yīng)值e^{-jω_i d_j}。步驟三求解系數(shù)這是一個超定線性方程組我們用最小二乘法求解C (A^H * A) \ (A^H * B)這里^H表示共軛轉(zhuǎn)置\是MATLAB中的左除運算符用于求解最小二乘解。步驟四驗證與評估求解出系數(shù)C后將其重新排列成N x (L1)的矩陣C_mat其中第k行第l列就是c_{k,l}。 然后我們可以固定一個d值如0.5用C_mat生成對應(yīng)的FIR系數(shù)h C_mat * [1; d; d^2; ... d^L]繪制其頻率響應(yīng)并與理想延遲e^{-jωd}比較。掃描d從0到1觀察濾波器幅頻響應(yīng)應(yīng)接近1和群延遲應(yīng)接近d的變化是否平滑。實操心得在MATLAB中構(gòu)建矩陣A時一定要注意索引和維度的對應(yīng)關(guān)系這是最容易出錯的地方。一個技巧是先用循環(huán)寫一個清晰但低效的版本確保邏輯正確后再嘗試用向量化操作meshgrid,kron等進行加速。另外帶寬因子alpha不要設(shè)得太大如0.95以上因為逼近理想延遲器在頻帶邊緣非常困難強求會導(dǎo)致通帶內(nèi)紋波增大。通常0.8~0.9是穩(wěn)健的選擇。4. 超越分數(shù)延遲Farrow結(jié)構(gòu)的變體與應(yīng)用擴展雖然Farrow結(jié)構(gòu)最初是為可變分數(shù)延遲而生但其“用固定子濾波器組合實現(xiàn)參數(shù)可變”的思想可以被推廣到更廣泛的可變參數(shù)濾波器設(shè)計中。這才是它真正強大和有趣的地方。4.1 可變截止頻率濾波器這是另一個非常實用的場景。假設(shè)我們需要一個截止頻率fc可變的低通濾波器。理想情況下濾波器的系數(shù)應(yīng)該是fc的函數(shù)。我們可以借鑒Farrow思想將每個系數(shù)h[k]用關(guān)于歸一化截止頻率ff fc / fsfs為采樣率的多項式來近似h[k](f) ≈ Σ_{l0}^{L} c_{k,l} * f^l設(shè)計過程與分數(shù)延遲器類似但目標函數(shù)變了。此時我們需要讓設(shè)計出的濾波器在通帶[0, f]內(nèi)響應(yīng)接近1在阻帶[fΔ, 0.5]內(nèi)響應(yīng)接近0Δ是過渡帶。這同樣可以轉(zhuǎn)化為一個在(ω, f)二維區(qū)域上的最小二乘優(yōu)化問題只是誤差權(quán)重函數(shù)需要精心設(shè)計例如在通帶和阻帶賦予高權(quán)重在過渡帶賦予低權(quán)重。實現(xiàn)結(jié)構(gòu)和經(jīng)典Farrow結(jié)構(gòu)一模一樣只是輸入?yún)?shù)從延遲d變成了截止頻率f子濾波器系數(shù)c_{k,l}是針對截止頻率變化而優(yōu)化的。4.2 應(yīng)用于采樣率轉(zhuǎn)換多項式插值這是Farrow結(jié)構(gòu)最早、也是最成功的應(yīng)用之一。在異步采樣率轉(zhuǎn)換中我們需要計算輸入序列在非整數(shù)采樣點上的值這本質(zhì)上就是一個分數(shù)延遲問題。例如常用的三次拉格朗日插值器其系數(shù)就是關(guān)于分數(shù)間隔μ相當于d的三次多項式。你可以驗證這些多項式系數(shù)正好可以排列成一個4x4的C_mat矩陣從而完美地用Farrow結(jié)構(gòu)實現(xiàn)一個高效的三次插值器。更一般地任何基于多項式的插值核如分段拋物線、樣條插值都可以用Farrow結(jié)構(gòu)來實現(xiàn)。這使得它在數(shù)字上下變頻、軟件無線電接收機中成為了標準模塊。4.3 其他可變參數(shù)理論上只要濾波器性能是某個參數(shù)p的平滑函數(shù)并且可以用多項式較好地近似就可以嘗試Farrow結(jié)構(gòu)。例如可變帶寬的帶通濾波器中心頻率固定帶寬可變。可變Q值的諧振器諧振頻率固定品質(zhì)因數(shù)Q可變。可變滾降因子的升余弦濾波器。設(shè)計經(jīng)驗在將這些擴展時最大的挑戰(zhàn)在于多項式近似的有效性。如果濾波器系數(shù)隨參數(shù)p的變化非常劇烈或非線性那么可能需要很高的多項式階數(shù)L才能較好近似這會急劇增加子濾波器的數(shù)量L1個和計算量。因此在決定采用Farrow結(jié)構(gòu)前一定要先分析系數(shù)隨參數(shù)變化的曲線是否相對平滑。通常在參數(shù)變化范圍較小如d在0~1之間時Farrow結(jié)構(gòu)表現(xiàn)最佳。5. 硬件實現(xiàn)考量與實戰(zhàn)中的“坑”將Farrow濾波器從算法模型搬到FPGA或ASIC上是另一個充滿細節(jié)的戰(zhàn)場。這里分享幾個我趟過的雷區(qū)。5.1 計算復(fù)雜度與資源優(yōu)化一個N抽頭、L階多項式的Farrow濾波器需要(L1)個并行的N抽頭FIR子濾波器。直接實現(xiàn)的乘法器數(shù)量是N*(L1)這看起來很大。但我們可以利用結(jié)構(gòu)特點進行優(yōu)化子濾波器合并計算觀察結(jié)構(gòu)每個輸入采樣x[n]需要與所有子濾波器的第一抽頭系數(shù)c_{0,0}, c_{0,1}, ..., c_{0,L}相乘然后分別延遲、再與下一組系數(shù)相乘。我們可以將同一抽頭位置k上的、屬于不同子濾波器的系數(shù)c_{k,0} ... c_{k,L}視為一個向量。當計算該抽頭的貢獻時我們實際上是在計算這個系數(shù)向量與權(quán)重向量[1, d, d^2, ..., d^L]^T的點積。這可以轉(zhuǎn)化為一個先乘累加、再乘以輸入信號x[n-k]的過程有時能減少乘法器數(shù)量。多項式求值優(yōu)化計算權(quán)重d^l可以使用霍納法則將Σ c_{k,l} * d^l的計算轉(zhuǎn)化為y_k c_{k,0} d*(c_{k,1} d*(c_{k,2} ... d*c_{k,L})...)這樣只需要L次乘法和L次加法而不是L次冪運算和L次乘法。系數(shù)對稱性利用如果目標頻率響應(yīng)是對稱的如線性相位那么設(shè)計出的子濾波器系數(shù)c_{k,l}也可能呈現(xiàn)某種對稱性。例如對于可變分數(shù)延遲濾波器其沖激響應(yīng)關(guān)于中心點近似對稱。利用這種對稱性可以將乘法器數(shù)量幾乎減半。5.2 動態(tài)參數(shù)更新的時序問題參數(shù)d或f是動態(tài)更新的。在硬件中這需要仔細處理更新速率d更新的速度不能超過數(shù)據(jù)處理流水線的“吞吐量”。如果d每個時鐘周期都變那么權(quán)重d^l需要每個周期重新計算并且必須確保在用到新權(quán)重的時刻數(shù)據(jù)流水線中對應(yīng)的是新的輸入數(shù)據(jù)。通常d的更新速率遠低于采樣率。同步新的d值必須與輸入數(shù)據(jù)流正確同步。一個常見的做法是使用一個“參數(shù)有效”信號當d更新時該信號拉高一個周期標志著從此之后進入濾波器的數(shù)據(jù)將使用新的d值。這需要在數(shù)據(jù)路徑上插入相應(yīng)的控制邏輯。5.3 有限字長效應(yīng)與精度管理這是硬件實現(xiàn)中最容易出問題的地方尤其是在定點設(shè)計中。系數(shù)量化設(shè)計得到的c_{k,l}通常是高精度浮點數(shù)。我們需要將其量化為定點數(shù)如Q格式。量化會引入誤差可能導(dǎo)致頻率響應(yīng)偏離設(shè)計目標甚至不穩(wěn)定。必須進行充分的仿真掃描不同的量化位寬如12位、16位、18位觀察通帶紋波、阻帶衰減等關(guān)鍵指標的變化找到滿足性能要求的最小位寬。中間結(jié)果位寬擴展在計算d^l以及子濾波器乘累加的過程中中間結(jié)果的動態(tài)范圍會擴大。例如計算d^2假設(shè)d是16位小數(shù)可能需要32位來保證精度不損失。在加法樹中位寬擴展更明顯。必須為每條數(shù)據(jù)路徑仔細規(guī)劃位寬防止溢出和精度過度損失。一個安全的方法是先做全精度仿真記錄中間結(jié)果的最大最小值再據(jù)此確定定點位寬。權(quán)重計算的非線性d^l的計算尤其是高次冪對d的量化誤差非常敏感。當d接近0或1時d^l的值可能非常小量化誤差會占據(jù)主導(dǎo)?？梢钥紤]采用查找表LUT來存儲d^l的值或者使用分段線性近似等方法來計算高次冪以平衡精度和資源消耗。踩坑實錄在一次通信接收機項目中我們使用Farrow結(jié)構(gòu)做定時恢復(fù)。仿真時性能完美但上板后誤碼率總是差一點。用邏輯分析儀抓取中間數(shù)據(jù)發(fā)現(xiàn)當分數(shù)延遲d在0.1以下時輸出信號有明顯的失真。最終定位到問題計算d^3的模塊為了節(jié)省乘法器我們采用了(d*d)*d的順序計算。由于d很小d*d的結(jié)果在定點量化后直接變成了0導(dǎo)致三次項完全失效。教訓是對于可能接近零的小數(shù)運算必須保證中間結(jié)果的精度或者改變計算順序如先計算高精度浮點再量化或者采用查找表。后來我們改用一個小型的、針對d的LUT來存儲d^2和d^3的值問題得以解決。6. 性能評估與設(shè)計實例一個可變延遲線讓我們用一個具體的設(shè)計實例把前面所有的點串起來。目標設(shè)計一個用于音頻處理的可變分數(shù)延遲線要求延遲d在0到1個采樣周期內(nèi)連續(xù)可調(diào)音頻帶寬為20kHz采樣率fs48kHz因此alpha ≈ 0.833。設(shè)計參數(shù)選擇主濾波器長度N32。這是一個折中能提供較好的阻帶抑制。多項式階數(shù)L3。三次多項式足以在0-1區(qū)間平滑近似。帶寬因子alpha0.85略高于需求留有余量。采用最小二乘準則在d[0:0.02:1]和ω[0:π/1000:0.85π]的網(wǎng)格上設(shè)計。設(shè)計結(jié)果 使用MATLAB的farrow設(shè)計函數(shù)或按前述步驟自編代碼得到32x4的系數(shù)矩陣C_mat。性能評估固定延遲響應(yīng)取d0.5生成系數(shù)h C_mat * [1; 0.5; 0.25; 0.125]。繪制其頻率響應(yīng)。實測在0-20kHz通帶內(nèi)幅度起伏 0.01 dB群延遲波動 0.001個樣本非常接近理想的0.5樣本延遲。動態(tài)掃描讓d從0線性變化到1觀察濾波器幅頻響應(yīng)。在整個過程中通帶內(nèi)幅度保持平坦阻帶20kHz衰減始終大于80dB。群延遲值緊密跟隨d值變化誤差在±0.005個樣本以內(nèi)。聽感測試用一段正弦波掃頻信號通過該延遲線同時用低頻三角波調(diào)制d在0-1之間變化。輸出信號聽起來應(yīng)該是純凈的、音高無變化因為只有延遲無頻移但帶有周期性相位調(diào)制的效果。如果設(shè)計不良可能會引入可聞的諧波失真或振幅調(diào)制。本例中聽感干凈。硬件實現(xiàn)概要FPGA系數(shù)存儲32*4128個系數(shù)每個量化為18位有符號定點數(shù)存儲在ROM中。數(shù)據(jù)處理路徑輸入音頻數(shù)據(jù)x[n]為24位。設(shè)計4個并行的32抽頭FIR濾波器H0-H3。利用系數(shù)的近似對稱性每個FIR采用轉(zhuǎn)置結(jié)構(gòu)乘法器數(shù)量減半約需(32/2)*4 64個乘法器。參數(shù)d為16位無符號小數(shù)0x0000~0xFFFF 對應(yīng) 0~1。權(quán)重計算模塊使用一個時鐘周期計算d2 d*d下一個周期計算d3 d2*d均為32位中間結(jié)果。同時使用一個三級流水線的霍納法則單元計算每個抽頭k的加權(quán)和sum_k c_k0 d*(c_k1 d*(c_k2 d*c_k3))。最終輸出y[n]為24位由所有sum_k * x[n-k]的累加和得到。資源與性能在Xilinx Artix-7 FPGA上占用約70個DSP slice最高運行時鐘頻率可達150MHz遠高于48kHz的音頻采樣率有充足資源進行多通道處理。通過這個實例可以看到一個精心設(shè)計的Farrow濾波器能夠在參數(shù)連續(xù)變化時保持卓越且穩(wěn)定的性能同時滿足實時硬件處理的要求。它不再是教科書上一個抽象的數(shù)學結(jié)構(gòu)而是一個能解決實際工程難題的強力工具。