Июл02

Трансформация частотной оси и создание полноценного БИХ фильтра

by admin on 02.07.2015 at 19:46
Опубликовано: Цифровая обработка аудио

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

Что же такое «frequency prewarping» ? Начну издалека. Когда мы использовали Билинейное преобразование, была осуществлена замена:

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

тем самым мы аппроксимировали непрерывную частотную ось в её дискретный аналог:

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

где \omega_a — угловая частота на непрерывной (аналоговой) частотной оси.

Поскольку z = e^{j \omega_d T}, тогда можно записать:

s = j\omega_a = \frac{2}{T} \frac{1 - e^{-j \omega_d T}}{1 + e^{-j \omega_d T}} = \frac{2}{T} \frac {e^{\frac{j\omega_dT}{2}} - e^{\frac{-j\omega_dT}{2}}} {e^{\frac{j\omega_dT}{2}} + e^{\frac{-j\omega_dT}{2}}} = \frac{2}{T} tanh(\frac{j\omega_dT}{2})

где \omega_d — угловая частота на дискретной частотной оси.

В результате получаем ключевое выражение:

\omega_a = \frac{2}{T}tan(\frac{\omega_dT}{2})

Представляю, каково вам сейчас… Ничего не понятно. Формулы действительно пугают, особенно гиперболический тангенс но, увы, без него никак. Для чего всё это? Дело в том, что если не брать во внимание искривление частотной оси при использовании Билинейного преобразования, то может, случится следующее: «Мы настроим частоту среза нашего с вами фильтра на 2.7 кГц при частоте дискретизации 8 кГц и будем ожидать, что пропустив синусоиду, с частотой 2.7 кГц и амплитудой в 1 (например, Вольт) на выходе получим синусоиду, на 3 дБ меньшей амплитуды, чем была на входе, т. е. ожидаем получить синусоиду с амплитудой 0.707 (Вольт).». Если проверить это на разработанном ранее фильтре, то результат огорчит. Дело как раз в трансформации частотной оси. Мы все сделали для частоты 2.7 кГц, но в процессе Билинейного преобразования не было учтено то, что данное значение искажается.

Как использовать полученное выражение при разработке реального фильтра мы увидим на примере старого доброго RC фильтра, который мы разрабатывали ранее. В статье «Как получить цифровой БИХ фильтр из аналогового или как самому «изобрести» цифровой БИХ фильтр» была получена зависимость образа выходного сигнала от входного:

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

Дабы не повторяться, продолжим с этого места. Выше было сказано, что мы должны использовать компенсацию частотной шкалы:

\omega_a = \frac{2}{T}tan(\frac{\omega_dT}{2})

Действительно, не понятно, куда можно подставить данное выражение. Вместо s подставляем \frac{2}{T}\frac{z-1}{z+1}, но никаких угловых частот \omega_a тут нет :-( . Зато мы имеем в наличии непонятное \tau (постоянная времени RC цепи) :-) . Предлагаю получить зависимость угловой частоты от постоянной времени фильтра (такова зависимость есть, уж поверьте мне). Заодно узнаем, откуда она берётся. Ведь в статье Википедии, на которую я ссылаюсь, этого не указано. Вы получите бесценные знания. Берётся данное значение с выражения, которым описывают АЧХ системы. Нам безумно повезло, ведь это выражение именно для этого типа фильтра мы получили ранее в статье «Амплитудно-частотная и фазо-частотная характеристики RC фильтра низких частот (метод преобразования Лапласа)»:

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

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

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

Отсюда получаем (а не из пальца, кстати):

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

Напомню, мы вывели зависимость постоянной времени от угловой частоты среза \omega_{cutoff} по уровню в -3 дБ. Подставляем полученную зависимость:

\frac{Y(s)}{X(s)} = \frac{1}{1+s\tau} = \frac{1}{1+\frac{s}{\omega_{cutoff}}} = \frac{1}{1+\frac{s}{\omega_a}}

Ну, это уже совсем другое дело. Теперь не составит труда осуществить подстановку для s и \omega_a:

\frac{Y(s)}{X(s)} = \frac{1}{1+\frac{s}{\omega_a}} = \frac{1}{1+\frac{\frac{2}{T}\frac{z-1}{z+1}}{\frac{2}{T}tan(\frac{\omega_dT}{2})}} = \frac{1}{1+\frac{\frac{z-1}{z+1}}{tan(\frac{\omega_dT}{2})}}

Давайте, для удобства, положим:

K = \frac{1}{tan(\frac{\omega_dT}{2})}

Тогда:

\frac{Y(z)}{X(z)} = \frac{1}{1+\frac{z-1}{z+1}K} = \frac{z+1}{z+1+(z-1)K} = \frac{1+z^{-1}}{1+z^{-1}+K-Kz^{-1}} = \frac{1+z^{-1}}{1+K+z^{-1}(1-K)}

Преобразовываем выражение:

X(z)+X(z)z^{-1} = Y(z)+KY(z)+(1-K)Y(z)z^{-1}

Пришло время вернуться в реальный мир, Нео! Делаем обратное z преобразование:

x[i]+x[i-1] = y[i]+Ky[i]+(1-K)y[i-1] = (1+K)y[i]+(1-K)y[i-1]

