Цифровой фильтр 2 порядка

[destination_sirius]

Строю цифровой фильтр 2-го порядка по этому руководству: http://yehar.com/blog/?p=121 (начиная с заголовка "From poles and zeros to filter coefficients")

Уравнение его такое:

$$out(t) = a_0in(t) + a_1in(t-1) + a_2in(t-2) +b_1out(t-1) + b_2out(t-2) \eqno(1)$$

где in(t) и out(t) - входные и выходные сигналы соответственно. Автор пишет, что любой фильтр 2-го порядка можно записать таким образом.

Передаточная функция: $H(z)=\frac{a_0(z-z_1)(z-z_2)}{(z-p_1)(z-p_2)}$,

где z - частота в виде комплексного числа.

$z_1=z_x+iz_y, z_2=z_x-iz_y$ - позиции двух сопряжённых "нулей", где фильтр подавляет частоты
$p_1=p_x+ip_y, p_2=p_x-ip_y$ - позиции двух сопряжённых "пиков" где фильтр подчёркивает частоты.

Мне нужен фильтр низких частот. Он имеет следующую схему:
($2\pi$ - частота дискретизации, соответственно, $\pi$ - частота Найквиста)

             __---|---__        два совпадающих нуля на (-1;0)
           _/     |     \_
          /       |    x  \     Пик и его сопряжённая пара на:
         |        |        |    (px=r*cos(2pi f/SR); py=r*sin(2pi f/SR)), где r-расстояние пика от центра единичного круга, f - частота среза фильтра
        -o--------+--------|-
         |        |        |    использую в нормализации 0 Гц.
          \_      |    x _/
            \__   |   __/
               ---|---

Значит, $z_1=z_2=-1+i\cdot0=-1$
Подставив и преобразовав, получил такое:

$H(z)=\frac{a_0(z-(-1))^2}{(z-p_1)(z-p_2)}=\frac{a_0(z+1)^2}{(z^2-p_1z-p_2z+p_1p_2)}=\frac{a_0(z+1)^2}{z^2 + (-(p_x + ip_y) -(p_x - ip_y))z + (p_x - ip_y)(p_x + ip_y)}$ $=\frac{a_0z^2+2a_0z+a_0}{z^2-2p_xz + p_x^2 + p_y^2}$

Автор упомянутого блога в статье разбирает похожий пример, только он делает фильтр без "нулей", и в числителе у него только $a_0$. Он предлагает дальше сделать так:

"Давайте посмотрим на знаменатель.

Пусть $z^2-2p_xz+p_x^2+p_y^2=0$ (помножим на $z^{-2}$)

$z^0-2p_xz^{-1}+(p_x^2+p_y^2)z^{-2}=0$

$z^0=2p_xz^{-1}-(p_x^2+p_y^2)z^{-2}$

Коэффициенты при разных степенях z - это на самом деле коэффициенты к выходу фильтра:

$out(t-0)=2p_xout(t-1)-(p_x^2+p_y^2)out(t-2)$"

Вопрос 1: откуда он это взял? (его фильтр работает, ошибок нет, я потестил)

Но у меня ещё есть чиститель. По аналогии я беру числитель: $a_0z^2+2a_0z+a_0$, умножаю на $z^{-2}$:
$a_0z^0+2a_0z^{-1}+a_0z^{-2}$.
Значит, коэффициенты ко входным значениям будут: $a_0in(t)$, $2a_0in(t-1)$, $a_0in(t-2)$, и я могу переписать уравнение (1) так:

$out(t)=a_0in(t)+2a_0in(t-1)+a_0in(t-2)+2p_xout(t-1)-(p_x^2+p_y^2)out(t-2)$ ?

Теперь надо сделать нормализацию, т.е. найти $a_0$. Для фильтра низких частот она вроде как делается по 0Гц, то есть $Z=(1;0)$

$H(z)=\frac{a_0(z-z_1)(z-z_2)}{(z-p_1)(z-p_2)}=\frac{a_0(1-(-1))^2}{(1-p_1)(1-p_2)}=\frac{4a_0}{(1-2p_x+p_x^2 + p_y^2)}=1$
$a_0=0.25(1-2p_x+p_x^2 + p_y^2)$
Поскольку $p_x=rcos(\frac{2\pi f}{F_s})$, $p_y=rsin(\frac{2\pi f}{F_s})$, $a_0=0.25(1-rcos(\frac{2\pi f}{F_s})+r^2)$.
Выше мы нашли:
$a_1=2a_0,$
$a_2=a_0,$
$b_1=2p_x=2rcos(\frac{2\pi f}{F_s})$
$b_2=-(p_x^2+p_y^2)=-r^2$

Подставляем в наше уравнение:
out(t) = a_0in(t) + a_1in(t-1) + a_2in(t-2) +b_1out(t-1) + b_2out(t-2)$

Вопрос 2: Нифига не работает. Кто видит ошибку?

[profrotter]

destination_sirius в сообщении #868623 писал(а):

откуда он это взял?

Обозначим $x(n)$ - сигнал на входе цифрового фильтра, $y(n)$-сигнал на выходе. Разностное уравнение цифрового фильтра 2-го порядка имеет вид:

$$y(n)=a_0x(n)+a_1x(n-1)+a_2x(n-2)+b_1y(n-1)+b_2y(n-2),$$

