Трансформация частотной оси и создание полноценного БИХ фильтра
on 02.07.2015 at 19:46В этой статье мы наконец-то разработаем полноценный фильтр низких частот. Я расскажу о том, что такое трансформация частотной оси (frequency prewarping), откуда она берётся и зачем нам её использовать. В предыдущей статье «Как получить цифровой БИХ фильтр из аналогового или как самому «изобрести» цифровой БИХ фильтр» этого блога мы с вами рассмотрели основы основ в плане создания реально работающего цифрового дискретного фильтра. Прошли путь от физической модели RC цепочки до получения разностного выражения, которое можно использовать в реальной жизни, запрограммировав выражение в память какого-либо устройства. В конце упомянутой статьи есть послесловие:
«P.P.S. Обратите внимание, трансформация частотной шкалы «frequency prewarping» была преднамеренно упущена мною, дабы не отпугнуть читающего, тем более что эта версия фильтра имеет право на существование и без этой трансформации. Т. е. если сравнивать версию фильтра С использованием «frequency prewarping» и БЕЗ него, то АЧХ фильтра будут схожи в диапазоне частот от 0 до частоты дискретизации делённой на 8 (приблизительное значение взято из моей памяти!). Детальнее об «frequency prewarping» я собирался написать и напишу в следующих статьях!»
Включать в и так длинную статью, наполненную массой сложного и непонятного материала для начинающего, было бы нецелесообразно. Тем, кто не согласен, предлагаю разозлиться на автора и продолжить чтение данной статьи .
Что же такое «frequency prewarping» ? Начну издалека. Когда мы использовали Билинейное преобразование, была осуществлена замена:
тем самым мы аппроксимировали непрерывную частотную ось в её дискретный аналог:
где — угловая частота на непрерывной (аналоговой) частотной оси.
Поскольку , тогда можно записать:
где — угловая частота на дискретной частотной оси.
В результате получаем ключевое выражение:
Представляю, каково вам сейчас… Ничего не понятно. Формулы действительно пугают, особенно гиперболический тангенс но, увы, без него никак. Для чего всё это? Дело в том, что если не брать во внимание искривление частотной оси при использовании Билинейного преобразования, то может, случится следующее: «Мы настроим частоту среза нашего с вами фильтра на 2.7 кГц при частоте дискретизации 8 кГц и будем ожидать, что пропустив синусоиду, с частотой 2.7 кГц и амплитудой в 1 (например, Вольт) на выходе получим синусоиду, на 3 дБ меньшей амплитуды, чем была на входе, т. е. ожидаем получить синусоиду с амплитудой 0.707 (Вольт).». Если проверить это на разработанном ранее фильтре, то результат огорчит. Дело как раз в трансформации частотной оси. Мы все сделали для частоты 2.7 кГц, но в процессе Билинейного преобразования не было учтено то, что данное значение искажается.
Как использовать полученное выражение при разработке реального фильтра мы увидим на примере старого доброго RC фильтра, который мы разрабатывали ранее. В статье «Как получить цифровой БИХ фильтр из аналогового или как самому «изобрести» цифровой БИХ фильтр» была получена зависимость образа выходного сигнала от входного:
Дабы не повторяться, продолжим с этого места. Выше было сказано, что мы должны использовать компенсацию частотной шкалы:
Действительно, не понятно, куда можно подставить данное выражение. Вместо подставляем
, но никаких угловых частот
тут нет
. Зато мы имеем в наличии непонятное
(постоянная времени RC цепи)
. Предлагаю получить зависимость угловой частоты от постоянной времени фильтра (такова зависимость есть, уж поверьте мне). Заодно узнаем, откуда она берётся. Ведь в статье Википедии, на которую я ссылаюсь, этого не указано. Вы получите бесценные знания. Берётся данное значение с выражения, которым описывают АЧХ системы. Нам безумно повезло, ведь это выражение именно для этого типа фильтра мы получили ранее в статье «Амплитудно-частотная и фазо-частотная характеристики RC фильтра низких частот (метод преобразования Лапласа)»:
Напомню, что для нас значит данное выражение. При постоянном значении мы изменяем частоту
и смотрим на коэффициент передачи системы (получившееся число). На каких-то частотах мы будем получать значения, близкие к единице (фильтр пропускает данную частоту), на каких-то — близкие к нулю (фильтр подавляет данную частоту), а на всех остальных — значения между 0 и 1
(фильтр частично подавляет частоту). Обычно, фильтры разрабатывают для частоты среза
, которая составляет -3 дБ (
от амплитуды входной гармонической составляющей, например, синусоиды). Что такое децибел — читаем тут. В таком случае можно записать:
Отсюда получаем (а не из пальца, кстати):
Напомню, мы вывели зависимость постоянной времени от угловой частоты среза по уровню в -3 дБ. Подставляем полученную зависимость:
Ну, это уже совсем другое дело. Теперь не составит труда осуществить подстановку для и
:
Давайте, для удобства, положим:
Тогда:
Преобразовываем выражение:
Пришло время вернуться в реальный мир, Нео! Делаем обратное z преобразование:
Получаем результирующее разностное выражение для выходного отчёта:
Давайте чуть-чуть перепишем выражение. Если положить, что:
Тогда:
Кроме этого, выражение для можно расписать и упростить:
где и
— частота среза по уровню -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. Ну, а если и это сложно , то просто посмотрите на результирующий график:
Как видим, фильтр работает! Более интересно посмотреть в действии работу фильтра в режиме подавления по уровню -3 дБ на конкретной частоте. Для этого тестовым сигналом должна выступать, например, синусоида и её частота должна совпадать с частотой среза
(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() |
Вот что у нас получается на графике:
Как видим, мы действительно получаем свои -3 дБ (0.707). И да — это синусоида
.
P.S. Внимательный читатель заметит, что полученное разностное выражение точностью совпадает с полученным ранее и это не совпадение. Я постарался привести результат именно к такому виду. Ну, а разница заключается всего лишь в расчете коэффициента .
Удачи в познании мира ЦОС!
Ну, в принципе, статья вполне годная.
А не проясните такой момент — чем же является Ваш «полноценный БИХ- фильтр» при значении параметра К=1?
Спасибо, Santik!
, а значит
. В этом случае мы получим буферное устройсво, которое не будет искажать входной сигнал. Если перенести на реальную RC цепочку, то этот случай можно представить, например, как
.
для RC фильтров разного порядка.
При
P.S. Буквально на днях собирался выложить статью, описывающую логику работы с
Но это не так! При
!
— это не
!!! Если
, то и
и тогда действительно получаем всепропускающий фильтр с
.
.
Т.е. я хочу сказать, что случай
Но меня, я повторю вопрос, интересует случай
Но всё же, это так
, если взять во внимание то, что символом К обозначена АЧХ системы. Если же взять во внимание тот факт, что тем же символом обозначено временное выражение, то Ваш вопрос звучит совсем по-другому, и я его понимаю.
.
становится единицей в том случае, если частота среза будет равняться четверти от частоты дискретизации. В этом частном случае разностное уравнение приобретает вид скользящего среднего (в данном случае — разновидность КИХ фильтра).
.
Честно говоря меня удивляет, что Вы сами не смогли разобраться с этим вопросом
Santik, если Вам ещё что-то непонятно, пожалуйста — обращайтесь! Как говорится: «Век живи — век учись!»