藉助 dsp-collection 庫實現 java 數字濾波器。

在網上搜索“使用 java 實現數字濾波器”，搜尋結果大多是利用 MATLAB 生成濾波器係數 a 、 b,然後在 java 中進行濾波，但是我的需求是整個過程完全在 java 中完成。之後搜尋“ java 實現 DSP Filter ”則在stackoverflow中找到了解決辦法。

dsp-collection 庫 支援多種 IIR 濾波器，包括巴特沃斯，切比雪夫和貝塞爾濾波器。筆者在這裡使用的是巴特沃斯濾波器。其他的請參考 API 文件 。

模擬訊號:
$f(t)=sin(10 \pi t)+sin(30 \pi t)$

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 系統的表示式：
$y(n)+\sum _1^N a(k)y(n-k)=\sum _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;
    }

濾波前後的幅頻影象（左圖為濾波前的幅頻圖，右圖為濾波後的幅頻圖）：

如何對訊號進行頻譜分析以及將頻譜圖在 Android 中顯示請參見：
java 實現離散時間傅立葉變換

注意：

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 生成濾波器係數的是一致的。

參考資料：