реализация IIR фильтра на целых числах Опции

[baralgin]

 По мотивам топика 75610. Значит из матлаба получаем некоторые коэффициенты во второй форме. Для простоты возмём односекционный фильтр, матлаб выдаёт следующее:

CODE

const int NL[MWSPT_NSEC] = { 1,3,1 };
const int16_T NUM[MWSPT_NSEC][3] = {
{
303, 0, 0
},
{
16384, 0, -16384
},
{
16384, 0, 0
}
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const int16_T DEN[MWSPT_NSEC][3] = {
{
16384, 0, 0
},
{
16384, -32152, 15778
},
{
16384, 0, 0
}
};

Написал функцию(точнее класс, но не важно):

Код

class BiquadSection
{
public:
    BiquadSection(){ W1 = 0; W2 = 0; }
    int operator << (int x0)
    {
        x0 *= 303;
        int W = x0 - ( -32152 * W1 + 15778 * W2 );

        W /= 16384;
        int ret = W - W2;

        W2 = W1;
        W1 = W;
        return ret;
    }
protected:
    int W1, W2;
};

Фильтр работает, но результат получается немного "квадратным"(рисунок). Такой же фильтр но в первой форме(direct form I) работает отлично(без видимых искажений). Где у меня косяк в вычислениях?

На скриншоте верхня синусоида это исходный сигнал, а нижняя фильтрованный. Апмлитуда исходного 460 единиц. Фильтр убирает постоянную составляющую.

Эскизы прикрепленных изображений

 

[baralgin]

Мне пришлось использовать "4096" (не знаю как научно это число обозвать...) - только так нет перепелонений. Я пробовал тоже запоминать остаток, но вышло быстрее вычислять его на ходу(т.е. W-стэйты вообще не делить). Вот моя реализация фильтра, который мне требовался несколько постов выше:

Код

/*
    Необходимо:
        входной гэйн делать заранее,
        выходное значение соответственно делить
*/
template <int Fraction, int K1, int K2>
class OneBiquadSection
{
public:
    OneBiquadSection(int init0 = 0)
    {
        W1 = init0;
        W2 = init0;
    }
    int operator << (int x0)
    {
        int W = x0, WI, WR;

        WI = (W1 / Fraction);
        WR = (W1 % Fraction);
        W -= K1 * WI + (K1 * WR) / Fraction;

        WI = (W2 / Fraction);
        WR = (W2 % Fraction);
        W -= K2 * WI + (K2 * WR) / Fraction;

        int ret = (W - W2);

        W2   = W1;
        W1   = W;
        return ret;
    }
protected:
    int W1;
    int W2;
};


// Касад из двух секций
class FilterBP4
{
public:
    FilterBP4(int init0 = 0)
        : Section1(init0), Section2(init0)
    {
    }
    int operator << (int X0)
    {
        X0 *= PreGain1;
        X0 = Section1 << X0; // первая секция

        X0 *= PreGain2;
        X0 /= Fraction;
        X0 = Section2 << X0; // вторая секция

        X0 /= Fraction;
        return X0;
    }
protected:
    static const int Fraction = (1 << 12); // 4096
    static const int PreGain1 = 74;
    static const int PreGain2 = 75;

    OneBiquadSection<Fraction, -8022, 3934> Section1;
    OneBiquadSection<Fraction, -8138, 4043> Section2;
};

На обработку одного канала уходит около около 90 тактов(без wait-стэйтов) - много, но терпимо(на float'х будет явно больше). Собственно самое "интересное" происходит тут:

Код

    WI = (W1 / Fraction);
    WR = (W1 % Fraction);
    W -= K1 * WI + (K1 * WR) / Fraction;

Страшно, но зато полученный фильтрованный сигнал гладкий даже на сверхмалых амплитудах.