digital-filters — библиотека синтеза и анализа цифровых фильтров

Python-библиотека для синтеза и анализа цифровых фильтров: реализует три классических метода — метод инвариантной импульсной характеристики, прямой синтез по полюсам системной функции и билинейное z-преобразование с частотными предыскажениями (prewarping) — поверх NumPy, с автоматическим выбором порядка, проверкой устойчивости и диагностической визуализацией для каждого получившегося фильтра.

Методы и инструменты: метод инвариантной импульсной характеристики · прямой синтез по полюсам системной функции · билинейное z-преобразование, prewarping · разностные уравнения, структуры FIR/IIR (трансверсальная/каноническая) · анализ в z-плоскости: нули, полюса, критерий устойчивости |z|<1 · синтез фильтра Баттерворта произвольного порядка · NumPy (свёртка полиномов, np.roots) · SciPy (freqz, lfilter) · Python (dataclasses, типизация) · Matplotlib

Highlights

• Три метода синтеза — модифицированный метод инвариантной импульсной характеристики (RC-цепи), прямой синтез по полюсам (резонатор), билинейное z-преобразование с prewarping (ФНЧ Баттерворта) — реализованы на NumPy, без вызовов scipy.signal.butter/iirdesign. SciPy используется отдельно, ниже по пайплайну, только для анализа уже рассчитанных коэффициентов (freqz, lfilter).
• В методических пособиях явно приведена формула передаточной функции фильтра Баттерворта только для порядка N=2 (пример «расчёт фильтра второго порядка»). В butterworth.py та же схема расчёта — полюса аналогового прототипа → раскрытие полинома по этим полюсам → билинейная подстановка — записана для произвольного N, а сам порядок вычисляется автоматически по частоте среза $f_0$, контрольной частоте $f_1$ и требуемому затуханию в дБ. При N=2 обобщённая формула численно сходится к явному виду из методических указаний.
• Устойчивость каждого IIR-фильтра проверяется вычислением корней знаменателя (np.roots) и условием $|z_i|<1$ — включая цифровой резонатор, где устойчивость и так следует из метода синтеза (полюса задаются внутри единичной окружности по построению), но всё равно проходит через ту же функцию is_stable(), что и остальные фильтры.
• Переход FIR → IIR — не ручное решение по каждому конкретному τ, а следствие порогового условия на затухание импульсной характеристики (5% от первого отсчёта); переключение проверено по обе стороны границы N=10 на нескольких значениях τ/T для дифференциатора и интегратора.
• Коэффициенты, полюса и результат проверки устойчивости для всех реализованных фильтров сверены с расчётом варианта 20 (МГТУ им. Н.Э. Баумана, курс «Цифровая обработка сигналов», преп. Г.В. Круглов) — совпадение с ручными вычислениями отчёта до 4–7 значащих цифр (таблица ниже).

Оглавление

• Происхождение и цель проекта
• Архитектура репозитория
• Теоретическая база
• Результаты и верификация
• Инженерные решения
• Быстрый старт
• Возможные улучшения
• См. также
• Литература

Происхождение и цель проекта

Репозиторий вырос из необходимости расчёта нескольких классов цифровых фильтров разными методами синтеза — ЦФ, эквивалентный дифференцирующей RC-цепи; ЦФ, эквивалентный интегрирующей RC-цепи; цифровой резонатор; ФНЧ Баттерворта, для дальнейшей проверки, при помощи макетирования, расчитанных фильтров.

Расчёт оформлен не набором одноразовых скриптов под конкретное ТЗ, а небольшой типизированной библиотекой: единые модели данных (FIRFilter/IIRFilter), единая функция проверки устойчивости, единый модуль визуализации — так, что любой новый набор входных параметров ($\tau$, $f_0$, $f_1$, $B$, $L$, $F_д$) проходит через уже проверенный код, а не через копию скрипта с изменёнными числами.

Архитектура репозитория

