MATLAB/Scipy filter有理传递函数解析及C#实现问询

一、先搞懂X(z)到底是什么 #

你提到的传递函数里，X(z)是输入信号dataBefore的Z变换，Y(z)是输出信号dataAfter的Z变换。这里的z是Z变换的复变量，而z^(-1)其实是延迟算子，对应时域里的“前一个样本”——比如z^(-1)X(z)对应的就是输入信号的前一个样本x[n-1]，z^(-2)X(z)对应x[n-2]，以此类推。

把传递函数转成时域的递归公式（这才是MATLAB filter函数实际执行的逻辑）会更直观：

假设a[0] = 1（MATLAB的butter函数返回的a数组首元素默认就是1），那么：
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]
其中：

• x[n]是输入信号dataBefore的第n个样本
• y[n]是输出信号dataAfter的第n个样本
• 当n<0时，x[n]和y[n]都默认为0（这是MATLAB filter的默认零初始条件）

二、C#等效实现MATLAB filter函数 #

既然你已经实现了butter函数拿到了b和a数组，下面这个C#方法完全对应MATLAB的filter(b,a,dataBefore)行为，包括零初始条件的处理：

public static double[] Filter(double[] b, double[] a, double[] dataBefore)
{
    // 先确保a[0]为1（如果不是，就归一化b和a，和MATLAB行为一致）
    if (Math.Abs(a[0] - 1.0) > 1e-8)
    {
        double a0 = a[0];
        for (int i = 0; i < b.Length; i++)
            b[i] /= a0;
        for (int i = 0; i < a.Length; i++)
            a[i] /= a0;
    }

    int inputLen = dataBefore.Length;
    double[] dataAfter = new double[inputLen];

    // 缓存输入的延迟样本（x[n-1], x[n-2], ...）
    double[] xDelays = new double[b.Length - 1];
    // 缓存输出的延迟样本（y[n-1], y[n-2], ...）
    double[] yDelays = new double[a.Length - 1];

    for (int n = 0; n < inputLen; n++)
    {
        // 计算输入项的加权和
        double inputSum = b[0] * dataBefore[n];
        for (int i = 0; i < xDelays.Length; i++)
        {
            inputSum += b[i + 1] * xDelays[i];
        }

        // 计算输出项的加权和
        double outputSum = 0;
        for (int i = 0; i < yDelays.Length; i++)
        {
            outputSum += a[i + 1] * yDelays[i];
        }

        // 当前输出样本
        dataAfter[n] = inputSum - outputSum;

        // 更新输入延迟缓存：把当前x[n]移到缓存头部，旧数据后移
        for (int i = xDelays.Length - 1; i > 0; i--)
        {
            xDelays[i] = xDelays[i - 1];
        }
        if (xDelays.Length > 0)
            xDelays[0] = dataBefore[n];

        // 更新输出延迟缓存：把当前y[n]移到缓存头部，旧数据后移
        for (int i = yDelays.Length - 1; i > 0; i--)
        {
            yDelays[i] = yDelays[i - 1];
        }
        if (yDelays.Length > 0)
            yDelays[0] = dataAfter[n];
    }

    return dataAfter;
}

使用示例（对应你的MATLAB代码参数） #

// 你的Butterworth滤波器参数
double cutoff = 0.05;
double fs = 5000;
// 假设通过你的Butter函数得到以下b和a数组（和MATLAB输出一致）
double[] b = new double[] { 0.99996859, -0.99996859 };
double[] a = new double[] { 1.0, -0.99993717 };
// 你的输入信号数组
double[] dataBefore = new double[] { /* 填充你的输入数据 */ };
// 执行滤波
double[] dataAfter = Filter(b, a, dataBefore);

这个实现完全遵循MATLAB filter的逻辑，零初始条件，归一化处理（如果a[0]不是1的话），输出结果和MATLAB会完全一致（误差在浮点精度范围内）。

内容的提问来源于stack exchange，提问作者Bjqn