Господа, посоветуйте чайнику библиотеку расчета (синтеза) цифровых фильтров, совместимую с дельфями.
На входе задал параметры фильтра: тип, частота среза, порядок - на выходе коэффициенты.

Применение. Пишу прогу для писюка, работающую со звуковушкой. Частота дискретизации звуковушки может переустанавливаться - частота среза ранее посчитанного фильтра естестьвенно сдвинется, надо пересчитывать. Также будет реализация фильтров с параметрами, задавемыми пользователем - например, частота резонансного пика - тоже надо рассчитывать на лету.

Intel Integrated Performance Primitives

Чайники обычно занимаются чайниками.

http://www.basegroup.ru/download/freecompo.../filteringbase/

Поищите по словам FIRDsgn, в исходниках на С - функция расчёта фильтров до десятка строк на каждый вид. И без .dll обойдётесь.
Не найдёте - в аттачменте возьмите.

Спасибо. Все понятно, код легко можно выдернуть и т.п.
Только там КИХ-фильтры, а нужны и КИХ, и БИХ. В первую очередь БИХ, забыл конкретизировать.
Попробую погуглить IIRDsgn.



Спасибо, не то.
В той компоненты реализованы:
1) Фильтр Калмана, в настоящий момент абсолютно мне не нужный.
2) Какой-то сложный универсальный ФНЧ, базирующийся на ПФ.

Мне же нужна процедура расчета классических КИХ и БИХ фильтров - такая, например, как в Matlab\FDATool, только реализованная в виде dll-ки без всяких гуёвых интерфейсов.




Тут, если я правильно понял краткое описание, есть лишь реализация самих алгоритмов _обработки_. Мне же нужна процедура расчета. Поправьте pls, если я понял неправильно, а то там наверное копать - не перекопать.

http://electronix.ru/forum/index.php?s=&am...st&p=856478

Спасибо. Руки дошли, качнул FDes.rar (f_design.cpp), прикрутил.
Только результаты не совпадают с моими контрольными данными. Может, есть дока или кто юзал эту библиотеку?

Какие еще контрольные данные? Задается шаблон фильтра. По рассчитанным к-там строится ачх и смотрится факт ее вписывания/не вписывания в шаблон. На моей памяти, все было нормально с этой либой.

Ну я вообще сравнивал по коэффициентам. Очень уж сильный разнобой. В случае с FDesign вообще получил коэффициенты больше единицы, для БИХ-фильтра разве возможно? Наверное, не так её юзаю.

Контрольный пример.

В WinFilter задаю ФНЧ с БИХ, порядок 2, Баттерворт, частота дискретизации 44100, частота среза 9кГц.
Получаю коэффициенты:
float ACoef[NCoef+1] = {
0.21320719693741449000,
0.42641439387482899000,
0.21320719693741449000
};

float BCoef[NCoef+1] = {
1.00000000000000000000,
-0.33915118494314400000,
0.19197997269398712000
};
Строю АЧХ. Соответствует.

Теперь по FDesign:

В переменной типа Filt_Params задаю:
apass1:=-1;
apass2:=-1;
astop1:=-80; //В децибелах? Иначе ругается.
astop2:=-80;
wpass1:=9000; //В герцах или круговая частота?
wpass2:=9000;
wstop1:=9100;
wstop2:=9100;
fsamp:=44100;
gain:=0;
acoefs:=@aCoeff;
bcoefs:=@bCoeff;
order:=2;
select:=#0;
approx:=#0;
implem:=#0;

Это задавалось из дельфей, но переменные передаются правильно, цеплялся дебаггером.
Вызываю функцию Calc_Butter_Coefs (оно?).
Получаю ACoeffs: (0, 0, 1.96522672836027), BCoeffs: (1, 1.98253712618972, 1.96522672836027)

Может, глянешь, есть под рукой примеры применения этой чудо-либы?

Если фильтр полосовой, надо
select = 'P' - bandPass
implem = 'I' - IIR
approx='B' - Butterworth
order = 2 - 0 - сам оценит порядок под селективные требования
fsamp=44100
apass1=-1 - в дБ
astop1 = -80
astop2 = -80
wpass1 = - круговая
wpass2 = - круговая
wstop1 = - круговая
wstop2 = - круговая
Calc_Filter_Coefs (оно).

Нужен ФНЧ.
Пробовал:
apass1:=-3; //дБ
apass2:=-3;
astop1:=-5; //дБ. Нужен фильтр с низким спадом, приблизительный аналог аналогового ФНЧ 2 порядка
astop2:=-5;
wpass1:=9000*2*3.14; //9кГц
wpass2:=9000*2*3.14;
wstop1:=10000*2*3.14;
wstop2:=10000*2*3.14;
fsamp:=44100;
gain:=0;
acoefs:=nil;
bcoefs:=nil;
order:=0;
select:='L';
approx:='B';
implem:='I';
end;
RetVal:=Calc_Filter_Coefs(@Params);

Короче, выдал коэффициенты:
Порядок=3
acoeffs 1 1 0 1
bcoeffs 1 -0.145210396025875 0 1

Такой фильтр выдаёт у меня клиппинг.

Чтобы не наслаивать проблемы стыковки с дельфями, залезал дебаггером в сишную dll-ку.
Ниже трассировка:

Calc_DigIIR_Coefs
Calc_Filter_Order
order=3

