使用java实现数字滤波器

 

本文系转载，原文地址：https://blog.csdn.net/weijifen000/article/details/80347973

 

使用java实现数字滤波器

2018年05月17日 11:06:21 weijifen000 阅读数 1813

借助 dsp-collection 库实现 java 数字滤波器。

在网上搜索“使用 java 实现数字滤波器”，搜索结果大多是利用 MATLAB 生成滤波器系数 a 、 b,然后在 java 中进行滤波，但是我的需求是整个过程完全在 java 中完成。之后搜索“ java 实现 DSP Filter ”则在stackoverflow中找到了解决办法。

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]);
}

• 1
• 2
• 3
• 4
• 5
• 6

• 1
• 2
• 3
• 4
• 5
• 6

生成滤波器：

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

• 1
• 2
• 3

• 1
• 2
• 3

查看系数的值：

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]);
        }
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13

输出为：

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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22

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

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

• 1

• 1

IIR 系统的表达式：
y(n)+∑N1a(k)y(n−k)=∑M0b(r)x(n−r)

y(n)+1∑N​a(k)y(n−k)=0∑M​b(r)x(n−r)
封装的滤波器函数,里面的数据类型可以改为 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
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34

滤波前后的幅频图像（左图为滤波前的幅频图，右图为滤波后的幅频图）：

滤波之后的幅频图像：

如何对信号进行频谱分析以及将频谱图在 Android 中显示请参见：
java 实现离散时间傅里叶变换

注意：

1. 滤波器系数：

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

• 1
• 2

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

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

• 1
• 2

• 1
• 2

结果：

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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22

可以看到，滤波器系数和之前是一样的。
当使用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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8

输出为：

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]

• 1
• 2
• 3
• 4
• 5
• 6
• 7

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

参考资料：

1. IIR 滤波器设计（ java ） https://blog.csdn.net/believe646875970/article/details/53909318
2. https://stackoverflow.com/questions/636686/signal-processing-library-in-java?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa