MATLAB 过滤器或 Scipy.signal 过滤器下的有理传递函数的含义

Question

我有一些使用 filter:

过滤输入信号的 MATLAB 代码

CUTOFF = 0.05;
FS = 5000;
[b, a] = butter(1, CUTOFF / (FS / 2), 'high');
% b = [0.99996859, -0.99996859]
% a = [1.0, -0.99993717]
dataAfter = filter(b, a, dataBefore);

我正在尝试将此代码转换为 C#。我已经让 butter 函数运行得非常快，但现在我无法转换 filter 函数。

我已阅读 MATLAB filter documentation and Python Scipy.signal filter 文档，但传递函数定义中有一个术语我不理解。

这是链接文档中的 "rational transfer function" 定义：

        b[0] + b[1]z^(-1)  + ... + b[M]z^(-M)
Y(z) = _______________________________________ X(z)
        a[0] + a[1]z^(-1)  + ... + a[N]z^(-N)

如果我错了请纠正我，但是 z 是输入数据的当前元素，而 Y(z) 是输出？

如果以上为真，这个等式中的X(z)是什么？

我想了解这个在C#中实现它，如果有等效的选项请赐教。

Answer 1

正如您所指出的，在 matlab 文档的 More About 部分中，他们描述了：

The input-output description of the filter operation on a vector in the Z-transform domain is a rational transfer function. A rational transfer function is of the form,

        b[0] + b[1]z^(-1)  + ... + b[M]z^(-M)
Y(z) = _______________________________________ X(z)
        a[0] + a[1]z^(-1)  + ... + a[N]z^(-N)

重新排列：

       Y(z)   b[0] + b[1]z^(-1)  + ... + b[M]z^(-M)
H(z) = ____ = _______________________________________
       X(z)   a[0] + a[1]z^(-1)  + ... + a[N]z^(-N)

因此，X(z) 是 z-domain transform of the input vector x (seeDigital Filter). It is important to mention that, also in the docs they give an alternate representation of the transfer function as a difference equation

哪个更适合移植到代码中。 C# 中的一种可能实现可能是 (using this answer as reference)

public static double[] Filter(double[] b, double[] a, double[] x)
{
   // normalize if a[0] != 1.0. TODO: check if a[0] == 0
   if(a[0] != 1.0)
   {
       a = a.Select(el => el / a[0]).ToArray();
       b = b.Select(el => el / a[0]).ToArray();
   }

   double[] result = new double[x.Length];
   result[0] = b[0] * x[0];
   for (int i = 1; i < x.Length; i++)
   {
       result[i] = 0.0;
       int j = 0;
       if ((i < b.Length) && (j < x.Length))
       {
           result[i] += (b[i] * x[j]);
       }
       while(++j <= i)
       {
            int k = i - j;
            if ((k < b.Length) && (j < x.Length))
            {
                result[i] += b[k] * x[j];
            }
            if ((k < x.Length) && (j < a.Length))
            {
                result[i] -= a[j] * result[k];
            }
        }
    }
    return result;
}

Driver:

static void Main(string[] args)
{
    double[] dataBefore = { 1, 2, 3, 4 };
    double[] b = { 0.99996859, -0.99996859 };
    double[] a = { 1.0, -0.99993717 };

    var dataAfter = Filter(b1, a, dataBefore);
}

输出

Matlab dataAfter = [0.99996859 1.999874351973491  2.999717289867956  3.999497407630634]
CSharp dataAfter = [0.99996859 1.9998743519734905 2.9997172898679563 3.999497407630634]

更新
如果系数向量 a 和 b 的长度固定为 2，则滤波函数可以简化为：

public static double[] Filter(double[] b, double[] a, double[] x)
{
    // normalize if a[0] != 1.0. TODO: check if a[0] == 0
    if (a[0] != 1.0)
    {
        a = a.Select(el => el / a[0]).ToArray();
        b = b.Select(el => el / a[0]).ToArray();
    }

    int length = x.Length;
    double z = 0.0;
    double[] y = new double[length];    // output filtered signal

    double b0 = b[0];
    double b1 = b[1];
    double a1 = a[1];
    for (int i = 0; i < length; i++)
    {
        y[i] = b0 * x[i] + z;
        z = b1 * x[i] - a1 * y[i];
    }
    return y;
}

MATLAB 过滤器或 Scipy.signal 过滤器下的有理传递函数的含义

Meaning of rational transfer function underlying MATLAB filter or Scipy.signal filter

c#

matlab

filter

signal-processing

transfer-function