wisetop/tcad-bodeplot
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
9182220c1daa52b44f3d19285dfe26dd6086c1fe
tcad-bodeplot/INTEGER_FILTER_IMPLEMENTATION.md
T

19 KiB
Raw Blame History

Difference Equation Analyzer (DEA)

定點數整數差分方程式實作規格說明書

本文件旨在為嵌入式韌體、DSP 與硬體描述語言（HDL）開發工程師，提供本專案「定點數整數濾波模擬（integer_lfilter）」的完整實作細節與移植指引。

本系統所採用的定點數架構為高精度狀態變數架構 (Extended Precision State Variables)，並支援 1-Clock 高效能硬體四捨五入 (Rounding)，能有效消除直流偏移（DC Bias）並最大化訊噪比（SNR）。

1. 定點數 Q 格式 (Q-Format) 定義

在實作中，所有浮點數（實數）均以帶符號的整數進行二進位縮放表示：

\text{整數值} = \text{round}\left( \text{浮點數值} \times 2^{Q} \right)

本專案定義了以下四個獨立的 Q 格式位移參數：

參數名稱 代號 說明 格式描述
shift_in Q_{in} 輸入訊號 x[n] 的量化位元數 x_{\text{int}}[n] \in Q_{in}
shift_out Q_{out} 最終輸出訊號 y[n] 的量化位元數 y_{\text{out}}[n] \in Q_{out}
shift_b Q_b 前饋（Feedforward）係數 b 的量化位元數 b_{\text{int}}[i] \in Q_b
shift_a Q_a 回授（Feedback）係數 a 的量化位元數 a_{\text{int}}[j] \in Q_a

2. 核心架構：高精度狀態變數 (Extended Precision State Variables)

為了極小化乘加運算中的截斷誤差 (Truncation Error) 並防止 IIR 濾波器極點附近因量化產生的極限環震盪（Limit Cycles），本架構採用前饋完全不右移的策略：

  1. 狀態變數 y_{\text{hist}} 精度保持：歷史狀態變數 y_{\text{hist}}[n-j] 並非儲存低精度的最終輸出 $y_{\text{out}}$，而是直接儲存前饋與回授累加後、尚未進行輸出位移的高精度狀態值。其格式固定為：

    y_{\text{hist}} \in Q_{in + b}

  2. 免除中間位移：在前饋運算中，b_{\text{int}} \times x_{\text{int}} 的結果為 $Q_b \times Q_{in} = Q_{in+b}$，完全不需進行任何位移，即可直接累加至 sum_b 中。

  3. 回授對齊縮放（歸一化 $Q_a$）：回授運算中，係數 $a_{\text{int}} \in Q_a$，歷史狀態 $y_{\text{hist}} \in Q_{in+b}$，兩者相乘後格式為 $Q_{in+b+a}$。

    💡 重要前提與物理意義：假設 $A_0 = 1$
    濾波器公式通常已正規化使得回授係數 $A_0 = 1$。當我們使用 Q_a 將 a 係數放大時，A_0 的整數值實際上就是 $2^{Q_a}$。

    在運算完畢後，這項「因為 a 係數被放大了 Q_a 倍」而產生的額外增益必須被歸一（Normalize），否則數值會以指數級別爆炸！因此，我們將回授乘加總和除以 $A_0$（在二進位中等同於算術右移 shift_a 位元），藉此消去 Q_a 的放大效應，使其完美對齊回 Q_{in+b} 格式，才能與前饋結果相減。

    💡 關於輸入訊號 x 與前饋係數 b 的歸一化討論
    許多工程師會問：既然回授部分需要被歸一以防爆炸，那輸入訊號 $x$（放大了 Q_{in} 倍）與前饋係數 $b$（放大了 Q_b 倍）在運算過程中是否也需要被歸一化？
    答案是：完全不需要，而且我們也故意沒有做！

    這是因為：

    1. 避免資訊截斷（Truncation Noise）：如果我們在累加到 sum_b 時就進行右移歸一，會直接把低位元的小數精度給丟棄，造成無法挽回的量化雜訊。
    2. 最優化運算效能：在迴圈內部，sum_b（$Q_{in+b}$）與被歸一後的 sum_a_scaled（$Q_{in+b}$）可以直接進行減法運算。整個迴圈內部不需要對 x 與 b 進行 any 移位，這在硬體實作中極大提升了運算速度。
    3. 延遲至最終輸出對齊：輸入與前饋所累積下來的放大倍率（也就是 $Q_{in+b}$），將會完完整整地保留在狀態變數中，一路以高精度參與未來的回授迭代。直到濾波運算結束、準備將訊號寫入硬體輸出端（如 DAC）時，我們才會一次性地移去這段增益（即右移 S_{\text{out}} = Q_{in} + Q_b - Q_{out} 位元）。

