| |
|
  |
 реализация IIR фильтра на целых числах |
|
|
|
|
 May 20 2010, 08:37
|
Частый гость
Группа: Участник
Сообщений: 92
Регистрация: 23-12-08
Из: Кишинёв
Пользователь №: 42 680
|
По мотивам топика 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 единиц. Фильтр убирает постоянную составляющую.
Эскизы прикрепленных изображений
|
|
|
|
|
|
|
|
May 21 2010, 03:34
|
Группа: Участник
Сообщений: 11
Регистрация: 7-09-07
Из: Омск
Пользователь №: 30 350
|
Косяка на самом деле вроде как и нет. Просто отношение частоты дискретизации к сигналу внесло свою лепту. Выход - либо увеличить частоту дискретизации, либо увеличить разрядность.
|
|
|
|
|
|
|
|
May 21 2010, 08:09
|
Частый гость
Группа: Участник
Сообщений: 92
Регистрация: 23-12-08
Из: Кишинёв
Пользователь №: 42 680
|
Цитата(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гц. Что можете посоветовать по такой задаче?
|
|
|
|
|
|
|
|
May 21 2010, 09:15
|
Группа: Участник
Сообщений: 11
Регистрация: 7-09-07
Из: Омск
Пользователь №: 30 350
|
То есть этот код не ваш? (Если да, то чем не устраивает) Код 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;
} А ваш тот, который в составе класса приведенного выше? Я правильно понимаю? Потому как проверив на отсчетах результат разный => функция, объявленная в классе, переписана с ошибкой.
Сообщение отредактировал IWG - May 21 2010, 09:26
|
|
|
|
|
|
|
|
Jun 21 2010, 06:05
|
Частый гость
Группа: Участник
Сообщений: 92
Регистрация: 23-12-08
Из: Кишинёв
Пользователь №: 42 680
|
Мне пришлось использовать "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; Страшно, но зато полученный фильтрованный сигнал гладкий даже на сверхмалых амплитудах.
|
|
|
|
|
|
|
|
Jun 23 2010, 14:05
|
Участник
Группа: Свой
Сообщений: 462
Регистрация: 2-04-07
Из: Иркутск
Пользователь №: 26 695
|
Цитата(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. Одним целочисленным делением сразу получим и целую и дробную часть.
|
|
|
|
|
|
|
|
Jun 23 2010, 15:40
|
Частый гость
Группа: Участник
Сообщений: 92
Регистрация: 23-12-08
Из: Кишинёв
Пользователь №: 42 680
|
Цитата(ae_ @ Jun 23 2010, 17:05) причём целая и дробная часть - это разные переменные. Зачем так сложно? Так ведь в начале(первый пост) не было сложно, но работало неудовлетворительно, поэтому пришлось сознательно усложнять. Я собственно поэтому и обратился за помощью, т.к. представляю что делаю скорее всего что-то не правильно. Скорее всего действительно все беды из-за большой разницы между частотй сэмплирования и полезным сигналом, ... но решать то проблему надо  . Цитата(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-процессоров.
|
|
|
|
|
|
|
|
Jun 24 2010, 07:24
|
Частый гость
Группа: Участник
Сообщений: 92
Регистрация: 23-12-08
Из: Кишинёв
Пользователь №: 42 680
|
xemul С переполнениями понятно, коэффициенты пришлось уменьшать. Разобраться с fdatool'ом так и не удалось - иногда он отказывается реагировать на изменение параметров квантизации, да и контекстный хэлп там отвратный. Насчёт округления: Код y = ( x + fraction/2) / fraction; есть один момент: со знаковыми данными этот фокус не пройдёт и нужно проверять знак перед операцией. Можно и это добавить, но скорость ещё больше приблизится к плавающим собратьям .
|
|
|
|
|
|
|
|
Jun 24 2010, 09:55
|
Группа: Свой
Сообщений: 1 928
Регистрация: 11-07-06
Пользователь №: 18 731
|
Цитата(baralgin @ Jun 24 2010, 11:24) Насчёт округления: Код y = ( x + fraction/2) / fraction; есть один момент: со знаковыми данными этот фокус не пройдёт и нужно проверять знак перед операцией. Знак можно проверять после операции. Код y = ( x + fraction/2) / fraction;
if(y<0) y--;
|
|
|
|
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|
