Ранее, я рассказывал Вам о том, как получить простой цифровой фильтр нижних частот на основе аналогового RC фильтра. Полученные выражения обладали определёнными свойствами. Как помните, некоторые частоты он хорошо пропускал, некоторые — хорошо подавлял. Если посмотреть на реализацию примера фильтра, написанного мною на Питоне, то можно заметить, что в результате получается не совсем чистая синусоида. Т. е. в этом конкретном случае уровень фильтрации был недостаточным. Что же сделать для увеличения эффекта фильтрации? Существуют разные варианты. Первое, что приходит на ум, это пропустить уже отфильтрованный сигнал через тот же самый фильтр ещё раз. Существует и более элегантный способ. Для этого нам с вами придется изобрести фильтр высшего порядка. В основе данного фильтра как раз и лежит идея включения того же самого элементарного фильтра несколько раз подряд. Но есть один нюанс. Эти действия необходимо выполнить в математической форме с использованием преобразования Лапласа (т. е. в частотной области) и обратного Z-преобразования. В результате получим элементарную формулу для фильтрации сигнала вида «что на что умножить и что с чем просуммировать чтобы получить выходной, отфильтрованный сигнал».

Промежуточным результатом вводной статьи по цифровой фильтрации было выражение, описывающее фильтр нижних частот первого порядка (в частотной области коплексной переменной s):

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

Напомню, для того чтобы отфильтровать входной сигнал в частотной области, необходимо умножить его на передаточную функцию (в частотной области, коплексной переменной s:)). Следовательно, чтобы отфильтровать результатирующий сигнал ещё раз, нужно опять же умножить его на ту же передаточную функцию. Т. е. передаточной функцией фильтра 2-го порядка будет возведённая в степень 2 передаточная функции фильра 1-го порядка:

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

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

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

Настал черёд преобразовать передаточную функцию непрерывного аргумента s в передаточную функцию дискретного аргумента z. Для этого, воспользуемся Билинейным преобразованием:

s=\frac{2}{T}\frac{z-1}{z+1}

Используем данную замену для выражения, описывающего фильтр 2-го порядка:

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

Выполняем подстановку:

\frac{Y(z)}{X(z)}=\frac{1}{1+2(\frac{2}{T}\frac{z-1}{z+1})\tau+((\frac{2}{T}\frac{z-1}{z+1})\tau)^2}=\newline =\frac{(T(z+1))^2}{T^2(z+1)^2+4T\tau(z-1)(z+1)+4\tau^2(z-1)^2}=\newline =\frac{T^2(z^2+2z+1)}{T^2(z+1)^2+4T\tau(z^2-1)+4\tau^2(z^2-2z+1)}=\newline =\frac{T^2(z^2+2z+1)}{T^2(z+1)^2+4T\tau(z^2-1)+4\tau^2(z^2-2z+1)}=\newline =\frac{T^2(z^2+2z+1)}{T^2(z^2+2z+1)+4T\tau(z^2-1)+4\tau^2(z^2-2z+1)}=\newline =\frac{T^2z^2+2T^2z+T^2}{z^2(T^2+4T\tau+4\tau^2)+z(2T^2-8\tau^2)+(T^2-4T\tau+4\tau^2)}\mid\frac{/z^2}{/z^2}

Получаем передаточную функцию дискретного аргумента z для низкочастотного фильтра 2-го порядка:

\frac{Y(z)}{X(z)}=\frac{T^2+2T^2z^{-1}+T^2z^{-2}}{(T^2+4T\tau+4\tau^2)+(2T^2-8\tau^2)z^{-1}+(T^2-4T\tau+4\tau^2)z^{-2}}

Следующим шагом будет получение конечно-разностного уравнения, которое уже можно непосредственно использовать в реальном мире (на практике, например, в програме, микроконтроллере и т. д. и т. п.). Перепишем полученное выражение:

T^2X[z]+2T^2z^{-1}X[z]+T^2z^{-2}X[z]=\newline =(T^2+4T\tau+4\tau^2)Y[z]+(2T^2-8\tau^2)z^{-1}Y[z]+(T^2-4T\tau+4\tau^2)z^{-2}Y[z]

От полученного выражения берём обратное Z-преобразование (пользуясь таблицей преобразований, например, из Википедии):

T^2X[i]+2T^2X[i-1]+T^2X[i-2]=\newline =(T^2+4T\tau+4\tau^2)Y[i]+(2T^2-8\tau^2)Y[i-1]+(T^2-4T\tau+4\tau^2)Y[i-2]

Переписываем выражение относительно Y[i]:

Y[i]=\frac{T^2X[i]+2T^2X[i-1]+T^2X[i-2]-(2T^2-8\tau^2)Y[i-1]-(T^2-4T\tau+4\tau^2)Y[i-2]}{(T^2+4T\tau+4\tau^2)}

Полученное конечно-разностное уравнение и есть результат нашего с вами труда. Именно это выражение можно смело записывать в контроллер и осуществлять фильтрацию. Но мы не ищем лёгких путей:). Полученное выражение можно привести к более красивому и простому виду:

Y[i]=\frac{T^2(X[i]+2X[i-1]+X[i-2])-2(T^2-4\tau^2)Y[i-1]-(T-2\tau)^2Y[i-2]}{(T+2\tau)^2}=\newline =\frac{T^2(X[i]+2X[i-1]+X[i-2])-2(T-2\tau)(T+2\tau)Y[i-1]-(T-2\tau)^2Y[i-2]}{(T+2\tau)^2}

Положим, что:

\alpha=\frac{T-2\tau}{T+2\tau}

Тогда:

Y[i]=(\frac{1+\alpha}{2})^2(X[i]+2X[i-1]+X[i-2])-2{\alpha}Y[i-1]-{\alpha}^2Y[i-2]

Это и есть конечный результат! Наступил тот момент, когда можно продемонстрировать преимущества фильтра низких частот 2-го порядка над аналогичным фильтром 1-го порядка. Как и раньше, воспользуемся услугами прекрасного програмного продукта Anaconda, который поддерживает язык програмирования Python. Ниже, представлен пример кода, демонстрирующий работу фильтров:

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import copy
import numpy as np
import matplotlib.pyplot as plot
 
Fs = 8000 # Sampling frequency (Hz)
Fc = 2000 # Cutoff frequency (Hz)
 
T = 1.0/Fs
 
N = Fs # Number of samples per 1 s
 
x = range(0, N)
y = range(0, N)
y1 = range(0, N)
y2 = range(0, N)
 
f1 = 100 # Hz
f2 = 3500 # Hz
 
N_TO_DISPLAY = 2*Fs/f1
 
# Generating test signal
for t in range(0, N):
    x[t] = (np.sin(2.0*np.pi*f1*t/Fs) + np.sin(2.0*np.pi*f2*t/Fs))/2.0
 
# Calculating necessary values for filter
dt = 1.0/Fs
tau = 1.0/(2.0*np.pi*Fc)
alpha = dt / (dt + 2*tau)
 
# Applying filter
for i in range(1, N):
    y[i] = alpha*(x[i]+x[i-1]) + (1-2*alpha)*y[i-1]
 
#copying the filtered result after 1st order LPF
y1 = copy.copy(y)
 
#direct filtering formula for the 2nd order RC LPF without any optimization
#for i in range(2, N):
#    y[i] = (pow(T,2)*x[i] + 2*pow(T,2)*x[i-1] + pow(T,2)*x[i-2] - (2*pow(T,2)-8*pow(tau,2))*y[i-1] - (pow(T,2)-4*T*tau+4*pow(tau,2))*y[i-2]) / (pow(T,2)+4*T*tau+4*pow(tau,2))
 
alpha = (T-2*tau)/(T+2*tau)
 
#filtering formula for the 2nd order RC LPF with some sort of optimization
for i in range(2, N):
    y[i] = pow((1+alpha)/2,2)*(x[i]+2*x[i-1]+x[i-2]) - 2*alpha*y[i-1] - pow(alpha,2)*y[i-2]
 
#copying the filtered result after 2nd order LPF
y2 = copy.copy(y)
 
# Drawing original and filtered signals
plot.clf()
plot.xlabel("Time")
plot.ylabel("Sample value")
plot.plot(x[0:N_TO_DISPLAY], 'r')
plot.plot(y1[0:N_TO_DISPLAY], 'g')
plot.plot(y2[0:N_TO_DISPLAY], 'b')
plot.draw()

А вот и результат работы програмы:

LPF_1st_order_vs_LPF_2nd_orderНа этом рисунке, красным цветом показан входной, зашумленный сигнал (сумма синусоид с частотами 100 Гц (полезный сигнал) и 3500 Гц (помеха)). Зелёным цветом показан результат фильтрации с использованием фильтра 1-го порядка, синим — результат работы фильтра 2-го порядка. Можно сделать выводы о том, что фильтр 2го порядка, при одинаковой частоте среза, фильтрует зашумленный сигнал лучше, чем фильтр 1го порядка. На этом, пожалуй, все. О методике оценки амплитудно-частотной и фазо-частотной характеристики я писал в контексте рассмотрения фильтра первого порядка.

P.S. Вот такие они, простые (на первый взгляд) цифровые фильтры. Естественно, в данном обзоре я не осветил все подробности конструирования цифровых фильтров, но повторюсь, что цель статьи не рассказать обо всем сразу, а дать толчёк для лёгкого старта. Надеюсь, у меня это получилось! Удачи!:)

P.P.S. Обратите внимание, трансформация частотной оси «frequency-prewrapping» была преднамеренно упущена мною, дабы не отпугнуть читающего. Тем более что эта версия фильтра имеет право на существование. Т. е. если сравнивать версию фильтра С использованием «frequency-prewrapping» и БЕЗ него, то АЧХ фильтра будут схожи в диапазоне частот от 0 до частоты дискретизации делённой на 8 (приблизительное значение взято из моей памяти!). Детальнее об «frequency-prewrapping» я собирался написать и напишу в следующих статьях!

P.P.P.S. А вот и долгожданная статья о трансформации частотной оси! В качестве примера, рассмотрен фильтр низких частот первого порядка.