3. 核心演算法步驟與數學推導

對於每個輸入採樣 $x_{\text{float}}[n]$，整數濾波器的計算流程如下：

步驟一：輸入量化 (Input Quantization)

將浮點數輸入訊號量化為整數 $x_{\text{int}}[n]$：

x_{\text{int}}[n] = \text{round}\left( x_{\text{float}}[n] \times 2^{Q_{in}} \right)

步驟二：前饋累加 (Feedforward Accumulation - 處理 b 係數)

計算前饋部分，累積至高精度累加器 sum_b（規格使用 32-bit 帶符號整數 int32_t）：

sum\_b = \sum_{i=0}^{N_b-1} b_{\text{int}}[i] \cdot x_{\text{int}}[n-i]

此時 $sum_b \in Q_{in+b}$。

步驟三：回授累加與歸一化 (Feedback Accumulation - 處理 a 係數)

  1. 計算高精度的回授乘加值：

    sum\_a = \sum_{j=1}^{N_a-1} a_{\text{int}}[j] \cdot y_{\text{hist}}[n-j]

    此時 $sum_a \in Q_{in+b+a}$。

  2. 歸一化 Q_a 放大倍數：將回授項除以 $A_0$（即 $a_{\text{int}}[0]$，通常為 $2^{Q_a}$）以對齊至 $Q_{in+b}$：

    • Floor 模式（無條件捨去，非常不建議）： sum\_a\_scaled = \lfloor \frac{sum\_a}{A_0} \rfloor \approx sum\_a \gg Q_a
    • Round 模式（四捨五入，專案預設與推薦）： round\_offset\_a = A_0 \gg 1 sum\_a\_scaled = \lfloor \frac{sum\_a + round\_offset\_a}{A_0} \rfloor \approx \left( sum\_a + (1 \ll (Q_a - 1)) \right) \gg Q_a

步驟四：更新狀態變數 (State Variable Update)

從前饋累加值中減去對齊後的回授值，並直接存入歷史狀態變數：

acc = sum\_b - sum\_a\_scaled y_{\text{hist}}[n] = acc

歷史狀態變數 y_{\text{hist}} 格式為 $Q_{in+b}$。

步驟五：最終輸出縮放 (Output Quantization & Scale-Out)

因為 y 的實際物理數值與系統增益（System Gain）主要由 b 和 x 決定，這時內部變數已經穩定保持在 Q_{in+b} 格式。 shift_out 的作用是將內部的高精度數值轉換為符合實體硬體接口（如 12-bit DAC）的位元寬度。 計算所需的總位移量 $S_{\text{out}}$：

S_{\text{out}} = Q_{in} + Q_b - Q_{out}
  • 當 $S_{\text{out}} > 0$（需要右移）：
    • Round 模式： round\_offset\_out = 1 \ll (S_{\text{out}} - 1) y_{\text{out}}[n] = (acc + round\_offset\_out) \gg S_{\text{out}}
  • 當 $S_{\text{out}} < 0$（需要左移）： y_{\text{out}}[n] = acc \ll (-S_{\text{out}})
  • 當 $S_{\text{out}} == 0$（不需位移）： y_{\text{out}}[n] = acc

4. 消除直流偏移的 1-Clock Rounding 演算法秘技

傳統的算術右移（Arithmetic Right Shift, SRA）在二進位中代表的是 Floor（向負無窮大取整）。 這會帶來平均 -0.5\text{ LSB} 的系統性截斷誤差，在時域訊號中會累積成顯著的直流偏移（DC Bias），在 IIR 濾波器中甚至會導致輸出不斷漂移。

在 C 語言中，千萬不要呼叫 round() 等浮點數函式，請使用以下純整數運算：

// 完美的 1-clock 四捨五入右移寫法：
// (1 << (shift - 1)) 在編譯時期會被編譯器直接優化為一個立即常數，
// 整體運算只比普通右移多消耗一個 ADD 指令，完美消除 DC Bias！
y_out = (acc + (1 << (shift - 1))) >> shift;