digital_filters/
├── digital_filters/
│   ├── fir/
│   │   ├── __init__.py
│   │   └── rc.py                  # differentiating_rc_fir, integrating_rc_fir
│   ├── iir/
│   │   ├── __init__.py
│   │   ├── rc.py                  # differentiating_rc_iir, integrating_rc_iir
│   │   ├── resonator.py           # resonator_iir (прямой синтез по полюсам)
│   │   └── butterworth.py         # butterworth_lowpass (билинейное преобразование, произвольный N)
│   ├── modeles.py                 # dataclass-модели FIRFilter / IIRFilter
│   ├── analysis.py                # is_stable() — проверка полюсов по |z|<1
│   ├── plot.py                    # ИХ / АЧХ / карта нулей-полюсов, сводная панель
│   └── __init__.py
├── images/                    # сохранённые графики по каждому фильтру
└── Test_filters_calc_and_plot.ipynb   # расчёт и визуализация для всех фильтров

FIRFilter/IIRFilter — @dataclass(slots=True)-модели: коэффициенты числителя/знаменателя, метаданные синтеза (design), и вычисляемые (не хранимые) свойства zeros/poles — is_stable() всегда отражает текущие коэффициенты, а не сохранённое ранее значение.

Теоретическая база

1. ЦФ, эквивалентный дифференцирующей RC-цепи

Дифференцирующая RC-цепь — фильтр верхних частот, поэтому прямой метод инвариантной импульсной характеристики (ИХ) здесь не работает: полоса пропускания выходит выше половины частоты дискретизации, и АЧХ ЦФ искажается наложением спектров (алиасинг складывает АЧХ дифференциатора с постоянной составляющей). Используется модификация метода — через переходную характеристику (ПХ) $g(t)=e^{-t/\tau}$; ИХ ЦФ получается как разность отсчётов ПХ, сдвинутых на один такт:

$$h(kT) = g(kT) - g((k-1)T) = e^{-kT/\tau}\left(1 - e^{T/\tau}\right)$$

Коэффициенты нерекурсивного (FIR) фильтра: $a_0=1,\ a_k = e^{-k/(\tau/T)}\left(1-e^{1/(\tau/T)}\right)$ для $k\ge1$.

Число звеньев N определяется условием $\forall k\ge N:\ |a_k/a_1| < 0.05$; при N>10 — переход к рекурсивной (IIR) форме:

$$H(z)=\frac{1-z^{-1}}{1-e^{-T/\tau}z^{-1}} \quad\Rightarrow\quad a_0=1,\ a_1=-1,\ b_1=e^{-T/\tau}$$

Так как $0 < e^{-T/\tau} < 1$ при любых $\tau, T>0$, рекурсивный дифференциатор устойчив при любых входных параметрах — полюс всегда строго внутри единичной окружности.

2. ЦФ, эквивалентный интегрирующей RC-цепи

Та же схема расчёта, симметрично: $g(t)=1-e^{-t/\tau}$,

$$h(kT)=e^{-kT/\tau}\left(e^{T/\tau}-1\right),\ k>0;\qquad h(0)=0$$

FIR: $a_0=0,\ a_k = e^{-k/(\tau/T)}\left(e^{1/(\tau/T)}-1\right)$. При N>10 — IIR-форма:

$$H(z)=\frac{\left(1-e^{-T/\tau}\right)z^{-1}}{1-e^{-T/\tau}z^{-1}} \quad\Rightarrow\quad a_0=0,\ a_1=1-e^{-T/\tau},\ b_1=e^{-T/\tau}$$

Устойчивость гарантирована по той же причине, что и для дифференциатора.

3. Цифровой резонатор — прямой синтез по полюсам

Здесь используется не метод по прототипу, а прямой синтез: полюса системной функции задаются напрямую в нужной точке z-плоскости —

$$s = e^{-\alpha-j\omega_0},\quad s^*=e^{-\alpha+j\omega_0},\qquad \omega_0=\frac{2\pi f_0}{f_д},\quad \alpha=\frac{\pi B}{f_д}$$

$$H(z)=\frac{1}{(1-sz^{-1})(1-s^*z^{-1})}=\frac{1}{1-2e^{-\alpha}\cos\omega_0,z^{-1}+e^{-2\alpha}z^{-2}}$$

$$a_0=1,\qquad b_1=2e^{-\alpha}\cos\omega_0,\qquad b_2=-e^{-2\alpha}$$

