pchang718
14 KiB
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），本架構採用前饋完全不右移的策略：
-
狀態變數
y_{\text{hist}}精度保持：歷史狀態變數y_{\text{hist}}[n-j]並非儲存低精度的最終輸出 $y_{\text{out}}$，而是直接儲存前饋與回授累加後、尚未進行輸出位移的高精度狀態值。其格式固定為：y_{\text{hist}} \in Q_{in + b} -
免除中間位移：在前饋運算中，
b_{\text{int}} \times x_{\text{int}}的結果為 $Q_b \times Q_{in} = Q_{in+b}$，完全不需進行任何位移，即可直接累加至sum_b中。 -
回授對齊縮放（歸一化 $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倍）在運算過程中是否也需要被歸一化？
答案是：完全不需要，而且我們也故意沒有做！這是因為：
- 避免資訊截斷（Truncation Noise）：如果我們在累加到
sum_b時就進行右移歸一，會直接把低位元的小數精度給丟棄，造成無法挽回的量化雜訊。 - 最優化運算效能：在迴圈內部，
sum_b（$Q_{in+b}$）與被歸一後的sum_a_scaled（$Q_{in+b}$）可以直接進行減法運算。整個迴圈內部不需要對x與b進行 any 移位，這在硬體實作中極大提升了運算速度。 - 延遲至最終輸出對齊：輸入與前饋所累積下來的放大倍率（也就是 $Q_{in+b}$），將會完完整整地保留在狀態變數中，一路以高精度參與未來的回授迭代。直到濾波運算結束、準備將訊號寫入硬體輸出端（如 DAC）時，我們才會一次性地移去這段增益（即右移
S_{\text{out}} = Q_{in} + Q_b - Q_{out}位元）。
- 避免資訊截斷（Truncation Noise）：如果我們在累加到
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 係數)
-
計算高精度的回授乘加值：
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}$。
-
歸一化
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 1sum\_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
- Floor 模式（無條件捨去，非常不建議）：
步驟四：更新狀態變數 (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}}
- Round 模式：
- 當 $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)