使用java实现数字滤波器

dsp-collection 库 支持多种 IIR 滤波器，包括巴特沃斯，切比雪夫和贝塞尔滤波器。笔者在这里使用的是巴特沃斯滤波器。其他的请参考 API 文档 。

仿真信号: 

f(t)=sin(10πt)+sin(30πt)f(t)=sin(10πt)+sin(30πt)

其中包含 5Hz 和 15Hz 两个频率分量,信号的采样率为 50Hz ：

double[] time = new double[150];
double[] valueA = new double[150];
for (int i = 0; i < 50 * 3; i++) {
    time[i] = i / 50.0;
    valueA[i] = Math.sin(2 * Math.PI * 5 * time[i])+Math.sin(2 * Math.PI * 15 * time[i]);
}

生成滤波器：

IirFilterCoefficients iirFilterCoefficients;
iirFilterCoefficients = IirFilterDesignExstrom.design(FilterPassType.lowpass, 10,
        10.0 / 50, 10.0 / 50);

查看系数的值：

public class Filter {
    public static void main(String[] args) {
        IirFilterCoefficients iirFilterCoefficients;
        iirFilterCoefficients = IirFilterDesignExstrom.design(FilterPassType.lowpass, 10,
                10.0 / 50, 10.0 / 50);
        for (int i=0;i<iirFilterCoefficients.a.length;i++) {
            System.out.println("A["+i+"]:"+iirFilterCoefficients.a[i]);
        }
        for (int i=0;i<iirFilterCoefficients.b.length;i++) {
            System.out.println("B["+i+"]:"+iirFilterCoefficients.b[i]);
        }
    }
}

输出为：

A[0]:1.0
A[1]:-1.9924014816014122
A[2]:3.0194828633553836
A[3]:-2.8185224264945132
A[4]:2.038720637062525
A[5]:-1.0545446210956793
A[6]:0.41444626875039875
A[7]:-0.1157186252368282
A[8]:0.022498509272218276
A[9]:-0.0026689123535761005
A[10]:1.4876445217776215E-4
B[0]:4.994540782331007E-4
B[1]:0.004994540782331008
B[2]:0.022475433520489533
B[3]:0.05993448938797209
B[4]:0.10488535642895115
B[5]:0.12586242771474138
B[6]:0.10488535642895115
B[7]:0.05993448938797209
B[8]:0.022475433520489533
B[9]:0.004994540782331008
B[10]:4.994540782331007E-4

使用巴特沃斯滤波器，滤波器设置为低通，滤波器阶数为 10 ，下截止频率为 10Hz 。则 “iirFilterCoefficients.a“ 为滤波器系数 a,“iirFilterCoefficients.b“ 为滤波器系数 b 。 使用滤波器处理信号：

valueA = IIRFilter(valueA, iirFilterCoefficients.a, iirFilterCoefficients.b);

IIR 系统的表达式：

封装的滤波器函数,里面的数据类型可以改为 float ：

public synchronized double[] IIRFilter(double[] signal, double[] a, double[] b) {

        double[] in = new double[b.length];
        double[] out = new double[a.length-1];

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

        for (int i = 0; i < signal.length; i++) {

            System.arraycopy(in, 0, in, 1, in.length - 1);
            in[0] = signal[i];

            //calculate y based on a and b coefficients
            //and in and out.
            float y = 0;
            for(int j = 0 ; j < b.length ; j++){
                y += b[j] * in[j];

            }

            for(int j = 0;j < a.length-1;j++){
                y -= a[j+1] * out[j];
            }

            //shift the out array
            System.arraycopy(out, 0, out, 1, out.length - 1);
            out[0] = y;

            outData[i] = y;


        }
        return outData;
    }

注意： 
1. 滤波器系数：

iirFilterCoefficients = IirFilterDesignExstrom.design(FilterPassType.lowpass, 10,
        10.0 / 50, 10.0 / 50);

··design()“ 方法的四个参数，第一个表示选择滤波器类型（低通，高通，带通，带阻），第二个参数表示滤波器的阶数，第三个参数表示下截止频率，第四个参数表示上截止频率，在低通和高通滤波器中只需要一个截止频率，所以在低通和高通中上截止频率是没有作用的。只有在带通或者带阻滤波器中才需要两个截止频率。 
例如当设置为：

iirFilterCoefficients = IirFilterDesignExstrom.design(FilterPassType.lowpass, 10,
        10.0 / 50, 13.0 / 50);

结果：

A[0]:1.0
A[1]:-1.9924014816014122
A[2]:3.0194828633553836
A[3]:-2.8185224264945132
A[4]:2.038720637062525
A[5]:-1.0545446210956793
A[6]:0.41444626875039875
A[7]:-0.1157186252368282
A[8]:0.022498509272218276
A[9]:-0.0026689123535761005
A[10]:1.4876445217776215E-4
B[0]:4.994540782331007E-4
B[1]:0.004994540782331008
B[2]:0.022475433520489533
B[3]:0.05993448938797209
B[4]:0.10488535642895115
B[5]:0.12586242771474138
B[6]:0.10488535642895115
B[7]:0.05993448938797209
B[8]:0.022475433520489533
B[9]:0.004994540782331008
B[10]:4.994540782331007E-4

可以看到，滤波器系数和之前是一样的。 
当使用python的时候：

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
b, a = signal.butter(10, 0.4, 'low')
print "A:"
print a
print "B:"
print b

输出为：

A:
[  1.00000000e+00  -1.99240148e+00   3.01948286e+00  -2.81852243e+00
   2.03872064e+00  -1.05454462e+00   4.14446269e-01  -1.15718625e-01
   2.24985093e-02  -2.66891235e-03   1.48764452e-04]
B:
[ 0.00049945  0.00499454  0.02247543  0.05993449  0.10488536  0.12586243
  0.10488536  0.05993449  0.02247543  0.00499454  0.00049945]

可以看到，与 java 生成滤波器系数的是一致的。