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函数。我的结果非常准确。