Реализация БИХ-фильтра нижних частот в C++

У меня есть коэффициенты фильтра для фильтра нижних частот Баттерворта, полученные из функции Matlab [b, a] = butter(3, 0.4, 'low'), и я реализовал те же вычисления, что и Matlab, согласно их документацию для filter(b, a, X). Результаты фильтрации постоянного сигнала 5.0, например, такие же, но только для первых 10 значений!?

Я предполагаю, что мой кольцевой буфер неправильный, но я не могу найти никаких проблем. Значения в x прописаны правильно методом фильтра, массивы инициализированы нулями, указатель циклического буфера n имеет правильные значения... Есть идеи?

// Interface
class LowpassFilter {
private:
    double x[10]; // input vector
    double y[10]; // output vector
    int n;    // pointer to the current array index

public:
    LowpassFilter(); 
    double filter(double sample);
};


// Implementation

// filter coefficients a and b
const double a[] = {1.0, -0.577240524806303, 0.421787048689562, -0.056297236491843};
const double b[] = {0.098531160923927, 0.295593482771781, 0.295593482771781, 0.098531160923927};
static int c = 0;

LowpassFilter::LowpassFilter() : x{0.0}, y{0.0}, n(0) { } // Constructor

double LowpassFilter::filter(double sample)
{
    x[n] = sample;
    y[n] = b[0] * x[n] + b[1] * x[(n-1)%10] + b[2] * x[(n-2)%10] + b[3] * x[(n-3)%10]
                       - a[1] * y[(n-1)%10] - a[2] * y[(n-2)%10] - a[3] * y[(n-3)%10];

    std::cout << c++ << ": " << y[n] << std::endl; // for checking the result with the Matlab results

    double result = y[n];
    n = (n + 1) % 10; // new pointer index 
    return result;
}

person Michael Dorner    schedule 20.03.2014    source источник
comment
Что происходит с модулем, когда, скажем, n == 0 и (n-2)% 10? Проверьте это.   -  person emsr    schedule 20.03.2014
comment
Что вы выводите? В чем проблема, которую вы видите.   -  person emsr    schedule 20.03.2014
comment
@MichaelDorner: модуль действительно имеет отрицательные значения, и это причина проблемы. Замените, например. (n-1)%10 с (n+10-1)%10, чтобы всегда получать ожидаемое положительное значение.   -  person Mike Seymour    schedule 20.03.2014
comment
Спасибо большое, вы совершенно правы!   -  person Michael Dorner    schedule 21.03.2014


Ответы (1)


Спасибо Майку Сеймуру и emsr проблема заключалась в отрицательных индексах при вычислении y[n]. Для решения задачи нужно принять только одну строку:

y[n] = b[0] * x[n] + b[1] * x[(n-1+m)%m] + b[2] * x[(n-2+m)%m] + b[3] * x[(n-3+m)%m]
                   - a[1] * y[(n-1+m)%m] - a[2] * y[(n-2+m)%m] - a[3] * y[(n-3+m)%m];

чтобы убедиться, что индекс положительный. Теперь он работает нормально. Большое спасибо!

person Michael Dorner    schedule 21.03.2014