5. 完整 C 語言 (MCU 32-bit 最佳化) 實作程式碼範例

在 32-bit 微處理器（如 ARM Cortex-M 或是 RISC-V RV32I）上，int32_t 運算是效能最高、最節省時鐘週期的核心格式（避免了昂貴的 64-bit 軟體乘除模擬）。

#include <stdint.h>
#include <stdbool.h>

#define FILTER_ORDER_MAX 4  // 最大支援濾波器階數

typedef struct {
    int32_t b[FILTER_ORDER_MAX + 1];  // Q_b 格式的前饋係數
    int32_t a[FILTER_ORDER_MAX + 1];  // Q_a 格式的回授係數 (a[0] 通常為 1<<Q_a)
    int32_t nb;                       // b 係數個數
    int32_t na;                       // a 係數個數
    
    int16_t shift_in;   // Q_in
    int16_t shift_out;  // Q_out
    int16_t shift_b;    // Q_b
    int16_t shift_a;    // Q_a
    bool use_round;     // 是否啟用四捨五入補償
    
    // 歷史狀態變數 (循環緩衝區)
    int32_t x_hist[FILTER_ORDER_MAX + 1]; // Q_in 格式的輸入歷史
    int32_t y_hist[FILTER_ORDER_MAX + 1]; // Q_{in+b} 格式的高精度狀態歷史
    int32_t x_index;                      // x 歷史緩衝區指針
    int32_t y_index;                      // y 歷史緩衝區指針
} FixedFilter_t;

/**
 * @brief 初始化濾波器狀態與係數
 */
void FixedFilter_Init(FixedFilter_t *filter, 
                      const int32_t *b_coeffs, int32_t nb,
                      const int32_t *a_coeffs, int32_t na,
                      int16_t shift_in, int16_t shift_out,
                      int16_t shift_b, int16_t shift_a,
                      bool use_round) 
{
    filter->nb = (nb > FILTER_ORDER_MAX + 1) ? FILTER_ORDER_MAX + 1 : nb;
    filter->na = (na > FILTER_ORDER_MAX + 1) ? FILTER_ORDER_MAX + 1 : na;
    
    for (int32_t i = 0; i < filter->nb; i++) filter->b[i] = b_coeffs[i];
    for (int32_t j = 0; j < filter->na; j++) filter->a[j] = a_coeffs[j];
    
    filter->shift_in = shift_in;
    filter->shift_out = shift_out;
    filter->shift_b = shift_b;
    filter->shift_a = shift_a;
    filter->use_round = use_round;
    
    filter->x_index = 0;
    filter->y_index = 0;
    
    for (int32_t i = 0; i < FILTER_ORDER_MAX + 1; i++) {
        filter->x_hist[i] = 0;
        filter->y_hist[i] = 0;
    }
}

/**
 * @brief 核心定點數差分方程式執行 (32-bit 最佳化版本)
 * @param x_in Q_in 格式的整數輸入值
 * @return Q_out 格式的整數輸出值
 */
int32_t FixedFilter_Process(FixedFilter_t *filter, int32_t x_in)
{
    // 1. 將新輸入寫入循環緩衝區
    filter->x_hist[filter->x_index] = x_in;
    
    // 2. 計算前饋 (Feedforward): 處理 b 係數 (Q_b * Q_in -> Q_{in+b})
    int32_t sum_b = 0;
    int32_t x_ptr = filter->x_index;
    for (int32_t i = 0; i < filter->nb; i++) {
        sum_b += filter->b[i] * filter->x_hist[x_ptr];
        if (--x_ptr < 0) x_ptr = filter->nb - 1; // 循環指針回繞
    }
    
    // 3. 計算回授 (Feedback): 處理 a 係數 (Q_a * Q_{in+b} -> Q_{in+b+a})
    int32_t sum_a = 0;
    int32_t y_ptr = filter->y_index;
    for (int32_t j = 1; j < filter->na; j++) {
        // y_hist 儲存的是歷史中的高精度狀態變數 (Q_{in+b})
        sum_a += filter->a[j] * filter->y_hist[y_ptr];
        if (--y_ptr < 0) y_ptr = filter->na - 1;
    }
    
    // 4. 歸一化 a 係數的放大倍數 (除以 A0，即算術右移 shift_a)
    int32_t A0 = filter->a[0];
    if (A0 <= 0) A0 = 1; // 防呆
    
    int32_t sum_a_scaled = 0;
    if (filter->use_round) {
        int32_t round_offset_a = A0 >> 1;
        sum_a_scaled = (sum_a + round_offset_a) / A0;
    } else {
        sum_a_scaled = sum_a / A0;
    }
    
    // 如果 A0 是 2 的冪次方 (例如 2^shift_a)，上述除法可被編譯器優化為算術右移：
    // sum_a_scaled = (sum_a + (1 << (shift_a - 1))) >> shift_a;
    
    // 5. 更新狀態變數並寫入歷史緩衝區
    int32_t acc = sum_b - sum_a_scaled;
    
    if (++(filter->y_index) >= filter->na) filter->y_index = 0;
    filter->y_hist[filter->y_index] = acc;
    
    // 6. 計算最終輸出量化 Q_{in+b} -> Q_{out}
    int32_t out_shift = filter->shift_in + filter->shift_b - filter->shift_out;
    int32_t y_out = 0;
    
    if (out_shift > 0) {
        if (filter->use_round) {
            int32_t round_offset_out = (1 << (out_shift - 1));
            y_out = (acc + round_offset_out) >> out_shift;
        } else {
            y_out = acc >> out_shift;
        }
    } else if (out_shift < 0) {
        y_out = acc << (-out_shift);
    } else {
        y_out = acc;
    }
    
    if (++(filter->x_index) >= filter->nb) filter->x_index = 0;
    
    return y_out;
}