Отдельная проверка устойчивости не требуется: полюса изначально размещаются строго внутри единичной окружности при любой полосе $B>0$ — устойчивость встроена в сам метод синтеза.

4. ФНЧ Баттерворта — билинейное z-преобразование

Аналоговый прототип задаётся не ИХ, а АЧХ, поэтому переход к ЦФ выполняется билинейным z-преобразованием, а не методом инвариантной ИХ:

$$|K(\omega)|^2=\frac{1}{1+(\omega/\omega_0)^{2N}}$$

Порядок фильтра по заданным $f_0$, $f_1$ и затуханию $L_{дБ}$ на частоте $f_1$:

$$N \ge \frac{\lg(2L_{раз}-1)}{2\lg(f_1/f_0)}, \qquad L_{раз}=10^{L_{дБ}/20}$$

Полюса нормированного прототипа: $s_i = \exp!\left[-j\pi\left(\tfrac12+\tfrac{2i-1}{2N}\right)\right],\ i=1..N$.

Частота среза «искажается» под билинейное преобразование (prewarping), чтобы аналоговая и цифровая частоты среза совпали после подстановки:

$$\omega_{0,1}^{бп}=2f_д\tan!\left(\frac{\pi f_{0,1}}{f_д}\right)$$

Методичка курса раскрывает передаточную функцию только для N=2. В butterworth.py та же схема — с тем же prewarping и той же билинейной подстановкой $p=\frac{2}{T}\cdot\frac{z-1}{z+1}$ — записана в общем виде для произвольного порядка N:

$$H(z)=\frac{(z+1)^N}{\displaystyle\sum_{n=0}^{N} c_n, k^n,(z-1)^n(z+1)^{N-n}}, \qquad k=\frac{2f_д}{\omega_0^{бп}}$$

где $c_n$ — коэффициенты полинома $D(x)=\prod_{i=1}^N(x-s_i)=\sum_n c_n x^n$, полученного раскрытием произведения по полюсам прототипа (np.poly). При N=2 выражение численно сворачивается ровно к явной формуле из методички — этим соответствием оно и проверялось.

Результаты и верификация

Коэффициенты, полюса/нули и проверка устойчивости, полученные библиотекой для ($F_д=10\text{ кГц}$, $\tau_{диф}/T=1.28$, $\tau_{инт}/T=20.5$, $f_0=1520\text{ Гц}$, $B=136\text{ Гц}$, $f_1=2f_0=3040\text{ Гц}$, $L=15\text{ дБ}$), сверенные с ручным расчётом:

Фильтр	Метод	N	Реализация	Коэффициенты	Устойчив
Дифференцирующая RC (τ/T=1.28)	Модиф. инвариантная ИХ	5	FIR (трансверсальная)	a=[1, -0.5422, -0.2482, -0.1136, -0.0520, -0.0238]	не требуется (FIR)
Интегрирующая RC (τ/T=20.5)	Модиф. инвариантная ИХ	63→IIR	IIR (каноническая, 1 порядок)	a=[0, 0.04761], b=[1, -0.95239]	$b_1=0.9524<1$
Цифровой резонатор (f0=1520, B=136)	Прямой синтез по полюсам	2 (фикс.)	IIR (каноническая)	b=[1, -1.10683, 0.91810]	полюса $0.5534\pm0.7822j$
ФНЧ Баттерворта (f0=1520, f1=3040, L=15дБ)	Билинейное преобразование + prewarping	2 (автовыбор)	IIR (каноническая)	a=[0.13391, 0.26783, 0.13391], b=[1, -0.73238, 0.26804]	полюса $0.3662\pm0.3660j$, $|z|=0.518<1$

Значения коэффициентов дифференцирующего/интегрирующего FIR, b₁/b₂ резонатора и всех пяти коэффициентов Баттерворта совпадают с ручным расчётом отчёта до 4–7 значащих цифр.

Графики (импульсная характеристика · АЧХ · карта нулей-полюсов для IIR) — в digital_filters/images/:

Структурные схемы (трансверсальная/каноническая), построенные по рассчитанным коэффициентам:

За пределами варианта 20. Каждый метод дополнительно прогонялся на других наборах параметров:

