Цифровой фильтр Бесселя

Все что касается фильтрации
INFERION
Сообщения: 34
Зарегистрирован: 20 окт 2015, 13:06
Откуда: Украина

Re: Цифровой фильтр Бесселя

Сообщение INFERION »

Цифровой БИХ фильтр с близкой к линейной ФЧХ. Минимально-фазовый. В природе он имеет аппроксимацию Бесселя:
Изображение
Этот фильтр прост и позволяет так же просто скомпенсировать его ГВЗ. Так же он себя отлично зарекомендовал в обработке аудио. Я хочу сравнить его реализацию классическим методом (каноническая форма БИХ-фильтра) с уже имеющейся у меня, и взять более оптимальный вариант. К тому же это поможет мне лучше разобраться в классическом методе их проектирования. Я удивился, когда в списке статей нашел кучу примеров фильтров Баттерворта и Чебышева, но не встретил третий популярный вариант аппроксимации - Бесселя... Он или не реализуем, или просто не востребован из-за каких-нибудь преимуществ КИХ фильтрации. Ну или просто статьи нет подходящей. Чтоб разобраться с этим моментом - я начал рыть данный форум.

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

...методы вычисления коэффициентов дискретной передаточной характеристики для меня подобны магии.
Возможно Вам помогут статьи: http://we.easyelectronics.ru/Theory/cif ... ast-1.html
Итак возвращаемся к обсуждению цифрового фильтра Бесселя.
Литература по фильтру Бесселя:
http://scask.ru/book_r_cos.php?id=81
Пример:
Аналоговый фильтр Бесселя 3 порядка:

Исходное выражение уже записано для некоторой частоты среза Ws=15^{1/3}. То есть преобразование ФНЧ-ФНЧ вида S=S/Ws уже сделано. Пересчитать на 1 рад/сек и потом на нужное значение beta легко :-)

Cделаем z- преобразование

Можно тупо это выражение подставить в формулу для H(S), но мне лично удобнее воспользоваться матрицей z- преобразования ФНЧ 3-го порядка:
Коэффициенты H(z)
CodeCogsEqn(16).gif
CodeCogsEqn(16).gif (2.09 КБ) 6882 просмотра

CodeCogsEqn(15).gif
CodeCogsEqn(15).gif (2.33 КБ) 6886 просмотров

Для Fs/FN=0.04 я получил вот такие коэффициенты:
Бессель таб2.JPG
Бессель таб2.JPG (13 КБ) 6886 просмотров
Бессель 3.JPG
Программа WinFilter даёт такие же результаты...

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

Давайте учтём ещё один член в разложении:


Пусть аналоговый ФНЧ-прототип имеет вид:

Найдем матрицу z - преобразования|K|:
CodeCogsEqn(17).gif
CodeCogsEqn(17).gif (2.89 КБ) 6874 просмотра
CodeCogsEqn(18).gif
CodeCogsEqn(18).gif (2.88 КБ) 6874 просмотра
CodeCogsEqn(19).gif
CodeCogsEqn(19).gif (2.53 КБ) 6874 просмотра

INFERION
Сообщения: 34
Зарегистрирован: 20 окт 2015, 13:06
Откуда: Украина

Re: Цифровой фильтр Бесселя

Сообщение INFERION »

Спасибо за статью, но, к сожалению, с самого начала она требует порог вхождения выше, чем имеется у меня:
1. Z=1/(jωC) - не понятно зачем влепили мнимую единицу вместо понятного и простого 2pi.
2. β - этот параметр мне не понятен. Ну определяет он частоту, а в чём выражается - непонятно. Как и частота, измеряющаяся не в Герцах, а почему-то в радианах/с. Это что, 2pi этих самых радианов равна одному периоду колебаний? Зачем такие сложности, чем понятные Герцы не угодили?
3. "Аналоговый ФНЧ Баттерворта 1 порядка" - ?.. Фильтр с критическим затуханием - только он может быть первого порядка, поскольку у одного реактивного компонента не может быть резонансной частоты. В составе такого фильтра нет колебательной системы, поэтому это просто RC или RL фильтр, с добротностью строго 0,5. Как к нему можно приписать аппроксимацию Баттерворта, имеющую добротность в 1/sqrt(2)=0,707?

В общем, я понял что мне потребуется время и разбирающийся в этих вопросах человек, чтоб объяснить понятным (не математику) языком, что из себя представляют те или иные единицы измерения и как их представлять в голове привыкшему к реальным физическим процессам человеку. И чем глубже я вникаю в данный вопрос, тем проще мне кажется текущая реализация системы. Сейчас нет необходимости пересчитывать коэффициенты под новую частоту среза такими сложными алгоритмами. Достаточно пересчитать номиналы двух реактивных компонентов под новую резонансную частоту, что делается через простую формулу расчёта резонансной частоты колебательного контура. А все эти матрицы, коэффициенты, преобразования... Нет, это знать, конечно же, полезно. Для реализации более сложных систем. Только вот сложно это для меня.

"Для Fs/FN=0.04 я получил вот такие коэффициенты:" - а в какую схему то они суются? Каноничную рекурсивную? Ей ведь требуется два набора коэффициентов. Почему их 8, если порядок 3-й (т.е. их у рекурсивной схемы должно быть 3 или 6)?
"Программа WinFilter даёт такие же результаты..." - а ГВЗ какая, относительно АЧХ? Она же решающая.

В общем, вопросов больше чем ответов...

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

"Для Fs/FN=0.04 я получил вот такие коэффициенты:" - а в какую схему то они суются? Каноничную рекурсивную? Ей ведь требуется два набора коэффициентов. Почему их 8, если порядок 3-й (т.е. их у рекурсивной схемы должно быть 3 или 6)?
Это коэффициенты H(z) - комплексной передаточной функции БИХ-фильтра, однозначно определяющая его АЧХ и ФЧХ. Соответственно коэффициентов должно быть именно 8. Коэффициенты H(z) также определяют рекуррентную формулу фильтра. В нашем случае эта формула будет иметь вид:
Y[n]=1/A0(B0 X[n]+B1 X[n+1]+B2 X[n+2]+B3 X[n+3]- A1 Y[n-1]- A2 Y[n-2]- A3 Y[n-3])
Z=1/(jωC) - не понятно зачем влепили мнимую единицу вместо понятного и простого 2pi.
Это комплексное (реактивное) сопротивление конденсатора. Для индуктивности Z=jwL. Как же тут без мнимой единицы?
"Аналоговый ФНЧ Баттерворта 1 порядка"

Можно считать, что на 1 порядке все типы фильтров (Баттерворта, Чебышева, Бесселя) совпадают.
...чтоб объяснить понятным (не математику) языком...
Ну не знаю... Я вроде в статьях специально старался выше 2- курса не подниматься...
Для начала советую по H(z) научиться строить АЧХ, ФЧХ и определять ГВЗ. Коэффициенты можно и в WinFilter определять. Хотя эта программа кажется все графики строит. Ну и конечно статьи на этом сайте почитать, хотя бы по цифровым фильтрам Баттерворта.

INFERION
Сообщения: 34
Зарегистрирован: 20 окт 2015, 13:06
Откуда: Украина

Re: Цифровой фильтр Бесселя

Сообщение INFERION »

Santik писал(а):Это комплексное (реактивное) сопротивление конденсатора. Для индуктивности Z=jwL. Как же тут без мнимой единицы?
Да вот так: Rl=2pi*f*L, Rc=1/(2pi*f*L). Те же сопротивления, только не учитывающие реактивный характер. Может быть, при выведении хитрых формул это и критично, но для расчёта физического фильтра удобнее сразу получать вещественный результат. Да и микроконтроллеру с ним работать проще, а характер сопротивления уже само собою в модели учитывается. Естественным образом.
Во всяком случае, мне удобнее было бы в статье видеть разницу в представлениях значений, чтоб понять что же это за графики с декартовой системой координат и какие-то кружочки, полюса, преобразования в что-то совершенно не физическое, что нельзя ни осциллографом пощупать, ни руками. Это какое-то математическое абстрагирование, что сразу в голову не лезет.
Можно считать, что на 1 порядке все типы фильтров (Баттерворта, Чебышева, Бесселя) совпадают.
Только в заблуждение вводит.
Для начала советую по H(z) научиться строить АЧХ, ФЧХ и определять ГВЗ.
Значит буду доставать своего друга, который уже магистр и которому тоже придётся вспоминать Z-преобразование. Чтоб объяснял мне непонятные моменты в горе закорлючек :). Да, так учёба и проходит. То, что забивают в головы - такой мусор, который надолго не задерживается без практического закрепления материала. Мне никогда раньше не требовалось это. Максимум - кватернионы, и-то только из-за гиростабилизированных платформ. Всё остальное рассчитывается куда более простыми методами.

Но это всё уже оффтоп, так что спасибо за внимание. Буду потихоньку въезжать в эту тему, почитывая статьи :). Мне бы интересно было ещё обсудить модель динамической головки, где-нибудь в соседней теме. Как её лучше построить и как вычислить необходимые параметры. Что скажете? Создавать тему?

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

но для расчёта физического фильтра удобнее сразу получать вещественный результат. Да и микроконтроллеру с ним работать проще, а характер сопротивления уже само собою в модели учитывается. Естественным образом.
Я вот этого не понял. Вы сейчас про аналоговые фильтры говорите или про цифровые?
Математику при записи цифрового фильтра проще использовать комплексную, но результаты (коэффициенты цифрового фильтра) всегда получаются действительными числами. Нет в конечной рекуррентной формуле фильтра никаких комплексных чисел!
z - преобразование при переходе от аналогового фильтра к цифровому -это обязательная операция. Вы какой-то другой способ построения цифрового фильтра знаете?
Рассмотрим RC- цепочку (ФНЧ) (или фильтр первого порядка)
Z= R+1/jwC; Uвых=Uc=Uвх/Z (1/jwC)=Uвх/(1+jRCw)
H(jw)=Uвых/Uвх=1/(1+jRCw)
Обозначим S=jw и Ws=1/RC
H(S)=Ws/(S+Ws)

Для перехода к цифровому фильтру надо сделать z-преобразование.
z=e^(jw)=e^S
Логарифмируем S=Ln(z) и ограничиваемся одним членом разложения Ln в ряд.
S=(z-1)/(z+1) =(1-z^{-1})/(1+z^{-1})
При z- преобразовании происходит деформация частотной шкалы. Всё это в статьях на сайте хорошо описано.

INFERION
Сообщения: 34
Зарегистрирован: 20 окт 2015, 13:06
Откуда: Украина

Re: Цифровой фильтр Бесселя

Сообщение INFERION »

Santik писал(а):Я вот этого не понял. Вы сейчас про аналоговые фильтры говорите или про цифровые?
Цифровые аналоги аналоговых фильтров. Ихние модели. Да, это немного иная реализация, поэтому мне и интересно понять классический способ их построения. У него больше возможностей.
Математику при записи цифрового фильтра проще использовать комплексную, но результаты (коэффициенты цифрового фильтра) всегда получаются действительными числами. Нет в конечной рекуррентной формуле фильтра никаких комплексных чисел!
Разница появляется, когда нужно реализовать перестраиваемый фильтр. Вот как мне пересчитывать коэффициенты на лету средствами МК? Так, как это показано в примерах с очень жирным Сишным кодом?
z - преобразование при переходе от аналогового фильтра к цифровому -это обязательная операция. Вы какой-то другой способ построения цифрового фильтра знаете?
Возможно. Но не думаю что это что-то новое. К примеру, вот это колебательный контур:

int K=256;
int C=32000;
int L=0;
while(1){
L+=C/K;
C-=L/K;

DAC=C;
}

Целочисленный, сразу дающий диффузию ошибки (скрывающий шум квантования псевдослучайным шумом последнего бита). Даёт готовый синус, на вычисление семпла уходит 4 арифметические операции (сложение, вычитание и два деления или умножения). Резонансная частота равна частоте дискретизации, делённой на 2pi*K (возможно, нужно ещё единичку куда-нибудь засунуть - я не проверял работу на близких к частоте дискретизации частотах). Как видите - я могу легко пересчитать коэффициент на месте. При чём амплитуда колебаний никак не меняется со временем несмотря на ошибки целочисленной системы (ошибка размазывается по спектру, просто так она не исчезает).
Хм, я у себя нашел фото тестирования аудио тракта трассоискателя на ATtiny25 таким кодом в нём:
P3120031.JPG
Теперь добавим к нему потери и получим RLC-фильтр с нужной нам аппроксимацией:

int K=256;
int R=K*Q
int C=0;
int L=0;
while(1){
L+=C/K+ADC;
C-=L/K+C/R;

DAC=C;
}

