2

我正在尝试按照本书中的说明编写一个简单的带通滤波器。我的代码创建了一个布莱克曼窗口，并结合了两个低通滤波器内核以使用频谱反转创建一个带通滤波器内核，如此处的第二个示例中所述（表 16-2）。

我正在通过将代码与我在 matlab 中获得的结果进行比较来测试我的代码。当我分别测试创建布莱克曼窗口和低通滤波器内核的方法时，我得到的结果接近于我在 matlab 中看到的结果（最多小数点后的一些数字 - 我将错误归因于 java 双变量舍入问题)，但我的带通滤波器内核不正确。

我运行的测试：

  • 创建了一个 blackman 窗口并将其与我在 matlab 中得到的进行比较 - 一切都很好。
  • fir1(N, Fc1/(Fs/2), win, flag);使用我的代码和matlab使用此窗口创建了一个低通滤波器（请参阅下面的完整代码）。我认为结果是正确的，尽管 Fc1 越大，我的错误就越大（为什么？）
  • 使用我的代码和在 matlab 中创建了一个泛通滤波器fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);- 结果完全关闭。
  • 使用我的代码和 matlab 生成的内核过滤我的数据——一切都很好。

那么 - 为什么我的带通滤波器内核关闭了？我做错什么了？我想我要么有错误，要么 fir1 使用了不同的算法，但我无法检查，因为其文档中引用的文章未公开。

这是我的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 代码（我相信 bug 在 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);
4

1 回答 1

1

我最终在 Java 中实现了 Matlab 的 fir1 函数。我的结果非常准确。

于 2013-05-10T22:24:42.430 回答