matlab — Внедрение фильтра низких частот IIR при переполнении стека

Хей.

У меня есть коэффициенты фильтра для фильтра низких частот Баттерворта, исходя из функции 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;
}

2

Решение

Благодаря Майк Сеймур а также 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];

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

2

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]