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


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

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

Про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()

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