6. 與 Python 模擬的一致性驗證

本專案後端 dea/csv_processing.py 中的 integer_lfilter 函數提供完全一致的參考實作。工程師在撰寫 C 語言硬體代碼時，可隨時將硬體輸出的暫存器數值與本專案的 Python 數值進行對齊測試：

# 專案後端 Python 核心運算對照
def integer_lfilter(b_int, a_int, x_float, shift_in, shift_out, shift_b, shift_a, use_round=False):
    x_int = np.round(x_float * (2**shift_in)).astype(np.int64)
    y_hist = np.zeros(len(x_int), dtype=np.int64) # Q_{in+b} 狀態變數
    y_out = np.zeros(len(x_int), dtype=np.int64)  # Q_out 輸出變數
    
    out_shift = shift_in + shift_b - shift_out
    A0 = int(a_int[0])
    
    round_offset_a = (A0 >> 1) if use_round else 0
    round_offset_out = (1 << (out_shift - 1)) if (use_round and out_shift > 0) else 0

    for n in range(len(x_int)):
        sum_b = 0
        # 前饋累積處理 b 係數 (Q_{in+b})
        for i in range(len(b_int)):
            if n - i >= 0:
                sum_b += b_int[i] * x_int[n - i]
                
        sum_a = 0
        # 回授累積處理 a 係數 (Q_{in+b+a})
        for j in range(1, len(a_int)):
            if n - j >= 0:
                sum_a += a_int[j] * y_hist[n - j]
                
        # 歸一化 a 係數放大倍數，對齊回 Q_{in+b}
        sum_a_scaled = (sum_a + round_offset_a) // A0
        
        acc = sum_b - sum_a_scaled
        y_hist[n] = acc # 存入高精度狀態變數
        
        # 最終輸出量化至 Q_out
        if out_shift > 0:
            y_out[n] = (acc + round_offset_out) >> out_shift
        elif out_shift < 0:
            y_out[n] = acc << (-out_shift)
        else:
            y_out[n] = acc

    return y_out.astype(float) / (2**shift_out)

7. 核心探討：高階直接型 (Direct Form) vs. 二階級聯 (Biquad / SOS) 的黃金法則

在實際將設計好的濾波器移植到定點數 MCU、DSP 或 FPGA 時，開發者常面臨一個核心抉擇：當需要高階濾波器時，應該將多個低階濾波器的係數做多項式相乘（摺積）合併為單一高階濾波器（直接型，Direct Form），還是將其拆解為多個獨立的二階濾波器進行串聯級聯（Cascade of Biquads / SOS）？

以下從數學限制、數值穩定度、計算效率等三個維度進行深度剖析，並說明為何**「二階級聯」是實務上無可爭議的黃金標準**。

7.1 數學限制：為什麼不能只用一階濾波器級聯？

一階 IIR 濾波器的轉移函數為：