Calc_Normal_Coefs
Calc_Butter_Coefs
acoeffs 0 0 1.00079180388367 0
bcoeffs 0 1 1.00079180388367 1

Unnormalize_Coefs
Unnorm_LP_Coefs
acoeffs 0 0 65832.8315322196 0
bcoeffs 0 1 65832.8315322196 1

Bilinear_Transform
acoeffs 1 1 0 1
bcoeffs 1 -0.145210396025875 0 1

UnWarp_Freqs

либа выдает коэффициенты биквадратных звеньев. иначе говоря, для 3-го порядка будет 2 по 3 к-та a и 2 по 3 к-та b. здесь a - к-ты числителей, b - к-ты знаменателей. ну и результат надо умножить на gain или умножить к-ты числителей на корень из gain.

Спасибо.
Похоже, проще самому эту кухню освоить и прописать. Чем и займусь.
Еще вопрос-офтопик. Не порекомендует ли кто прогу построения АЧХ/ФЧХ цифровых фильтров? Чтоб сунул ей коэффициенты, и вперед. С удобнымм масштабированием и т.п. WinFilter, Matlab/FDATool видел.

я не в курсе. всю жизнь пользуюсь матлабской функцией freqz.

Возвращаясь к либе. Не совсем понял эту фразу.
Не пояснишь ли тупому, можно ли из этих коэффициентов биквадратных звеньев перейти к обычным коэффициентам, что выдаёт например Matlab/FDATool, WinFilter.
Перемножить эти коэффициенты на gain? У меня, кстати, gain получился в этом примере 0.103367640320208, не обратил сразу внимания.
Перемножил, получилось что-то похожее, но не совсем то.

Вообще-то рекурсивные фильтры реализуют именно последовательным или параллельным включнием звеньев 2-го порядка.
Ну а если в лоб:


Все получается как нужно.

Проcти, не понял, как ты строишь матрицу sos.
В моём примере порядок = 3. векторы a и b индексируются от 0 до 3-х.
a = [1 1 0 1]
b = [1 -0.145210396025875 0 1]
gain=0.103367640320208

Откуда взялись индексы a(3:5)? На первую строку кода Матлаб ругается "Subscript indices must either be real positive integers or logicals.".

1. бих-фильтры на практике реализуются соединением звеньев 2-го порядка
2. фильтр 3-го порядка содержит 2 звена: 2-го и 1-го, который обычно реализуется как второго с нулевыми к-тами при 2-й степени. итого - 2 биквадратных звена.
3. биквадратное звено имеет 3 к-та числителя и 3 к-та знаменателя (к-т 1 при квадрате тоже обычно присутствует для наглядности) итого 2 по 3 к-та числителей и 2 по 3 к-та знаменателей.
4. матлаб индексирует массивы начиная с 1. я привел в скрипте си-индексацию, которая естественно даст ошибку.

Вот мой скрипт Матлаба к данному примеру:
clear;
acoeff=[1 1 0 1 0 0]; %старшие коэффициенты заполняем нулями
bcoeff=[1 -0.145210396025875 0 1 0 0];
gain=0.103367640320208;
sos=[acoeff(1:3) bcoeff(1:3); acoeff(4:6) bcoeff(4:6)]
[b,a]=sos2tf(sos,gain)

Результаты:
sos =
1.0000 1.0000 0 1.0000 -0.1452 0
1.0000 0 0 1.0000 0 0

b =
0.1034 0.1034 0 0

a =
1.0000 -0.1452 0 0

АЧХ не та. Получилось наоборот - подъем по частоте.

Кстати, код фильтра использую тот, что даёт WinFilter:

float iir(float NewSample) {
float ACoef[NCoef+1] = { 1.0000, -0.1452, 0, 0 };

float BCoef[NCoef+1] = { 0.1034, 0.1034, 0, 0 };

static float y[NCoef+1]; //output samples
static float x[NCoef+1]; //input samples
int n;

//shift the old samples
for(n=NCoef; n>0; n--) {
x[n] = x[n-1];
y[n] = y[n-1];
}

//Calculate the new output
x[0] = NewSample;
y[0] = ACoef[0] * x[0];
for(n=1; n<=NCoef; n++)
y[0] += ACoef[n] * x[n] - BCoef[n] * y[n];

return y[0];
}


Это код именно для передаточной функции или я что-то напутал?

к-ты не надо изобретать. их либа честно посчитает. нужно тупо скопировать в матлаб.

bcoeff=[ 1 -0.14487845197425167 0 1 -0.38365823972570179 0.35186127061953859];
acoeff=[1 1 0 1 2 1];


код есть разностное уравнение для всего фильтра. никакого каскадирования тут нет.

Скопировал эти коэффициенты в Matlab\FDATool через "Import filter".
Matlab сказал, что фильтр нестабильный, и АЧХ совсем не та. Почему например один из acoeff в этом примере 2? Явно лажовые коэффициенты.

Да, это работает.
Хорошо, а как эти коэффициенты перевести в формат Numerator/Denominator, что используется в FDATool?
У меня он не воспринимает ни acoeff/bcoeff, ни a/b.

1 b->numerator
2 a->denominator
3 44100 -> Fs
4 import filter
5 PROFIT...

Прошло...
Осталось софтверно реализовать sos2tf()

ну, если делать больше нечего...