а выражение для системной функции:

$$H(z)=\frac{a_0+a_1z^{-1}+a_2z^{-2}}{1-b_1z^{-1}-b_2z^{-2}}.$$

Если выражение для системной функции приведено к указанному (стандартному) виду, то коэффициенты разностного уравнения находятся при соответствующих степенях $z^{-1}$. Чтобы получить этот результат надо рассмотреть Z-преобразование от разностного уравнения при нулевых начальных условиях и учесть, что системная функция представляет собой отношение изображения сигнала на выходе фильтра к изображению сигнала на входе при нулевых начальных условиях.

destination_sirius в сообщении #868623 писал(а):

Нифига не работает

Никто, кроме вас, не знает что такое "не работает" и как Вы хотите чтобы это работало.

[svv]

Попробую ортогонально дополнить.

destination_sirius в сообщении #868623 писал(а):

Коэффициенты при разных степенях z - это на самом деле коэффициенты к выходу фильтра:
$out(t-0)=2p_xout(t-1)-(p_x^2+p_y^2)out(t-2)$"
Вопрос 1: откуда он это взял?

Понятно, что в таком фильтре два любых последовательных значения $y_{n-1}$ и $y_{n-2}$ определяют все последующие, поэтому с практической точки зрения он малополезен. А с теоретической он описывает систему, в которой, начиная с некоторого $y_n$, возбуждающий (входной) сигнал $x_n$ равен нулю, и система долго (до бесконечности) «звенит» на собственных частотах. Вопрос: на каких именно?

Пусть значение $y_n$ соответствует моменту времени $t_n=n\tau$. В случае непрерывного сигнала одной «частоте» $\omega$, которая может быть и комплексной, соответствовал бы сигнал $y(t)=e^{i\omega t}$. Дискретным аналогом будет
$y_n=e^{i\omega t_n}=e^{i\omega \tau n}=z^n$, где $z=e^{i\omega \tau}$.

Подставим это в разностное уравнение $y_n+a_1y_{n-1}+a_2y_{n-2}=0$, описывающее систему. Получим после сокращения квадратное уравнение относительно $z$:
$z^2+a_1 z+a_2=0$
Важный момент: итак, частоты собственных колебаний не могут быть любыми, они определяются подобными уравнениями.
Обозначим корни через $p_1$ и $p_2$ (обычно рассматриваются корни полинома от $z^{-1}$, но я подстраиваюсь под Ваш текст). Чтобы в системе был возможен вещественный сигнал, корни должны быть либо вещественными, либо образовывать пары сопряженных корней. Выбираем последний вариант исходя из цели построения фильтра. Тогда
$p_1=p_x+ip_y$
$p_2=p_x-ip_y$
По теореме Виета
$p_1+p_2=2p_x=-a_1$
$p_1 p_2=p_x^2+p_y^2=a_2$
Следовательно,
$z^2-2p_x z+p_x^2+p_y^2=0$

Как видите, всё логично, если рассуждать в обратном направлении, от разностного уравнения к уравнению для корней полинома в знаменателе.

[destination_sirius]

profrotter в сообщении #868690 писал(а):

destination_sirius в сообщении #868623 писал(а):

Нифига не работает

Никто, кроме вас, не знает что такое "не работает" и как Вы хотите чтобы это работало.

Пускаю белый шум, регулирую частоту, и вижу, что он не является фильтром низких частот. Может где-то ошибка у меня в логике рассуждений?
Если сделать как у автора в статье: убрать числитель с "нулями", превратив его в band-pass filter, оставив в числителе только $a_0$, и найти его исходя из нормализации, что H(z) на частоте "пика" равна 1, то всё ок, то есть сам метод верный, я его не так понял, и хотелось бы понять, где именно...

[profrotter]

destination_sirius в сообщении #868749 писал(а):

убрать числитель с "нулями",

Ваш фильтр можно рассматривать как последовательное соединение трёх звеньев. Числитель ваш аккурат соответствует ФНЧ, который представляет собой два последовательных звена фильтра скользящего среднего. Системная функция фильтра скользящего среднего первого порядка даётся выражением $H_{MAF}(z)=\frac{1}{2}(1+z^{-1})$. Два последовательно соединённых звена таких фильтров будут иметь системную функцию $H_{MAF2}=\frac{1}{4}(1+z^{-1})^2$. Это гарантировано ФНЧ. А дальше последовательно цифровой фильтр второго порядка с системной функцией $H_2(z)=\frac{4a_0}{1-2p_xz^{-1}+r^2}$. Он разобран в примере, на который Вы ссылаетесь. В зависимости от того, каким Вы выбираете модуль $r$ и аргумент $\varphi=\omega_0T$ ($T$ - период дискретизации, $\omega_0$ - параметр) полюса системной функции $p=p_x+ip_y=re^{i\varphi}$, можно получить полосовой фильтр, фильтр нижних частот или фильтр верхних частот. В случае комплексно-сопряжённых полюсов системной функции, как у вас, фильтр нижних частот получается, если ориентировочно $\omega_0T<\frac{\pi}{4}$ и $r<\frac{1}{2}$. В общем при малых $r$.

Для того, чтобы понять, что за фильтр Вы получили следует исследовать его амплитудно-частотную характеристику, а не пускать белый шум.