Как реализовать БИХ-фильтр в C?

Question

Я пытаюсь реализовать БИХ-фильтр.

Я пытался реализовать следующий фильтр, но

FIR : y(n) = b0(x[n]) + ... +bM-1(x[n-M+1])

IIR : y(n) = {b0(x[n]) + ... +bM-1(x[n-M+1])} - {a1(y[n-1]) + ... +aN(y[n-N}

Я не понимаю, как реализовать y[n-1]......

Вот мой код.

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *coeffa, *coeffb;
    double *inputp;
    double *outputp;
    int n,k;

    //filter length =7
    for (n = 0; n < length; n++) {
        coeffa = coeffs_A;
        coeffb = coeffs_B;
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) {
            bcc += (*coeffb++) * (*inputp--); //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 1; k < filterLength; k++) {
            acc += (*coeffa++) * (*output--); //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

Я не буду здесь копировать код функции memove для краткости.

Answer 1

БИХ-фильтр является рекурсивным, что означает, что он принимает прошлые значения вывода (эквивалентно, что он может быть выражен как бесконечная последовательность ввода).

Предположим, у вас есть фильтр

y[n] = b*x[n]-a*y[n-1]

Обычно первый выход инициализируется заданным значением, например 0.

Предположим, что вы хотите отфильтровать сигнал x длины N, тогда вы должны сделать что-то вроде:

double out[N]; //output vector
double a = 0.3, b=0.5; //assign the coefficients a value
out[0] = 0; //initialize the first element

for (int i=1; i<N; i++)
{
    out[i] = b*x[i] -a*[i-1];  
}

В случае вашего кода я не могу знать, что вы делаете в строке inputp = &insamp[filterLength - 1 + n];, и это может быть проблемой.

Я предполагаю, что inputp - это сигнал, который вы хотите отфильтровать.

Другая проблема: вы использовали фильтр filterLength, чтобы указать длину как входных, так и выходных элементов:обычно отсутствует в БИХ-фильтре.

Элементы вывода от 0 до filterlength должны быть как-то инициализированы, предположим, к 0. Затем во втором цикле вы сделали индекс цикла, начиная с 1, но массив коэффициентов должен начинаться с 0.

(DOUBT: Если самый старый элемент y[5], как длина может быть 7?)

Используя индекс вместо разыменования массива, ваш код должен выглядеть примерно так:

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *inputp;
    int n,k;


    for (int ii=0; ii<filterLength; ii++)
    {
        output[ii] = 0;
    }

    //filter length =7
    for (n = 0; n < length; n++) {
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) 
        {
            bcc += coeffb[k] * inputp[filterLength-k-1]; //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 0; k < filterLength; k++) 
        {
            acc += coeffa[k] * output[filterLength-k-1]; //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

РЕДАКТИРОВАТЬ Я думаю, что два for, имеющие одинаковую длину, могут быть объединены вместе:

for (k = 0; k < filterLength; k++) 
{
    output[n] += (coeffb[k] * inputp[filterLength-k-1] - coeffa[k] * output[filterLength-k-1]); 
}

Answer 2

Если вы действительно хотите сделать это с помощью указателей:

void filter1(const double *b, const double *a, size_t filterLength, const double *in, double *out, size_t length) {
    const double a0 = a[0];
    const double *a_end = &a[filterLength-1];
    const double *out_start = out;
    a++;
    out--;
    size_t m;
    for (m = 0; m < length; m++) {
        const double *b_macc = b;
        const double *in_macc = in;
        const double *a_macc = a;
        const double *out_macc = out;
        double b_acc = (*in_macc--) * (*b_macc++);
        double a_acc = 0;
        while (a_macc <= a_end && out_macc >= out_start) {
            b_acc += (*in_macc--) * (*b_macc++);
            a_acc += (*out_macc--) * (*a_macc++);
        }
        *++out = (b_acc - a_acc) / a0;
        in++;
    }
}

Я сравнил результат этого алгоритма с функцией фильтра MATLAB.

Примечание : Вы можетеполучить большой прирост производительности, если вы нормализуете свои коэффициенты (т.е. a0 == 1).Чтобы сделать это, вы просто делите свои a и b векторы на a0, и тогда вам не нужно делить b_acc - a_acc на a0 на каждой итерации.

Answer 3

y[n-1] - это просто результат предыдущего временного шага, например: output[n-1] (y[n-N] для N-го временного шага от текущего).Таким образом, вам нужно индексировать в выходной массив на сумму, соответствующую output[n-k] (?) (Будьте осторожны, чтобы избежать индексации за пределами массива, однако вам понадобится код для предотвращения этого).

Кроме того, вы неправильно индексируете свои коэффициенты coeffa / coeffb (я думаю).Чтобы получить желаемый результат, вам может потребоваться сделать это * (coeffa ++), чтобы убедиться, что вы увеличиваете указатель, прежде чем разыменовываете его.Однако я бы настоятельно рекомендовал вам сделать это путем непосредственного индексирования в массиве, например, coeff_A[k], поскольку гораздо легче увидеть, что происходит.

Последнее замечание: операция (* output--) перемещает указатель вывода способом, который я считаю непреднамеренным.Индекс в массив вместо.