与 numpy.filtfilt 相当的 C# 数字滤波器实现

Question

所以，我尝试在 C# 中实现 scipy.signal.filtfilt 函数，但是我没有得到我在 python 中得到的结果，要么我看到了大的瞬态，要么看到了非常弱的信号。我查看了 scipy 函数的来源并试图模仿它，但没有成功。

我直接从我的python脚本中获取了过滤器系数a和b以及初始状态，所以没有区别。与 filtfilt 中的方法类似，我通过将信号旋转 180° 将两端的数据填充为滤波器大小的 3 倍（也在这里完成MATLAB's filtfilt() Algorithm）。初始状态乘以填充数据的第一个元素。 过滤器本身计算新输出和新输入数据的所有状态 d[]。

那么，有人知道我的错误在哪里吗？

这是我的代码：

/// <remarks>
/// The filter function is implemented as a direct II transposed structure. This means that the filter implements:
///
///a[0]*y[n] = b[0]*x[n] + b[1]*x[n - 1] + ... + b[M]*x[n - M]
///              - a[1]*y[n - 1] - ... - a[N]*y[n - N]
///where M is the degree of the numerator, N is the degree of the denominator, and n is the sample number.It is implemented using the following difference equations(assuming M = N):
///
///a[0]* y[n] = b[0] * x[n] + d[0][n - 1]
///d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n - 1]
///d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n - 1]
///...
///d[N - 2][n] = b[N - 1]* x[n] - a[N - 1]* y[n] + d[N - 1][n - 1]
///d[N - 1][n] = b[N] * x[n] - a[N] * y[n]
///d[N] = 0</remarks>
public static double[] lfilter(double[] x, double[] a, double[] b, double[] zi)
{
    if (a.Length != b.Length)
        throw new ArgumentOutOfRangeException("A and B filter coefficents should have the same length");
    double[] y = new double[x.Length];
    int N = a.Length;
    double[] d = new double[N];
    zi.CopyTo(d, 0);
    for (int n = 0; n < x.Length; ++n)
    {
        y[n] = b[0] * x[n] + d[0];
        for(int f = 1; f < N; ++f)
        {
            d[f - 1] = b[f] * x[n] - a[f] * y[n] + d[f];
        }
    }
    return y;
}


/// <summary>
/// Apply a digital filter forwards and backwards to have a zero phase shift filter
/// </summary>
/// <param name="data">The data to filter</param>
/// <param name="a">Filter coefficents a</param>
/// <param name="b">Filter coefficents b</param>
/// <returns>The filterd data</returns>
/// <remarks>
/// In order to prevent transients at the end or start of the sequence we have to pad it
/// The padding is done by rotating the sequence by 180° at the ends and append it to the data
/// </remarks>
public static double[] FiltFilt(double[] data, double[] a, double[] b, double[] zi)
{
    //Pad the data with 3 times the filter length on both sides by rotating the data by 180°
    int additionalLength = a.Length * 3;
    int endOfData = additionalLength + data.Length;
    double[] x = new double[data.Length + additionalLength * 2];
    Array.Copy(data, 0, x, additionalLength, data.Length);
    for (int i = 0; i < additionalLength; i++)
    {
        x[additionalLength - i - 1] = (x[additionalLength] * 2) - x[additionalLength + i + 1];
        x[endOfData + i] = (x[endOfData - 1] * 2) - x[endOfData - i - 2];
    }
    //Calculate the initial values for the given sequence
    double[] zi_ = new double[zi.Length];
    for (int i = 0; i < zi.Length; i++)
        zi_[i] = zi[i] * x[0];
    double[] y = lfilter(x, a, b, zi_);
    double[] y_flipped = new double[y.Length];
    //reverse the data and filter again
    for (int i = 0; i < y_flipped.Length; i++)
    {
        y_flipped[i] = y[y.Length - i - 1];
    }
    for (int i = 0; i < zi.Length; i++)
        zi_[i] = zi[i] * y_flipped[0];
    y = lfilter(y_flipped, a, b, zi_);
    y_flipped = new double[data.Length];
    //rereverse it again and return
    for (int i = 0; i < y_flipped.Length; i++)
    {
        y_flipped[i] = y[endOfData - i - 1];
    }
    return y_flipped;
}

Answer 1

因此，事实证明 IIR 滤波器对给定的参数和数值精度非常敏感。我从一个文本文件中导入了我的过滤器系数，我认为这已经足够好了，但事实并非如此。从二进制文件导入会得到我想要的结果，我的算法会产生与 scipy filtfilt 相同的输出。

为了完整起见，用于从 python 导出参数并将它们导入 C# 的函数：

import sys
import os;
import numpy as np
import struct

import scipy.signal as signal

order = 4
lowF = 0.2
highF = 0.8

b, a = signal.butter(order, [lowF, highF], btype='band')

#b, a = signal.iirfilter(float(sys.argv[4]), float(sys.argv[5]), btype='lowpass')

path = os.path.dirname(os.path.abspath(sys.argv[0])) + '\\';

# Get the steady state of the filter's step response.
zi = signal.lfilter_zi(b, a)

f = open(path + 'parameterfile.bin',"wb") 

def writeToFile(array, file):
    file.write(struct.pack("<I", array.shape[0]))
    print('Byte length: ' + str(len(struct.pack("<I", array.shape[0]))))
    for i in range(array.shape[0]):
        file.write(bytes(array[i]));

writeToFile(a, f)
writeToFile(b, f)
writeToFile(zi, f)

f.close();

还有 C# 代码：

private static void GetFilterAndZiFromBin(string path, out double[] a, out double[] b, out double[] zi)
{
    try
    {
        List<double> a_ = new List<double>();
        List<double> b_ = new List<double>();
        List<double> zi_ = new List<double>();
        FileStream fs = new FileStream(path,
                            FileMode.Open,
                            FileAccess.Read);
        BinaryReader br = new BinaryReader(fs);
        int length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for a from file");
        while(length > 0)
        {
            a_.Add(br.ReadDouble());
            length--;
        }
        length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for b from file");
        while (length > 0)
        {
            b_.Add(br.ReadDouble());
            length--;
        }
        length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for zi from file");
        while (length > 0)
        {
            zi_.Add(br.ReadDouble());
            length--;
        }
        a = a_.ToArray();
        b = b_.ToArray();
        zi = zi_.ToArray();
    }
    catch (Exception e)
    {
        a = new double[0];
        b = new double[0];
        zi = new double[0];
        throw e;
    }

}