Если в качестве Q указать 0,707 (подставив вместо K*Q 181, т.к. компилятор сам не посчитает с такой переменной) - получится ФНЧ Баттерворта 2-го порядка. Я тут могу где-нибудь ошибиться т.к. сейчас создавать проект, чтоб проверить работу кода, не очень хочу, а забыть детали мог спокойно. Для фильтра Бесселя необходимо указать добротность в 0,57. Где-то у меня была пачка WAV файлов с тестовыми сигналами, пропущенными через такие фильтры... Наращивать порядок можно последовательным включением идентичных звеньев. Фазовые фильтры тоже реализуются. Нужно рыться в своих проектах, если пример нужен. Деление можно было бы заменить на умножение, но в целочисленной системе от этого толку нет. С высокой добротностью и набором таких колебательных систем можно получить не требующий никаких окон спектроанализатор (E=C^2+L^2) или узкополосный режекторный фильтр. От нечего делать можно и фортепиано синтезировать, которое спокойно потянет какой-нибудь 32-х битный МК. Устойчивость фильтра легко рассчитывается через добротность, максимальный уровень входного сигнала, и максимально допустимое значение переменных. При желании можно ограничить эти значения используя регистр статуса процессора, тогда просто получим клиппинг.
RC-фильтр: C-=(С-ADC)/R; - постоянная времени равна частоте дискретизации, делённой на R. Ну а из постоянной времени легко и частоту среза найти - Fcut=Fsmpl/(2pi*R)+1. Куда уж проще то?
Как видите - я обошелся без Z-преобразования...
Последний раз редактировалось INFERION 24 окт 2015, 14:04, всего редактировалось 2 раза.

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

Как видите - я обошелся без Z-преобразования...
Отлично!
Однако справедливости ради отмечу, что эти алгоритмы фильтрации не являются универсальными.
Вот универсальная рекуррентная формула фильтров 2-го порядка (Баттерворта, Чебышева, Кауэра, Бесселя и пр.)
Y[n]=b0 X[n]+b1 X[n-1]+b2 X[n-2] - a1 Y[n-1]-a2 Y[n-2]
Y[n] - выборки на выходе фильтра , X[n] - входные выборки
Формула требует 5 умножений.
Я бы Вам посоветовал на http://we.easyelectronics.ru/ статью написать на эту тему. Успех гарантирован ! :D
Кстати, там же можно обсудить модель динамической головки. Аудиофилов там полно...

Аватара пользователя
Santik
Сообщения: 609
Зарегистрирован: 28 дек 2010, 08:04
Откуда: Мирный (Якутия)
Контактная информация:

Re: Цифровой фильтр Бесселя

Сообщение Santik »

Рассмотрим подробнее программу, приведенную INFERION

Код: Выделить всё

int K=256;
int R=K*Q
int C=32000;
int L=0;
while(1){
L+=C/K+ADC;
C-=L/K+C/R;
DAC=C;
}
Хорошо было бы посмотреть какая АЧХ и ФЧХ у данного фильтра.
Пусть Q=0.707 По утверждению автора - это фильтр Баттерворта 2-го порядка.
Тогда R=181
... частоту среза найти - Fcut=Fsmpl/(2pi*R). Куда уж проще то?
Удобно с нормированными частотами работать Fcut/FN=1/(pi*R), где FN=Fsmpl/2 -частота Найквиста
Итак, при Q=0.707 получим Fcut/FN=0.00175951
Через чур низкочастотный фильтр получается. При Fsmpl=20 кГц частота среза 17.6 Гц
Ну ладно, проверим позже.
Возвратимся к тексту программы.
Обозначим X=ADC - входной сигнал Y=DAC=С - выходной сигнал.
L=L +Y[i-1]/K +X
Y=Y[i-1]-(L/K+Y[i-1]/R)=Y[i-1](1-1/R)-L/K
или
Y=Y[i-1](1-1/R)-(L[i-1]+Y[i-1]/K +X)/K
Продолжим последовательно определять L[i-2],L[i-3]...
Окончательно получим:
Y[i]=Y[i-1](1-1/R-1/K^2)- (Y[i-2]+Y[i-3]+...)/K^2 - (X[i]+X[2]+X[i-3]+....)/K
Из этой формулы следует, что комплексная передаточная характеристика фильтра Н(z) имеет вид:
... Блин, что-то редактор формул перестал работать. Потом допишу формулу.
Да и так понятно, что все коэффициенты Bi=1/K
A0=1;A1=1/R+1/K^2-1;A2=-1/K^2;A2=-1/K^2...
Зная Н(z) легко построить АЧХ и ФЧХ.
[attachment=0]Ф01.JPG[/attachment]
Вложения
Ф01.JPG

Ответить