strap89 0 17 сентября, 2018 Опубликовано 17 сентября, 2018 · Жалоба Здравствуйте! Помогите разобраться пожалуйста. Есть цифровой БИХ-фильтр 2-го порядка, режекторный. Есть алгоритм расчета его коэффициентов. И есть пара коэффициентов, определяющих параметры фильтра (это мои предположения). Не могли бы Вы подсказать мне, какой тип фильтра тут использован и что определяют его параметры. Сам алгоритм: void makefiltr2coeff(filtr2& flt,double Tsam,double f0,double f1,double ksi0,double ksi1) { if(!f0||!f1||!Tsam||!ksi0||!ksi1) return; Tsam*=scaleTsam; double om0=2.*M_PI*f0; double om1=2.*M_PI*f1; double tau0=.5/tan(Tsam*om0/2.); double tau1=.5/tan(Tsam*om1/2.); double d=1.+(4.*tau1*(tau1+ksi1));if(!d) return; double b0=(1.+(4.*tau0*(tau0+ksi0)))/d; double b1=2.*(1.-4.*tau0*tau0)/d; double b2=(1.+(4.*tau0*(tau0-ksi0)))/d; double a1=2.*(1.-4.*tau1*tau1)/d; double a2=(1.+(4.*tau1*(tau1-ksi1)))/d; a0=1; Сама реализация фильтра есть, фильтр работает, но формулу расчета коэффициентов мне найти не удалось. Как соответственно и информацию, что означают его параметры. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Самурай 12 17 сентября, 2018 Опубликовано 17 сентября, 2018 · Жалоба Здравствуйте! Помогите разобраться пожалуйста. Есть цифровой БИХ-фильтр 2-го порядка, режекторный. Есть алгоритм расчета его коэффициентов. И есть пара коэффициентов, определяющих параметры фильтра (это мои предположения). Не могли бы Вы подсказать мне, какой тип фильтра тут использован и что определяют его параметры. Сам алгоритм: void makefiltr2coeff(filtr2& flt,double Tsam,double f0,double f1,double ksi0,double ksi1) { if(!f0||!f1||!Tsam||!ksi0||!ksi1) return; Tsam*=scaleTsam; double om0=2.*M_PI*f0; double om1=2.*M_PI*f1; double tau0=.5/tan(Tsam*om0/2.); double tau1=.5/tan(Tsam*om1/2.); double d=1.+(4.*tau1*(tau1+ksi1));if(!d) return; double b0=(1.+(4.*tau0*(tau0+ksi0)))/d; double b1=2.*(1.-4.*tau0*tau0)/d; double b2=(1.+(4.*tau0*(tau0-ksi0)))/d; double a1=2.*(1.-4.*tau1*tau1)/d; double a2=(1.+(4.*tau1*(tau1-ksi1)))/d; a0=1; Сама реализация фильтра есть, фильтр работает, но формулу расчета коэффициентов мне найти не удалось. Как соответственно и информацию, что означают его параметры. Ваш IIR 2-ого порядка получен через билинейное преобразование аналогового фильтра-прототипа общего вида: т.е. через замену и Выводить не буду, там все тривиально, но больно много писанины... Для режекторного фильтра необходимо выполнение условий: f0 = f1 = частоте подавления, ksi0 < ksi1, оба эти коэффициента влияют на уровень затухания и ширину полосы подавления, если ksi0=0, подавление максимальное. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
strap89 0 18 сентября, 2018 Опубликовано 18 сентября, 2018 · Жалоба Ваш IIR 2-ого порядка получен через билинейное преобразование аналогового фильтра-прототипа общего вида: т.е. через замену и Выводить не буду, там все тривиально, но больно много писанины... Для режекторного фильтра необходимо выполнение условий: f0 = f1 = частоте подавления, ksi0 < ksi1, оба эти коэффициента влияют на уровень затухания и ширину полосы подавления, если ksi0=0, подавление максимальное. Спасибо за объяснение. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Ruslan1 16 18 марта, 2019 Опубликовано 18 марта, 2019 · Жалоба Вопрос про усиление: у меня по этим формулам получается та же форма АЧХ, но только с ослаблением сигнала, от 0 dB и вниз, а если через Матлабовскую функцию- то все как у оригинала H(s), где-то на +30dB. Какую константу я забыл? Имею очень похожую формулу: H(S) = s^2 + 31.1 s + 799.4 --------------------- s^2 + 7.109 s + 25.27 если делаю в Матлабе: b = 1 31.1 799.4 a = 1 7.109 25.27 fs = 1000 [bd,ad] = bilinear(b,a,fs) то имею коэффициенты фильтра Quote bd = 1.0121 -1.9925 0.9812 ad = 1.0000 -1.9929 0.9929 И по ним fvtool строит именно то, что у аналогового прототипа (та же форма и то же усиление). Но вот если иду через вышеприведенные формулы (чтобы совсем не трогать код, "вытащил" из моих констант эти кси и омега): % s^2 + 31.1 s + 799.4 % --------------------- % s^2 + 7.109 s + 25.27 clear all M_PI = pi om1 = sqrt(25.27); om0= sqrt(799.4) ksi0 = 31.1 / (2*om0); ksi1 = 7.109 / (2*om1); fs = 1000; % sample rate = 1000 SPS Tsam=1/fs; tau0=.5/tan(Tsam*om0/2.); tau1=.5/tan(Tsam*om1/2.); d=1.+(4.*tau1*(tau1+ksi1)); b0=(1.+(4.*tau0*(tau0+ksi0)))/d; b1=2.*(1.-4.*tau0*tau0)/d; b2=(1.+(4.*tau0*(tau0-ksi0)))/d; a1=2.*(1.-4.*tau1*tau1)/d; a2=(1.+(4.*tau1*(tau1-ksi1)))/d; a0=1; то имею коэффициенты: Quote bd = 0.0320 -0.0630 0.0310 ad = 1.0000 -1.9929 0.9929 вижу [bd], уменьшенный примерно в 31 раз если сравнить с результатом функции матлаба. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Самурай 12 19 марта, 2019 Опубликовано 19 марта, 2019 (изменено) · Жалоба Все дело в том, что в формулах ТС есть некоторые упрощения, которые не влияют на результат в частном случае, когда w0 = w1. Для общего случая, формулы нужно слегка поправить : double om0=2.*M_PI*f0; double om1=2.*M_PI*f1; double tau0=.5/tan(Tsam*om0/2.); double tau1=.5/tan(Tsam*om1/2.); double d=1.+(4.*tau1*(tau1+ksi1));double d0=1.+(4.*tau0*(tau0+ksi0)); double k=(d0/d)*(tau1/tau0)^2; double b0=k*(1.+(4.*tau0*(tau0+ksi0)))/d0; double b1=k*2.*(1.-4.*tau0*tau0)/d0; double b2=k*(1.+(4.*tau0*(tau0-ksi0)))/d0; double a1=2.*(1.-4.*tau1*tau1)/d; double a2=(1.+(4.*tau1*(tau1-ksi1)))/d; a0=1; Изменено 19 марта, 2019 пользователем Самурай Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Ruslan1 16 20 марта, 2019 Опубликовано 20 марта, 2019 · Жалоба 22 hours ago, Самурай said: Все дело в том, что в формулах ТС есть некоторые упрощения, которые не влияют на результат в частном случае, когда w0 = w1. Для общего случая, формулы нужно слегка поправить : ... Уважаемый Самурай, все получилось! Спасибо огромное!! С помощью Вашей подсказки и Википедии получилось сделать универсальный вариант, с любыми нумераторами-деноминаторами 2-й степени (в стартовом коде было a0=b0=1). Оставлю тут , может пригодится кому-нибудь: % second-order bilinear transform % input: % Ha(s) = ( b0*s^2 + b1*s + b2 ) / (a0*s^2 + a1*s + a2) % % output: (digital biquad filter, a0 = 1) % Hd(z) = (b0 + b1*z^-1 + b2*z^-2) / (a0 + a1*z^-1 + a2* z^-2) % % reference: % #1 (General second-order biquad transformation): https://en.wikipedia.org/wiki/Bilinear_transform % #2 (correction for case w0 != w1): https://electronix.ru/forum/index.php?app=forums&module=forums&controller=topic&id=148759 clear all %% input data b_in = [23 31.1 799.4]; a_in = [19 7.109 25.27]; fs = 1000; %% read data from input arrays b0 = b_in(1); b1 = b_in(2); b2 = b_in(3); a0 = a_in(1); a1 = a_in(2); a2 = a_in(3); %% calculation w0 = sqrt(b2); w1= sqrt(a2) T = 1/fs; K0 = w0/tan(w0*T/2); K1 = w1/tan(w1*T/2); %not normalled coefficients bn0 = b0*K0^2 + b1*K0+b2; bn1 = 2*b2-2*b0*K0^2; bn2 = b0*K0^2 - b1*K0 + b2; an0 = a0*K1^2 + a1*K1 + a2; an1 = 2*a2 - 2*a0*K1^2; an2 = a0*K1^2 - a1*K1 + a2; % correction for denominator k = (K1/K0)^2; %k=1 % Normally the constant term in the denominator must be normalized to 1 bd0 = k* bn0 / an0; bd1 = k* bn1 / an0; bd2 = k* bn2 / an0; ad0 = an0 / an0; ad1 = an1 / an0; ad2 = an2 / an0; %% output data bd = [bd0 bd1 bd2] ad = [ad0 ad1 ad2] %% test [bd_matlab,ad_matlab] = bilinear(b_in,a_in,fs); ad ad_matlab bd bd_matlab fvtool (bd_matlab,ad_matlab); fvtool (bd,ad) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться