Matlab filtfilt（）函数在Java中的实现

Question

有没有人尝试过filtfilt()用Java （或至少在C ++中）实现matlab的功能？如果你们有一个算法，那将有很大的帮助。

Answer 1

这是我filtfilt在MATLAB中实现的算法的C ++ 实现。希望这对您有所帮助。

Answer 2

好吧，我知道这个问题很古老，但也许我可以帮助到这里想知道filtfilt实际情况的其他人。

尽管从前filtfilt向后向（又名零相位）过滤的文档中可以明显看出，它如何处理填充和初始条件之类的事情对我来说并不那么明显。

由于我在此处（或其他地方）找不到任何其他答案以及有关 的这些实现细节的足够信息filtfilt，因此我基于其源代码和文档（因此，不是，也不是，但是）实现了Python's的简化版本。我相信该版本的工作方式与's相同。scipy.signal.filtfiltJavaC++PythonscipyMatlab

为了简单起见，下面的代码是为二阶专门编写IIR滤波器，并且它假定系数矢量a和b是已知的（例如，从获得的scipy.signal.butter，或通过手算）。

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

免责声明：此代码仅旨在提供对 的某些实现细节的一些见解filtfilt，因此目标是清晰而不是计算效率/性能。该scipy.signal.filtfilt执行的速度要快得多（如100X根据一个快速和肮脏的更快timeit我的系统上测试）。

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。此外，numpy通过写出解决方案zi并使用列表而不是 numpy 数组，它可以很容易地适应在纯 python 中工作，甚至没有。这甚至带来了显着的性能优势，因为在 python 循环中访问单个 numpy 数组元素比访问列表元素慢得多。

过滤器本身是在一个简单的Python循环中实现的。它使用状态空间表示，因为它无论如何都用于确定初始条件（请参阅 参考资料scipy.signal.lfilter_zi）。我相信scipy线性过滤器的实际实现（即scipy.signal.sigtools._linear_filter）在 中做了类似的事情C，可以在这里看到（感谢这个答案）。

下面是一些代码，提供了一个（非常基本的）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 位小数：