diff --git a/INTEGER_FILTER_IMPLEMENTATION.md b/INTEGER_FILTER_IMPLEMENTATION.md new file mode 100644 index 0000000..98a83fc --- /dev/null +++ b/INTEGER_FILTER_IMPLEMENTATION.md @@ -0,0 +1,293 @@ +# 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()` 等浮點數函式，請使用以下純整數運算： + +```c +// 完美的 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 軟體乘除模擬）。 + +```c +#include +#include + +#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<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 +# 專案後端 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) +``` diff --git a/dea/csv_processing.py b/dea/csv_processing.py index bea7aec..5849964 100755 --- a/dea/csv_processing.py +++ b/dea/csv_processing.py @@ -82,23 +82,23 @@ def integer_lfilter(b_int, a_int, x_float, shift_in, shift_out, shift_b, shift_a round_offset_out = (1 << (out_shift - 1)) if (use_round and out_shift > 0) else 0 for n in range(len(x_int)): - acc = 0 + sum_b = 0 # Feedforward: 前饋完全不位移，結果為 Q_{in + b} for i in range(nb): if n - i >= 0: - acc += b_int[i] * x_int[n - i] + sum_b += b_int[i] * x_int[n - i] # Feedback: 歷史紀錄為 Q_{in + b}，係數為 Q_a，乘積為 Q_{in + b + a} - fb_sum = 0 + sum_a = 0 for j in range(1, na): if n - j >= 0: - fb_sum += a_int[j] * y_hist[n - j] + sum_a += a_int[j] * y_hist[n - j] - # Feedback 縮放：除以 A0 (即 >> shift_a)，結果退回 Q_{in + b} 完美對齊 + # Feedback 縮放：前提 A0 = 1 (或者代表放大了 Q_a 倍的常數)，除以 A0 (即 >> shift_a) 歸一化 # Python 的 // 等同於硬體的 SRA (Arithmetic Right Shift)，會向負無窮大 Floor - fb_shifted = (fb_sum + round_offset_a) // A0 + sum_a_scaled = (sum_a + round_offset_a) // A0 - acc -= fb_shifted + acc = sum_b - sum_a_scaled # 將超高精度的 acc 直接存入歷史變數 y_hist[n] = acc diff --git a/vite.config.js b/vite.config.js index c72015c..2182c61 100644 --- a/vite.config.js +++ b/vite.config.js @@ -7,8 +7,9 @@ export default defineConfig({ server: { proxy: { '/api': { - target: 'http://127.0.0.1:8000', + target: 'https://127.0.0.1:8000', changeOrigin: true, + secure: false, }, }, },