• Дифференциатор/интегратор — на нескольких соотношениях τ/T по разные стороны порога N=10: τ/T=1.28→N=5 (FIR), τ/T=3.4→N=12 (IIR), τ/T=2.6→N=9 (FIR), τ/T=20.5→N=63 (IIR).
• Синтез Баттерворта — на параметрах, не совпадающих с вариантом 20 (например, $f_0=40\text{ кГц}$, $F_д=400\text{ кГц}$, $L=22\text{ дБ}$), что даёт порядок N, отличный от 2, и проверяет обобщённую формулу за пределами частного случая из методички.

Инженерные решения

• Расчёт коэффициентов — на NumPy, без scipy.signal.butter/bilinear/iirdesign. SciPy применяется ниже по пайплайну, только для freqz/lfilter — анализ уже готовых коэффициентов, не их расчёт.
• Первая версия differentiating_rc_fir использовала ту же формулу масштабирования, что и integrating_rc_fir, отличаясь только начальным коэффициентом $a_0$. Это давало АЧХ дифференциатора с тем же спадающим профилем, что у интегратора, вместо ожидаемого роста усиления к высоким частотам (нулевой коэффициент передачи на постоянном токе — характерное свойство дифференцирующей цепи). Расхождение стало видно при сравнении графика АЧХ с этим свойством и было исправлено сменой знака в масштабирующем множителе.
• В модели a — числитель, b — знаменатель, что обратно порядку имён в scipy.signal (там числитель передаётся первым позиционным аргументом, но называется b). Во всех вызовах freqz/lfilter порядок аргументов прокомментирован явно.
• is_stable() в analysis.py — общая точка проверки для всех фильтров: даже для резонатора, где устойчивость следует из метода синтеза, она пересчитывается по факту, а не принимается по умолчанию.

Быстрый старт

git clone <URL_репозитория>
cd digital_filters
pip install numpy scipy matplotlib
jupyter notebook Test_filters_calc_and_plot.ipynb

Пример использования библиотеки напрямую:

from digital_filters.fir import *
from digital_filters.iir import *

# RC fir дифференциатор
tau_T = 1.28
my_RC_filter_1 = differentiating_rc_fir(tau_T)

# ФНЧ Баттерворта, вариант 20: автоматический расчёт порядка по f0/f1/L
F_d = 10000
f_0 = 1520
f_1 =  2 * f_0
L_dB = 15
lp = butterworth_lowpass(f_0, F_d, f1 = f_1, attenuation_db = L_dB)
print(lp.a, lp.b, is_stable(lp))
plot_filter(lp, fd=10000)   # ИХ + АЧХ + карта нулей-полюсов на одной панели

# Цифровой резонатор: прямой синтез по полюсам
res = resonator_iir(center_frequency=1520, bandwidth=136, sample_rate=10000)
plot_filter(res, fd=10000)

Возможные улучшения

• Обобщить пайплайн Баттерворта (прототип → prewarping → билинейная подстановка) на фильтры верхних частот, полосовые и режекторные — три этапа уже разделены как отдельные функции, добавление нового типа не требует переписывать билинейное преобразование.
• Добавить прототипы Чебышева/эллиптический как альтернативные poles_prototype в тот же пайплайн билинейного преобразования.
• Формализовать текущую ручную кросс-проверку с вариантом 20 в pytest-тесты с эталонными значениями коэффициентов.
• Сравнить с scipy.signal.butter/bilinear как эталоном для произвольного порядка N.
• Добавить в plot.py панель ФЧХ/групповой задержки в дополнение к ИХ/АЧХ/нулям-полюсам.

См. также

ofdm-phy-simulation — полностью смоделированный в GNU Radio OFDM-приёмопередатчик, где decimating LPF после квадратурного смесителя (тот же класс задачи, что и ФНЧ здесь) убирает image-частоту и LO-leakage перед демодуляцией, а весь тракт TX → AWGN-канал → RX измерен по packet delivery rate на 30 уровнях шума.

Литература

Круглов Г.В. К выполнению домашнего задания и лабораторной работы по курсу «Цифровая обработка сигналов». — М.: МГТУ им. Н.Э. Баумана, 2015.

Лицензия

MIT — см. LICENSE.