В статье «Цифровой фильтр 2 порядка» есть упоминание о том, как спроектировать RC фильтр любого порядка и указана передаточная функция \frac{Y(s)}{X(s)} = \frac{1}{(1+s\tau )^{n}}, где n — порядок фильтра. Читателю будет полезно знать, как выглядит АЧХ (амплитудно-частотная характеристика) того, или иного фильтра. В данный момент, мне уже намного легче объяснить материал, так как на блоге имеются некоторые наработки и существует возможность направлять читателя по ссылкам, детально описывающим частные вещи. Так вот, как выглядит АЧХ RC фильтра первого порядка указано тут или тут. Зная выражение, описывающее АЧХ RC фильтра первого порядка, несложно записать АЧХ для RC фильтра любого порядка:

|K(j\omega)| = (\frac{1}{\sqrt{1+({\omega}RC})^2})^n = (\frac{1}{\sqrt{1+({\omega}\tau})^2})^n

В то время как выражение для фильтра первого порядка:

\tau = RC = \frac{1}{\omega_{cutoff}}

верно (выражение взято из статьи «Трансформация частотной оси и создание полноценного БИХ фильтра«), то для фильтра второго порядка и старше нужно придумать что-то другое. Напомню, идея создания фильтров старшего порядка состоит в том, что мы последовательно, друг за другом включаем RC фильтры первого порядка. А теперь давайте подумаем. Если один фильтр первого порядка на частоте среза будет уменьшать амплитуду, скажем, на 3 дБ (это приблизительно в 1.41 раза), то если включить 2 таких фильтров один за другим, вместе они будут уменьшать амплитуду на данной частоте среза уже на 6 дБ (это приблизительно в 2 раза)! Именно такими свойствами и обладает RC фильтр второго порядка, описанный в статье «Цифровой фильтр 2 порядка«. Скажу прямо — это не очень хорошо! Обычно, на частоте среза фильтра принято иметь ослабление уровня -3 дБ. В нашем случае, добиться именно такого результата проще простого. Давайте ещё раз посмотрим на формулу передаточной функции:

\frac{Y(s)}{X(s)} = \frac{1}{(1+s\tau )^{n}}

В этом выражении мы видим \tau, но чему оно равно мы не знаем. Это установить совсем не сложно. Именно это константное значение и определяет такие параметры фильтра как частота среза, и уровень ослабления на этой частоте. Вывести значение \tau можно из формулы:

|K(j\omega)| = (\frac{1}{\sqrt{1+({\omega}\tau})^2})^n

Стоит напомнить, что для нас значит данное выражение. При постоянном значении \tau мы изменяем частоту \omega и смотрим на коэффициент передачи системы (получившееся число). На каких-то частотах мы будем получать значения, близкие к единице (фильтр пропускает данную частоту), на каких-то — близкие к нулю (фильтр подавляет данную частоту), а на всех остальных — значения между 0 и 1 :-) (фильтр частично подавляет частоту). Как я уже писал, обычно, фильтры разрабатывают для частоты среза, которая составляет -3 дБ (\frac{1}{\sqrt{2}}=0.707 от амплитуды входной гармонической составляющей, например, синусоиды). Что такое децибел — читаем тут. В таком случае, для частоты среза \omega_{cutoff} можно записать:

\frac{1}{\sqrt{2}} = \frac{1}{\sqrt{1+({\omega_{cutoff}}\tau})^n}

Несложно переписать выражения относительно \tau:

\tau = \frac{\sqrt{2^{\frac{1}{n}}-1}}{\omega_{cutoff}}

Все карты на руках! Мы получили желаемую зависимость \tau от частоты среза \omega_{cutoff} для RC фильтра любого порядка. Переходим непосредственно к отображению АЧХ на графике. Существует несколько путей для отображения АЧХ. Если бы мы были ограничены математикой без комплексных чисел, мы бы воспользовались выражением:

|K(j\omega)| = (\frac{1}{\sqrt{1+({\omega}\tau})^2})^n

Но это очень банально. Вы должны уметь получать АЧХ сразу из передаточной функции:

\frac{Y(s)}{X(s)} = \frac{1}{(1+s\tau )^{n}}

Как его получить? Просто использовать подстановку s = j\omega:

\frac{Y(s)}{X(s)} = \frac{1}{(1+j\omega\tau )^{n}}

Python имеет в своём арсенале такую мощную библиотеку как NumPy, которая прекрасно умеет работать с комплексными числами. В результате, имеем несколько строк кода для отображения АЧХ фильтров разных порядков:

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
# -*- coding: utf-8 -*- 
""" 
Created on Tue Jul 02 21:52:15 2015 
@author:  Viktor Signaievskyi 
@blog:    www.HowToDoIt.com.ua 
@e-mail:  piligrim2007@meta.ua 
@details: Example program for evaluation and plotting of LPFs' frequency   
          response for different order filters
Copyrighted by ©Viktor Signaievskyi 
""" 
import numpy as np 
import matplotlib.pyplot as plot 
 
Fs = 200.0 # Sampling frequency (Hz)
Fc = 1.0 # Cutoff frequency (Hz) for -3 dB level
 
MAX_LPF_ORDER = 5
 
omega = np.arange(0, np.pi, 0.01)
omega_cutoff = 2.0*np.pi*Fc/Fs # Angular cutoff frequency
 
plot.clf()
 
# Calculating frequency response for LPFs of different order
for n in range(1, MAX_LPF_ORDER+1):
    tau = np.sqrt(np.power(2.0, 1.0/n) - 1.0) / omega_cutoff
    K = np.power(1.0/(1.0+(1j*omega)*tau), n)
    plot.plot((Fs/2)*omega/np.pi, 20*np.log10(np.abs(K)[:]))
    plot.xscale('log')
 
plot.xlabel("Frequency, Hz")
plot.ylabel("|K(jw)|, dB")
 
plot.draw()

Ниже, результат отработки скрипта в виде графика:

Frequency_Responses

Цифрами от 1 до 5 указан порядок фильтра. Мини-анализ окончен!

Удачи в познаниях мира ЦОС!