java - Matlab filtfilt()函数在Java中的实现

Question

有人试过用 Java (或至少用 C++) 实现 matlab 的 filtfilt() 函数吗？如果你们有算法，那将有很大帮助。

Answer 1

好吧，我知道这个问题很古老，但也许我可以帮助其他想知道什么的人 filtfilt确实如此。

尽管从文档中可以明显看出 filtfilt进行前向-后向(也称为零相位)过滤，我不太清楚它如何处理填充和初始条件之类的事情。

因为我无法在这里(或其他地方)找到任何其他答案以及关于 filtfilt 的这些实现细节的足够信息。 , 我实现了 Python 的简化版本的 scipy.signal.filtfilt ，基于其来源和文档(因此，不是 Java ，也不是 C++ ，而是 Python )。我相信scipy版本 works the same方式Matlab的。

为了简单起见，下面的代码是专门为第二个订单编写的IIR filter , 它假设系数 vector a和 b是已知的(例如从 scipy.signal.butter 或 calculated by hand 获得)。

它匹配 filtfilt默认行为，使用 odd填充长度 3 * max(len(a), len(b)) ，在前向传播之前应用。使用来自 scipy.signal.lfilter_zi 的方法找到初始状态(docs)。

免责声明:此代码仅旨在提供对 filtfilt 的某些实现细节的一些见解。 ，因此目标是清晰度而不是计算效率/性能。 scipy.signal.filtfilt实现速度要快得多(例如，根据在我的系统上进行的快速而肮脏的 timeit 测试，速度提高了 100 倍)。

import numpy


def custom_filter(b, a, x):
    """ 
    Filter implemented using state-space representation.

    Assume a filter with second order difference equation (assuming a[0]=1):

        y[n] = b[0]*x[n] + b[1]*x[n-1] + b[2]*x[n-2] + ...
                         - a[1]*y[n-1] - a[2]*y[n-2]

    """
    # State space representation (transposed direct form II)
    A = numpy.array([[-a[1], 1], [-a[2], 0]])
    B = numpy.array([b[1] - b[0] * a[1], b[2] - b[0] * a[2]])
    C = numpy.array([1.0, 0.0])
    D = b[0]

    # Determine initial state (solve zi = A*zi + B, see scipy.signal.lfilter_zi)
    zi = numpy.linalg.solve(numpy.eye(2) - A, B)

    # Scale the initial state vector zi by the first input value
    z = zi * x[0]

    # Apply filter
    y = numpy.zeros(numpy.shape(x))
    for n in range(len(x)):
        # Determine n-th output value (note this simplifies to y[n] = z[0] + b[0]*x[n])
        y[n] = numpy.dot(C, z) + D * x[n]
        # Determine next state (i.e. z[n+1])
        z = numpy.dot(A, z) + B * x[n]
    return y

def custom_filtfilt(b, a, x):
    # Apply 'odd' padding to input signal
    padding_length = 3 * max(len(a), len(b))  # the scipy.signal.filtfilt default
    x_forward = numpy.concatenate((
        [2 * x[0] - xi for xi in x[padding_length:0:-1]],
        x,
        [2 * x[-1] - xi for xi in x[-2:-padding_length-2:-1]]))

    # Filter forward
    y_forward = custom_filter(b, a, x_forward)

    # Filter backward
    x_backward = y_forward[::-1]  # reverse
    y_backward = custom_filter(b, a, x_backward)

    # Remove padding and reverse
    return y_backward[-padding_length-1:padding_length-1:-1]

请注意，此实现不需要scipy .此外，它可以很容易地适应纯 python，甚至不需要 numpy。 , 通过写出 zi 的解决方案并使用列表而不是 numpy 数组。这甚至带来了显着的性能优势，因为在 python 循环中访问单个 numpy 数组元素比访问列表元素慢得多。

过滤器本身是在一个简单的 Python 中实现的环形。它使用状态空间表示，因为它无论如何都用于确定初始条件(参见 scipy.signal.lfilter_zi )。我相信实际的scipy线性滤波器的实现(即 scipy.signal.sigtools._linear_filter )在 C 中做了类似的事情，可以看出here (感谢 this answer)。

这里有一些代码提供(非常基本的)检查 scipy 是否相等。输出和custom输出:

import numpy
import numpy.testing
import scipy.signal
from matplotlib import pyplot
from . import custom_filtfilt

def sinusoid(sampling_frequency_Hz=50.0, signal_frequency_Hz=1.0, periods=1.0,
             amplitude=1.0, offset=0.0, phase_deg=0.0, noise_std=0.1):
    """
    Create a noisy test signal sampled from a sinusoid (time series)

    """
    signal_frequency_rad_per_s = signal_frequency_Hz * 2 * numpy.pi
    phase_rad = numpy.radians(phase_deg)
    duration_s = periods / signal_frequency_Hz
    number_of_samples = int(duration_s * sampling_frequency_Hz)
    time_s = (numpy.array(range(number_of_samples), float) /
              sampling_frequency_Hz)
    angle_rad = signal_frequency_rad_per_s * time_s
    signal = offset + amplitude * numpy.sin(angle_rad - phase_rad)
    noise = numpy.random.normal(loc=0.0, scale=noise_std, size=signal.shape)
    return signal + noise


if __name__ == '__main__':
    # Design filter
    sampling_freq_hz = 50.0
    cutoff_freq_hz = 2.5
    order = 2
    normalized_frequency = cutoff_freq_hz * 2 / sampling_freq_hz
    b, a = scipy.signal.butter(order, normalized_frequency, btype='lowpass')

    # Create test signal
    signal = sinusoid(sampling_frequency_Hz=sampling_freq_hz,
                      signal_frequency_Hz=1.5, periods=3, amplitude=2.0,
                      offset=2.0, phase_deg=25)

    # Apply zero-phase filters
    filtered_custom = custom_filtfilt(b, a, signal)
    filtered_scipy = scipy.signal.filtfilt(b, a, signal)

    # Verify near-equality
    numpy.testing.assert_array_almost_equal(filtered_custom, filtered_scipy,
                                            decimal=12)

    # Plot result
    pyplot.subplot(1, 2, 1)
    pyplot.plot(signal)
    pyplot.plot(filtered_scipy)
    pyplot.plot(filtered_custom, '.')
    pyplot.title('raw vs filtered signals')
    pyplot.legend(['raw', 'scipy filtfilt', 'custom filtfilt'])
    pyplot.subplot(1, 2, 2)
    pyplot.plot(filtered_scipy-filtered_custom)
    pyplot.title('difference (scipy vs custom)')
    pyplot.show()

这个基本比较产生如下图，表明至少 14 位小数相等，对于这种特定情况(我猜是机器精度？):