Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: реализация IIR фильтра на целых числах
Форум разработчиков электроники ELECTRONIX.ru > Сайт и форум > В помощь начинающему > Программирование
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 единиц. Фильтр убирает постоянную составляющую.
IWG
Косяка на самом деле вроде как и нет. Просто отношение частоты дискретизации к сигналу внесло свою лепту. Выход - либо увеличить частоту дискретизации, либо увеличить разрядность.
baralgin
Цитата(IWG @ May 21 2010, 06:34) *
Косяка на самом деле вроде как и нет. Просто отношение частоты дискретизации к сигналу внесло свою лепту. Выход - либо увеличить частоту дискретизации, либо увеличить разрядность.

Меня смущает что этот код работает:
Код
    int operator << (int x0)
    {
        int ret = ( x0 - x2 ) - ( ( -32152 * y1 + 15778 * y2 ) / 16384 );
        x2 = x1;
        x1 = x0;
        y2 = y1;
        y1 = ret;
        return (ret * 303) / 16384;
    }

Коэффициенты теже, сигнал тот же, а результат красивый(скриншот не привожу - там всё гладко). Хотя везде пишут что каноническая форма лучше...

Ладно, зада такая: частота сэмплирования 12.5кГц(меньше ну очень не хочется), 12bit, полезный сигнал 50Гц. Нужно убрать постоянную составляющую и немного пригладить высокочастотные шумы. В примерах выше полосный фильтр Баттерворта: Fcut1=25Гц, Fcut2=100Гц. Желательно, но не обязательно, иметь фазовый сдвиг поменьше в районе 50гц. Что можете посоветовать по такой задаче?
IWG
То есть этот код не ваш? (Если да, то чем не устраивает)
Код
    int operator << (int x0)
    {
        int ret = ( x0 - x2 ) - ( ( -32152 * y1 + 15778 * y2 ) / 16384 );
        x2 = x1;
        x1 = x0;
        y2 = y1;
        y1 = ret;
        return (ret * 303) / 16384;
    }

А ваш тот, который в составе класса приведенного выше? Я правильно понимаю? Потому как проверив на отсчетах результат разный => функция, объявленная в классе, переписана с ошибкой.
baralgin
IWG Нет, код мой во всех случаях. Вымученный многими часами плясок. В первом случае что-то с промежуточными данными(Wx), я так понимаю сказываются ошибки оругления после деления на 16384. Пробовал делить позже, сперва вычисляя результат - всё равно "квадратит".
baralgin
В общем с горем пополам разобрался с проблемой из первого поста: чтобы результат был гладким пришлось учитывать при расчётах остаток от деления W на 16384. Что в итоге привело к двукратному падению темпа выполнения, но меня это пока устраивает. Возникла другая проблема: размерности int не хватает(2-3 разряда) для вычисления промежуточного значения W ("int W = x0 - ( -32152 * W1 + 15778 * W2 )"). Что нужно учитывать при квантовании коэффициентов чтобы этого избежать?

ps: мк cortex-m3.
Serj78
Цитата(baralgin @ Jun 17 2010, 14:52) *
В общем с горем пополам разобрался с проблемой из первого поста: чтобы результат был гладким пришлось учитывать при расчётах остаток от деления W на 16384. Что в итоге привело к двукратному падению темпа выполнения, но меня это пока устраивает. Возникла другая проблема: размерности int не хватает(2-3 разряда) для вычисления промежуточного значения W ("int W = x0 - ( -32152 * W1 + 15778 * W2 )"). Что нужно учитывать при квантовании коэффициентов чтобы этого избежать?

ps: мк cortex-m3.


Cтолкнулся с похожей проблемой при реализации iir фнч на STM32. Поборол тем что делил на "16384" (у меня 32768 ) каждое слагаемое отдельно и запоминал остаток. Получилось значительно быстрее чем во float и double, хоть и пара лишних переменных добавилась.
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;

Страшно, но зато полученный фильтрованный сигнал гладкий даже на сверхмалых амплитудах.
ae_
Цитата(baralgin @ Jun 21 2010, 15:05) *
Код
    WI = (W1 / Fraction);
    WR = (W1 % Fraction);
    W -= K1 * WI + (K1 * WR) / Fraction;

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

Название темы "... на целых числах" и данный фрагмент, где используется дробная арифметика в формате с фиксированной точкой вида 16.16, причём целая и дробная часть - это разные переменные. Зачем так сложно?
Если работать с дробным представлением числа, то проще сразу работать с Long Int (32 bit), представляя его как 16.16. Одним целочисленным делением сразу получим и целую и дробную часть.
baralgin
Цитата(ae_ @ Jun 23 2010, 17:05) *
причём целая и дробная часть - это разные переменные. Зачем так сложно?

Так ведь в начале(первый пост) не было сложно, но работало неудовлетворительно, поэтому пришлось сознательно усложнять. Я собственно поэтому и обратился за помощью, т.к. представляю что делаю скорее всего что-то не правильно. Скорее всего действительно все беды из-за большой разницы между частотй сэмплирования и полезным сигналом, ... но решать то проблему надо smile.gif .

Цитата(ae_ @ Jun 23 2010, 17:05) *
Если работать с дробным представлением числа, то проще сразу работать с Long Int (32 bit), представляя его как 16.16. Одним целочисленным делением сразу получим и целую и дробную часть.

Каюсь, нужно было явно указать размерность типов: в моём случае int=int32_t. Тоесть в моём коде все переменные это знаковые 32 бита. А в этом коде:
Код
    WI = (W1 / Fraction);
    WR = (W1 % Fraction);

делений(div) нет вообще, но несколько лишних команд конечно же добавляет. Так писать пришлось из-за нехватки 32 бит - переполнялось.
Собственно кусок кода от старших товарищей мог бы помочь, только без специфических инструкций и особенностей dsp-процессоров.
xemul
int32_t у Вас может переполниться: 12 бит с АЦП Вы разгоняете PreGain'ами на 6+, потом коэффициентами фильтров ещё на 13+ бит, и один бит даст сумма (K1 * W1 + K2 * W2).
Для корректного округления при целочисленном делении стОит делать как-то так:
y = ( x + fraction/2) / fraction ;
baralgin
xemul С переполнениями понятно, коэффициенты пришлось уменьшать. Разобраться с fdatool'ом так и не удалось - иногда он отказывается реагировать на изменение параметров квантизации, да и контекстный хэлп там отвратный.
Насчёт округления:
Код
y = ( x + fraction/2) / fraction;
есть один момент: со знаковыми данными этот фокус не пройдёт и нужно проверять знак перед операцией. Можно и это добавить, но скорость ещё больше приблизится к плавающим собратьям smile.gif .
xemul
Цитата(baralgin @ Jun 24 2010, 11:24) *
Насчёт округления:
Код
y = ( x + fraction/2) / fraction;
есть один момент: со знаковыми данными этот фокус не пройдёт и нужно проверять знак перед операцией.

Знак можно проверять после операции.
Код
y = ( x + fraction/2) / fraction;
if(y<0) y--;
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2025 Invision Power Services, Inc.