Библиотека (dll) расчета фильтров, Посоветуйте pls Опции

[ymv]

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

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

[thermit]

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

[ymv]

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

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

[thermit]

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

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

Код

sos=[a(0:2) b(0:2); a(3:5) b(3:5)];
[b,a]=sos2tf(sos,gain);

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

Сообщение отредактировал thermit - Feb 17 2015, 16:04

[ymv]

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

Код

sos=[a(0:2) b(0:2); a(3:5) b(3:5)];
[b,a]=sos2tf(sos,gain);

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

Про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.".

[thermit]

Про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. я привел в скрипте си-индексацию, которая естественно даст ошибку.

[ymv]

Вот мой скрипт Матлаба к данному примеру:
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];
}


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

[thermit]

Вот мой скрипт Матлаба к данному примеру:
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];


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

[ymv]

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

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


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

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

[thermit]

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

Код

bcoeff=[ 1 -0.14487845197425167 0 1 -0.38365823972570179 0.35186127061953859];
acoeff=[1 1 0 1 2 1];
sos=[acoeff(1:3) bcoeff(1:3); acoeff(4:6) bcoeff(4:6)];
gain=0.103367640320208;
[b,a]=sos2tf(sos,gain);

plot(20*log10(abs(freqz(b,a,1000))));
abs(roots(a))

[ymv]

Код

bcoeff=[ 1 -0.14487845197425167 0 1 -0.38365823972570179 0.35186127061953859];
acoeff=[1 1 0 1 2 1];
sos=[acoeff(1:3) bcoeff(1:3); acoeff(4:6) bcoeff(4:6)];
gain=0.103367640320208;
[b,a]=sos2tf(sos,gain);

plot(20*log10(abs(freqz(b,a,1000))));
abs(roots(a))

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