Значение рациональной передаточной функции, лежащей в основе фильтра MATLAB или фильтра Scipy.signal

У меня есть код MATLAB, который фильтрует входной сигнал, используя filter:

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);

Я пытаюсь преобразовать этот код в С#. У меня уже есть функция butter, которая работает довольно быстро, но теперь я застрял на преобразовании функции filter.

Я прочитал документацию по фильтру MATLAB и Python Scipy.signal filter, но есть термин, присутствующий в определении передаточной функции, которого я не понимаю.

Вот определение «рациональной передаточной функции» из связанной документации:

        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# filter matlab signal-processing transfer-function
person Bjqn    schedule 24.05.2018    source источник
comment
Я полагаю, что это может быть функция временной задержки, исх. эти заметки Калифорнийского технологического института, пример 6.2 и описание рисунка 6.8. Примечание. У меня нет опыта в этом, я только когда-либо использовал передаточную функцию в Simulink, которая не требует функции X.   -  person Wolfie    schedule 25.05.2018
comment
X — входной сигнал, X(z) — z-преобразование X, а z — оператор единичного сдвига.   -  person percusse    schedule 15.06.2018

Ответы (1)


arrow_upward
2
arrow_downward

В разделе More About документов Matlab, как вы указали, они описывать:

Описание ввода-вывода операции фильтрации вектора в области Z-преобразования представляет собой рациональную передаточную функцию. Рациональная передаточная функция имеет вид

        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 входного вектора x (см.Digital Filter). Важно отметить, что и в документах они дают альтернативное представление передаточной функции в виде 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;
}

Драйвер:

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;
}
person StaticBeagle    schedule 29.06.2018
comment
Спасибо за ответ и реализацию! Но это очень-очень МЕДЛЕННО. Мы говорим несколько минут (иш). - person Bjqn; 02.07.2018
comment
Размер моих входных данных (x) равен 240001, поэтому вложенный цикл while снижает производительность. - person Bjqn; 02.07.2018
comment
@Bjqn размер b и a фиксирован? - person StaticBeagle; 02.07.2018
comment
Да b и a всегда будут иметь длину 2 - person Bjqn; 03.07.2018
comment
@Bjqn Я добавил обновление к своему ответу, чтобы иметь оптимизированную версию, когда оба вектора коэффициентов имеют фиксированную длину 2. Если этот подход не работает, вам, возможно, придется попробовать другие подходы, см. этот ответ - person StaticBeagle; 03.07.2018
comment
Большое спасибо! Обновленная реализация работает действительно очень хорошо, и вывод максимально близок к MATLAB: D И я действительно хорошо понимаю обновленную версию и могу связать ее с измененной передаточной функцией. :D - person Bjqn; 04.07.2018
comment
@Bjqn круто! Рад, что смог помочь. - person StaticBeagle; 04.07.2018