Получаем результирующее разностное выражение для выходного отчёта:

y[i] = \frac{x[i]+x[i-1]-(1-K)y[i-1]}{1+K}

Давайте чуть-чуть перепишем выражение. Если положить, что:

\alpha = \frac{1}{1+K}

Тогда:

y[i] = \frac{x[i]+x[i-1]-(1-K)y[i-1]}{1+K} = \alpha(x[i]+x[i-1])+(1-2\alpha)y[i-1]

Кроме этого, выражение для \alpha можно расписать и упростить:

\alpha = \frac{1}{1+K} = \frac{1}{1+\frac{1}{tan(\frac{\omega_dT}{2})}} = \frac{tan(\frac{\omega_dT}{2})}{tan(\frac{\omega_dT}{2})+1} = \frac{tan(\frac{2{\pi}F_c\frac{1}{F_s}}{2})}{tan(\frac{2{\pi}F_c\frac{1}{F_s}}{2})+1} = \frac{tan(\pi\frac{F_c}{F_s})}{tan(\pi\frac{F_c}{F_s})+1} = 1 - \frac{1}{ tan(\pi\frac{F_c}{F_s})+1}

где F_c и F_s — частота среза по уровню -3 дБ и частота дискретизации (семплирования) соответственно.

Ну, всё! Я вас поздравляю! Теперь вы знаете, как создаются (или же изобретаются) полноценные фильтры. Не оставим этот пост без примера использования. Вот он на удобном и всеми любимом 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

# -*- coding: utf-8 -*-
"""
Created on Tue Jun 30 21:52:15 2015
 
@author:  Viktor Signaievskyi
@blog:    www.HowToDoIt.com.ua
@e-mail:  piligrim2007@meta.ua
 
@details: Example program of digital realization of 1st order RC filter with 
          using of frequency prewarping technique
          
Copyrighted by ©Viktor Signaievskyi
"""
 
import numpy as np
import matplotlib.pyplot as plot
 
Fs = 8000 # Sampling frequency (Hz)
Fc = 2001 # Cutoff (-3dB) frequency (Hz)
 
N = Fs # Number of samples per 1 s
 
x = range(0, N)
y = [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 alpha coefficient for the filter
alpha = 1.0 - 1.0 / (np.tan(np.pi*Fc/Fs)+1.0)
 
# Applying filter
for i in range(1, N):
    y[i] = alpha*(x[i]+x[i-1]) + (1-2*alpha)*y[i-1]
 
# Drawing original and filtered signals
plot.clf()
plot.xlabel("Time")
plot.ylabel("Sample value")
plot.plot(x[0:N_TO_DISPLAY], 'r')
plot.plot(y[0:N_TO_DISPLAY], 'g')
plot.draw()

 Тем, кому сложно устанавливать Python и запускать данный скрипт, могу посоветовать протестировать его с использованием Haskell. Ну, а если и это сложно :-) , то просто посмотрите на результирующий график:

FilterInActionКак видим, фильтр работает! Более интересно посмотреть в действии работу фильтра в режиме подавления по уровню -3 дБ на конкретной частоте. Для этого тестовым сигналом должна выступать, например, синусоида и её частота должна совпадать с частотой среза F_c (cutoff frequency). Пример кода:

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

# -*- coding: utf-8 -*-
"""
Created on Tue Jun 30 21:52:15 2015
 
@author:  Viktor Signaievskyi
@blog:    www.HowToDoIt.com.ua
@e-mail:  piligrim2007@meta.ua
 
@details: Example program of digital realization of 1st order RC filter with 
          using of frequency prewarping technique
          
Copyrighted by ©Viktor Signaievskyi
"""
 
import numpy as np
import matplotlib.pyplot as plot
 
Fs = 8000 # Sampling frequency (Hz)
Fc = 1400 # Cutoff (-3dB) frequency (Hz)
 
N = Fs # Number of samples per 1 s
 
x = range(0, N)
y = [0]*N
 
N_TO_DISPLAY = 10*Fs/Fc # Number of periods of Fc to display
 
# Generating test signal
for t in range(0, N):
    x[t] = np.sin(2.0*np.pi*Fc*t/Fs)
 
# Calculating alpha coefficient for the filter
alpha = 1.0 - 1.0 / (np.tan(np.pi*Fc/Fs)+1.0)
 
# Applying filter
for i in range(1, N):
    y[i] = alpha*(x[i]+x[i-1]) + (1-2*alpha)*y[i-1]
 
# Drawing original and filtered signals
plot.clf()
plot.xlabel("Time")
plot.ylabel("Sample value")
plot.plot(x[0:N_TO_DISPLAY], 'r')
plot.plot(y[0:N_TO_DISPLAY], 'g')
plot.draw()

Вот что у нас получается на графике:

FilterInActionAttenuationКак видим, мы действительно получаем свои -3 дБ (0.707). И да — это синусоида :-) .

P.S. Внимательный читатель заметит, что полученное разностное выражение точностью совпадает с полученным ранее и это не совпадение. Я постарался привести результат именно к такому виду. Ну, а разница заключается всего лишь в расчете коэффициента  \alpha.

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