FFTW 頻域濾波實作：C++ 代碼除以 N 導致波形異常的根因分析

問題現象

使用 FFTW 執行 R2C 正向轉換、頻域濾波、C2R 逆向轉換後，輸出波形與原始資料極其相似，而非預期的頻帶濾波結果。MATLAB 實作同邏輯則正確。

核心程式碼片段

// 正向 FFT
fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out1, FFTW_ESTIMATE);
fftw_execute(p);

// 頻域濾波
for (int m = 0; m < Nout-1; m++) {
    float K = m / (Nout / nSamples); // 問題行
    if (K < f1 || K > f2) { out2[m][0] = out2[m][1] = 0; }
    else { out2[m][0] = out1[m][0]; out2[m][1] = out1[m][1]; }
}

// 逆向 FFT
fftw_plan c2cP = fftw_plan_dft_c2r_1d(N, out2, out3, FFTW_ESTIMATE);
fftw_execute(c2cP);
// FFTW 逆變換不自動除以 N，需手動縮放
for (int i = 0; i < N; i++) out3[i] /= N;

根因定位

關鍵錯誤在於頻率索引計算：

float K = m / (Nout / nSamples); // 整數除法導致 K 恆為 0 或 1

C/C++ 中整數除法會捨去小數部分，導致 `K` 無法正確映射至實際頻率，濾波條件失效，幾乎所有頻率分量被保留，逆變換後自然接近原始波形。

修正方案

float K = 1.0f * m / (Nout / nSamples); // 強制浮點除法
// 或
float K = (float)m / (Nout / nSamples);

FFTW 正確使用要點補充

1. 記憶體對齊：使用 `fftw_malloc` 分配輸入/輸出陣列，確保 SIMD 優化生效
2. Plan 重用：相同尺寸多次轉換時重用 `fftw_plan`，避免重複規劃開銷
3. 標準化縮放：FFTW 正向+逆向未自動歸一化，必須手動除以 `N`（如代碼所示）
4. Hermitian 對稱：R2C 輸出僅含非負頻率，濾波時需同時處理對稱的負頻率部分（若使用 C2C 則需完整頻譜）
5. 頻率軸映射：第 `k` 個輸出對應頻率 `k * fs / N`（fs 為採樣率），而非簡單索引比例

正確頻率映射公式

float fs = nSamples; // 假設採樣率 = 樣本數/秒
for (int m = 0; m < Nout; m++) {
    float freq = m * fs / N; // 真實頻率 Hz
    if (freq < fmin || freq > fmax) { /* 零化 */ }
}

驗證建議

修正後，輸出波形應明顯不同於原始波形，僅保留目標頻帶成分。可對比 MATLAB `fft` -> 零化 -> `ifft` 結果確認一致性。