H(z) = \frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1}}
  • 實數極點限制：在實體系統中，濾波器係數 a_1, b_1 必須是實數 (Real Numbers)。這代表一階濾波器的極點（$z = -a_1$）與零點（$z = -b_0/b_1$）只能落在複數平面的實數軸上。
  • 無法實現共軛複數極點：現代陡峭濾波器（如 Butterworth, Chebyshev）或窄頻共振/陷波濾波器（High-Q Notch Filter）必須依賴共軛複數極點對 (Complex Conjugate Pole Pairs) 來提供系統共振，進而撐起頻率響應的陡峭邊緣。
  • 二階為最小基底：要產生一對共軛複數極點 p_1, p_2 = r e^{\pm j\theta} 且維持係數為實數，分母多項式至少必須是二階： (1 - p_1 z^{-1})(1 - p_2 z^{-1}) = 1 - 2r\cos(\theta) z^{-1} + r^2 z^{-2} 因此，二階濾波器 (Biquad) 是容納共軛複數極點並保持實數係數的最小數學單元。若只用一階濾波器級聯，將永遠無法設計出高選擇性的濾波器。

7.2 數值穩定度：威爾金森多項式效應 (Wilkinson's Phenomenon)

當我們硬將多個二階多項式相乘合併為一個高階直接型濾波器時，在定點數系統中會面臨致命的不穩定 RISK：

  1. 極點對係數的靈敏度指數級暴增： 多項式根（極點）對係數微小擾動的偏差量 \delta z_k 與該極點到其他所有極點的距離乘積成反比： \delta z_k \propto \frac{1}{\prod_{i \neq k} (z_k - z_i)} 高階濾波器的極點通常在 Z 平面上靠得非常近（尤其是在低頻或高 Q 值設計中）。極點之間的極小距離使得分母極小，這意味著即便定點數係數因為捨入誤差（Quantization Error）只改變了最後一個 LSB，都會導致極點發生巨大的位移，甚至直接飛出單位圓外 (|z| > 1) 導致濾波器發散崩潰。
  2. 係數動態範圍懸殊： 高階多項式合併後，各項係數值往往會產生好幾個數量級的差距（例如部分中間項係數為 $3.8$，而有些高階項僅為 $0.002$）。 為了配合同一個定點數 Q 格式表示，極小的係數會面臨嚴重下溢 (Underflow) 被直接抹零，導致頻率響應產生劇烈畸變。

在二階級聯中，每個 Stage 僅有 2 個極點，距離較遠，極點靈敏度極低，且每級係數大小相近，能將 Q 格式 the 二進位精度發揮到極致。

7.3 計算效率與資源開銷對比

直覺上，合併為高階直接型似乎可以省下很多計算量，但精確計算指令次數後會發現這個優勢極其微弱。 以 4 階濾波器 (Order N=4) 為例，對比兩個 2 階濾波器 (Biquads) 串聯（假設 $a_0 = 1$）：

指標 4 階直接型 (Single 4th-Order) 兩個 2 階級聯 (Cascaded 2 x Biquad) 差異與分析
係數個數 b_0 \sim b_4, a_1 \sim a_4 (共 9 個) 每級 5 個 (共 10 個) 直接型少 1 個係數
每樣本乘法次數 9 次 10 次 直接型僅節省 1 次乘法 (約 10%)
每樣本加法次數 8 次 8 次 完全相同
暫存器/狀態變數 4 個歷史狀態 (DF-II) 4 個歷史狀態 (2 x DF-II) 完全相同
  • 結論：合併成一個 4 階濾波器，僅僅只能省下 1 次乘法！在具備硬體乘加單元（MAC）的現代處理器上，這種效能差異微不足道，根本無法彌補其喪失數值穩定度的巨大代價。

7.4 實務設計與移植啟發 (Engineering Guidelines)

  • SOS 標準化：不論濾波器階數多高，工程實務上皆應使用 二階區域 (Second-Order Sections, SOS) 級聯實作，並在晶片端撰寫單一高度優化的 Biquad 核心處理單元以復用程式碼。
  • 級聯配對與排列最佳化：
    • 極零點配對 (Pole-Zero Pairing)：將最靠近的極點與零點配在同一個二階級聯 Stage，以防中間頻率增益峰值過高導致暫存器飽和溢位。
    • 順序排列 (Ordering)：將極點距離單位圓較遠（Q 值較低、響應較平緩）的 Stage 排在前面；將靠近單位圓（Q 值較高、共振較強）的 Stage 排在後面，以獲得最佳的動態範圍與信噪比 (SNR)。
Reference in New Issue View Git Blame Copy Permalink