[Pat 0]

Здравствуйте.

Вот предпринял очередную попытку сделать фильтрацию сигнала, уже не раз пытался и все мимо.

Ну никак не дается она моему пониманию ;)

Нашел вот здесь реализацию фильтра НЧ Баттеруорта 4 порядка

 

http://www.may.nnov.ru/mak/DSP/chBut4.shtml

 

Вот формула

 

Y = B[0]*X + ... + B[4]*X[i-4] - A[1]*Y[i-1] - ... - A[4]*Y[i-4]

 

Взял коэф. B и A те которые описаны там же.

Сделал небольшую программку на Visual C ( проект приложил).

В программе у меня формируется таблица синуса (256 точек).

Вообще то по задаче мне надо сделать ФНЧ 30 Гц . Необходимо выделить сигнал 25 Гц в котором присутствует 50 и 100 Гц.

 

Отсюда получается и таблица синуса 4 периода на 1 период 25 Гц = 100 Гц.

 

После чего преобразую эту таблицу (float) в целочисленные значения (int).

 

Дальше пропускаю этот массив через функцию float FIR_Filter(int val)

 

 

// Коэф. перевернуты т.е. в b[0] записано b[4] и наоборот

const float b[5] =

{

0.004824343358F,

0.019297373431F,

0.028946060146F,

0.019297373431F,

0.004824343358F

};

 

// Коэф. перевернуты т.е. в a[0] записано a[3] и наоборот

const float a[4] =

{

0.18737949237F,

-1.0546654059F,

2.3139884144F,

-2.3695130072F

};

 

float X[5] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };

float Y[5] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };

float res[NUM_SMP];

 

//-----------------------------------------------------------------------

// Main

//-----------------------------------------------------------------------

int main(int argc, char* argv[])

{

unsigned int i;

// Заполняем таблицу значениями синуса

// 256 выборок длиной 0.04 c (25 Гц) заполнена 100 Гц

fill_sin_table(100,sin_table);

// Преобразуем в int

fill_int_table(int_table,sin_table,100);

 

// Свертка

for(i = 0;i < NUM_SMP;i++)

{

 

res = FIR Filter(int_table);

 

printf("%d",int_table);

printf(" ");

printf("%f\n",res);

}

 

getchar();

return 0;

}// end main

 

//-----------------------------------------------------------------------

// Функция фильтрации

//-----------------------------------------------------------------------

float FIR_Filter(int val)

{

unsigned char i;

float sumb, suma;

sumb=0.0F;

suma=0.0F;

 

 

for (i=0;i<=3;i++)

{

X = X[i+1];

Y = Y[i+1];

sumb += X*b;

suma += Y*a;

}

X = val;

sumb += val*b;

 

Y[3] = sumb - suma;

 

return Y[3];

}// end FIR_Filter

 

Если я правильно понимаю принципы работы ЦФ то у меня на выходе res = FIR Filter(int_table);

Должны быть практически нули.

На самом деле в массиве res постоянно нарастают значения и где то на 96 шаге происходит переполнение.

 

Провел второй эксперимент взял Filter Solutions 10.0 и при помощи его сформировал ФНЧ с частотой выборки 6400 Гц и срезом 30 Гц.

Сгенерировал код и вставил полученную функцию и коэф. вместо своей функции float FIR_Filter(int val) и получил такое же переполнение.

 

Подскажите, пожалуйста, что не так я делаю. Вроде как все по теории.

 

Спасибо.

FIR.rar

[el34 0]

как я понимаю

Y = B[0]*X + ... + B[4]*X[i-4] - A[1]*Y[i-1] - ... - A[4]*Y[i-4]

не FIR, а IIR.....

 

а выделять 25Hz при помехе на 50Hz при помощи LowPass будет довольно накладно.....

большой порядок может потребоватся.....

думаю лучше поставить режекторы на 50 и 100 герц....

[Pat 0]

как я понимаю

Y = B[0]*X + ... + B[4]*X[i-4] - A[1]*Y[i-1] - ... - A[4]*Y[i-4]

не FIR, а IIR.....

 

Да не правильно дал название функции.

 

а выделять 25Hz при помехе на 50Hz при помощи LowPass будет довольно накладно.....

большой порядок может потребоватся.....

думаю лучше поставить режекторы на 50 и 100 герц....

 

А нельзя ли получить полную реализацию.

Дело в том что я пробовал и режекторный фильтр на 100Гц из Filter Solutions и та же ерунда с переполнением.

Тут где то одна ошибка идеологии построения или моего эксперимента.

[el34 0]

мой совет: сходите на

http://lord-n.narod.ru/walla.html#razdelDSP

http://dsp-book.narod.ru/

 

там много хороших книжек...

 

а вот тут :Расчет биквадратного фильтра

 

http://dsp-book.narod.ru/flt_calc.zip

[Pat 0]

мой совет: сходите на

http://lord-n.narod.ru/walla.html#razdelDSP

http://dsp-book.narod.ru/

 

там много хороших книжек...

 

Да ходил уже не раз, читал и даже чего то понял.

У меня реализация хромает.

То что вначале темы привел вроде все по теории, но не работает. Пусть там ФНЧ, но он же должен работать, а не переполнятся.

[Гость LordN]

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

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

нужно очень внимательно и вдумчиво продумать алгоритм и просчитать коэффициенты.

[Pat 0]

а вот тут :Расчет биквадратного фильтра

http://dsp-book.narod.ru/flt_calc.zip

 

За эту ссылку большой :a14:, вроде как работает по крайней мере в Visual C все получилось.

Счас буду пробовать реализацию в C8051F412.

 

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

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

нужно очень внимательно и вдумчиво продумать алгоритм и просчитать коэффициенты.

 

Наверное вы правы, меня просто сбивает с толку что в Visual C не хочет работать (там вроде как проблем с типами данных не должно быть) причем, очень странно что решение от Filter Solutions не работает.

 

А Вы не подскажите что же все таки лучше использовать в маленьких применениях КИХ или БИХ?

[el34 0]

Pat>А Вы не подскажите что же все таки лучше использовать в маленьких применениях КИХ или БИХ?

 

по этому поводу вспоминается фраза из анекдота:

"Вам шашечки или ехать?"

 

вообще везде пытаюсь применять FIR и только если никак

(по каким либо причинам)....переползаю на IIR....

и если есть ресурсы процессора (а они сейчас значительно дешевле чем - "думать") то конечно КИХ....

имхо

[Гость LordN]

только щас внимание обратил - я бы обязательно заменил float на double и в опциях компиллера включил бы максимальный размер дабла.

[Самурай 12]

float FIR_Filter(int val)

{

unsigned char i;

float sumb, suma;

sumb=0.0F;

suma=0.0F;

 

 

for (i=0;i<=3;i++)

{

X = X[i+1];

Y = Y[i+1];

sumb += X*b;

suma += Y*a;

}

X = val;

sumb += val*b;

 

Y[3] = sumb - suma;

 

return Y[3];

}// end FIR_Filter

 

 

 

Спасибо.

 

 

И еще наверно можно обратить внимание на строчку Y[3] = sumb - suma; Разве здесь индекс должен быть равен 3? Точно не 4? И наверно следствием этого является то, что когда i=3, то

suma += Y[3]*a[3], но при этом Y[3] =Y[4] = 0.