java中的簡單帶通濾波器

Question

我正試圖按照this book中的說明書寫一個簡單的帶通濾波器。我的代碼創建了一個blackman窗口，並將兩個低通濾波器內核組合起來，使用頻譜反轉創建一個帶通濾波器內核，如第二個示例here（表16-2）中所述。

我正在測試我的代碼，將它與我在matlab中得到的結果進行比較。當我測試分別創建blackman窗口和低通濾波器內核的方法時，我得到的結果與我在matlab中看到的接近（在小數點後最多有一些數字 - 我將錯誤歸因於java double變量舍入問題），但我的帶通濾波器內核不正確。

測試我跑：

• 創建一個布萊克曼窗，並比較其與我在MATLAB獲得 - 都很好。
• 創建一個低通濾波器使用這個窗口使用我的代碼和fir1(N, Fc1/(Fs/2), win, flag);在MATLAB（見下面的完整代碼）。我認爲結果是正確的，雖然我得到更大的誤差更大的Fc1是（爲什麼？）
• 在matlab中創建了一個使用我的代碼和fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);的pand pass filter - 結果完全關閉。
• 使用我的代碼和matlab生成的內核過濾我的數據 - 都很好。

所以 - 爲什麼我的帶通濾波器內核關閉？我做錯了什麼？ 我想我有一個bug或fir1使用不同的算法，但我無法檢查，因爲article referenced in its documentation不公開。

這是我的MATLAB代碼：

Fs = 200;  % Sampling Frequency 
N = 10;   % Order 
Fc1 = 1.5;   % First Cutoff Frequency 
Fc2 = 7.5;   % Second Cutoff Frequency 
flag = 'scale';  % Sampling Flag 

% Create the window vector for the design algorithm. 
win = blackman(N+1); 

% Calculate the coefficients using the FIR1 function. 
b = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag); 
Hd = dfilt.dffir(b); 
res = filter(Hd, data); 

這是我的Java代碼（我相信錯誤是在bandPassKernel）：

/** 
    * See - http://www.mathworks.com/help/signal/ref/blackman.html 
    * @param length 
    * @return 
    */ 
    private static double[] blackmanWindow(int length) { 

     double[] window = new double[length]; 
     double factor = Math.PI/(length - 1); 

     for (int i = 0; i < window.length; ++i) { 
      window[i] = 0.42d - (0.5d * Math.cos(2 * factor * i)) + (0.08d * Math.cos(4 * factor * i)); 
     } 

     return window; 
    } 

private static double[] lowPassKernel(int length, double cutoffFreq, double[] window) { 

    double[] ker = new double[length + 1]; 
    double factor = Math.PI * cutoffFreq * 2; 
    double sum = 0; 

    for (int i = 0; i < ker.length; i++) { 
     double d = i - length/2; 
     if (d == 0) ker[i] = factor; 
     else ker[i] = Math.sin(factor * d)/d; 
     ker[i] *= window[i]; 
     sum += ker[i]; 
    } 

    // Normalize the kernel 
    for (int i = 0; i < ker.length; ++i) { 
     ker[i] /= sum; 
    } 

    return ker; 
} 

private static double[] bandPassKernel(int length, double lowFreq, double highFreq) { 

    double[] ker = new double[length + 1]; 
    double[] window = blackmanWindow(length + 1); 

    // Create a band reject filter kernel using a high pass and a low pass filter kernel 
    double[] lowPass = lowPassKernel(length, lowFreq, window); 

    // Create a high pass kernel for the high frequency 
    // by inverting a low pass kernel 
    double[] highPass = lowPassKernel(length, highFreq, window); 
    for (int i = 0; i < highPass.length; ++i) highPass[i] = -highPass[i]; 
    highPass[length/2] += 1; 

    // Combine the filters and invert to create a bandpass filter kernel 
    for (int i = 0; i < ker.length; ++i) ker[i] = -(lowPass[i] + highPass[i]); 
    ker[length/2] += 1; 

    return ker; 
} 

private static double[] filter(double[] signal, double[] kernel) { 

    double[] res = new double[signal.length]; 

    for (int r = 0; r < res.length; ++r) { 

     int M = Math.min(kernel.length, r + 1); 
     for (int k = 0; k < M; ++k) { 
      res[r] += kernel[k] * signal[r - k]; 
     } 
    } 

    return res; 
} 

這就是我如何使用我的代碼：

double[] kernel = bandPassKernel(10, 1.5d/(200/2), 7.5d/(200/2)); 
double[] res = filter(data, kernel); 

Answer 1

我最終在Java中實現了Matlab的fir1函數。我的結果非常準確。