Июл19

Как получить цифровой БИХ фильтр из аналогового или как самому «изобрести» цифровой БИХ фильтр

by admin on 19.07.2014 at 17:27
Опубликовано: DSР, Цифровая обработка аудио

Несмотря на свою столь широкую распространенность, разработка цифровых БИХ фильтров (фильтры с Бесконечной Импульсной Характеристикой), даже в наше время, время высоких технологий, для многих остается загадкой. Время от времени мне попадаются статьи или же видео о том, каков принцип разработки БИХ фильтров. Становится печально, когда просмотрев очередное из них я осознаю, что благодаря именно такой подаче материала большинству учащихся действительно не остается ничего другого, как опустить руки и придти в унылое состояние. У многих моих знакомых имеются практические навыки в разработке БИХ фильтров, но они заключаются в банальном использовании, например пакета Matlab, для получения коэффициентов того или иного БИХ фильтра, иногда даже без понимания сути его работы. Лично я, пока полностью не разберусь с принципом работы того или иного участка кода, не могу с уверенностью заявлять о его надёжности. Возобновив свои университетские знания я решил поделиться с Вами своими познаниями в разработке цифрового БИХ фильтра «с нуля»! Дальнейшее изложение материала будет происходить в контексте разработки реального цифрового ФНЧ (фильтра низких частот) на основе аналогового прототипа RC — фильтра. Мы не будем вдаваться во все подробности преобразования Лапласа и Z-преобразования. Я поставил себе за задачу показать читателю сам скелет разработки и мне кажется для старта этого будет вполне достаточно. Все недостающие знания при желании можно без труда почерпнуть из Интернета. «Что нового в этой статье?» — спросите вы. Новшество заключается в том, что в данном относительно коротком посте будут собраны воедино такие вещи как:

  • история возникновения аналоговой, а позже и цифровой фильтрации;
  • описание схемы с помощью математической модели;
  • преобразование полученной модели в другую, более удобную, с помощью преобразования Лапласа;
  • будет объяснено, зачем нужен предыдущий пункт;
  • Z преобразование;
  • получение разностного уравнения (результат работы).

Итак, начнём! Откуда же взялись БИХ фильтры? Как часто бывает, они взяты с уже годами проверенной практикой схемотехнических решений. Исторически так сложилось, что такой электронный компонент как конденсатор был разработан задолго (1745 год) до разработки, скажем радиоприёмника (1895 год). Так вот, ещё на заре развития аудиотехники люди разрабатывали радиоприёмники. В любом из таких устройств в наличии конечно же имелись наушники или динамик. Разработчики заметили, что если подключать конденсатор к динамику (параллельно или же последовательно), то в звуке что-то меняется. Однозначно можно было подобрать конденсатор таким образом, что звук становился чище, т.е. конденсатор уменьшал количество шума (при параллельном подключении к наушникам). Время шло, наука стремительно развивалась и вскоре выросла целая физико-математическая машина, с помощью которой стало возможным объяснить подобные явления. Как мы знаем, цифровая фильтрация до недавнего времени не сильно использовалась из-за невозможности её физической реализации, ведь процессоры обладали низкими значениями тактовой частоты. В виду этого развивалась индустрия фильтрации в аналоговом виде. К слову, развивалась она довольно хорошо и что самое главное качественно. Прошло время и процессоры, с частотами выше 100 МГц, стали реальностью. Стало реальным производить фильтрацию используя механизм БИХ фильтров. Но для того, чтобы использовать модель проверенных годами БИХ фильтров в цифровой фильтрации их необходимо каким-то образом описать и преобразовать математической моделью с целью использования в цифровом виде. Дальше, мы с Вами пройдем путь от аналоговой RC цепочки до конечной реализации цифрового БИХ фильтра конечно же с примером программы в конце статьи, так как статья ориентированна на практиков, а не только теоретиков. Итак, имеем следующую схему:1.Scheme

Видим, в составе резистор и конденсатор. На вход подается напряжение, которое меняется во времени, на выходе будем иметь отличное от входного напряжение, которое, в свою очередь, тоже меняется во времени. Ток, который течет в цепи можно описать, используя следующее выражение для конденсатора:

2.Current

В таком случае можно записать:

3.Equation

Только ради удобства сделаем замену:

4.SubstitutionПолучаем следующее выражение: 3.DifferentialEquation На данном шаге нами было получено выражение, которые в полной мере описывает поведение нашей RC цепи. Но вместо простой зависимости выхода от входа, мы получили что-то, на первый взгляд страшное, ведь мы имеем производную выходного сигнала по времени. В общем главной нашей задачей является ответ на следующий вопрос: «Что нужно сделать с входным сигналом, чтобы получить выходной? Умножить, суммировать, сдвинуть во времени, что?«. Каким бы наивным не казался этот вопрос, но он является ключевым. Все последующие действия, связанные с математическими преобразованиями, будут нацелены именно на то, чтобы можно было получить чёткую запись, которая покажет зависимость выхода от входа. Иными словами, нам необходимо решить полученное дифференциальное уравнение. Конечно же, решать его есть возможным множественным количеством методов, но одним из мне известных, универсальным, является метод с использованием преобразования Лапласа. Называется этот метод в математике методом операционного исчисления. Что это такое? Я специально не буду писать формул, а постараюсь объяснить «на пальцах». Допустим мы имеем функцию (в дискретном виде ее можно представить, скажем, последовательностью чисел, осциллограммой). Так вот, под преобразованием Лапласа понимают разложение начальной функции в сумму затухающих, незатухающих и возрастающих осциллирующих экспонент. Вот как они выглядят: 5.ExponentialOscillations Некоторые из читающих могут провести некую параллель с преобразованием Фурье. Действительно, преобразование Фурье является частным случаем для преобразования Лапласа при константной (незатухающей) осциллограмме.

6.Sinewave

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

7.LaplaceTransform

Ну, уже совсем другое дело. Теперь не составит труда показать зависимость образа выходного сигнала от входного:

8.FrequencyResponse

где:

9.TimeConstant

Полученное отношение образов можно рассматривать как отношение выходного спектра (экспоненциальных осциллограмм) к входному спектру (экспоненциальных осциллограмм). Что же такое переменная S? Тут посложнее, хотя это всего лишь комплексная переменная. Откуда она тут взялась? Из преобразования Лапласа. Тут можно оценить достижения математики, так как благодаря операциям с комплексными числами стало возможным получать решение в более элегантном представлении. В общем я обещал не вдаваться в детали. Только скелет (ну, с минимальными отступлениями)! Должен признать, дальше будет ещё сложнее. Но если не вдаваться в подробности, то всё довольно таки просто. Следующим шагом будет переход от преобразования Лапласа к Z-преобразованию (оно же преобразование Фишера или же преобразование Лорана). Z-преобразование это то же самое преобразование Лапласа, но вот только для дискретного (а не непрерывного) сигнала. К счастью для нас, существует простой способ перехода от непрерывного преобразования Лапласа к дискретному Z-преобразованию и называется оно «Билинейное преобразование». Пользоваться им очень просто. Достаточно в преобразовании Лапласа произвести замену:

10.BilinearTransformFormula

После подстановки, получим следующее:

11.Z-Transform

Напомню, пока мы работаем в частотной области, не во временной. Наша задача получить зависимость выхода от входа именно во временной области, а не в частотной. Нужно что-то сделать для перехода от частотной во временную область. Этим «что-то» является обратное Z-преобразование. Оно также имеет свою таблицу преобразований и нашей задачей является приведение уравнения до вида «табличных составляющих».  Если переписать полученное выражение в виде:

12.DifferenceEquation1то получить обратное Z-преобразование очень просто:

12.DifferenceEquation2

Чтобы читатель понимал, мы уже находимся во временной области. Перепишем полученное выражение:

12.DifferenceEquation3

Полученное выражение можно считать конечным результатом, но если произвести замену:

13.DifferenceEquationSubstitutionто получим более элегантную запись:

14.DifferenceEquationFinalПоздравляю нас с Вами! Мы получили выражение, пользуясь которым возможно фильтровать дискретный сигнал. Для тех кто не понимает как применить данное выражение на практике привожу пример кода, написанного на языке программирования 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

import numpy as np
import matplotlib.pyplot as plot
 
Fs = 8000 # Sampling frequency (Hz)
Fc = 2000 # Cutoff frequency (Hz)
 
N = Fs # Number of samples per 1 s
 
x = range(0, N)
y = 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]
 
# Drawing original and filtered spectrums
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()

 

В этом коде начальным сигналом является сумма двух синусоид с частотой 100Гц и 3500Гц. Предположим, что синусоида с частотой 3500 Гц — это помеха и мы хотим, с помощью созданного нами цифрового фильтра, отфильтровать её. Частота сэмплирования в данном случае равна 8кГц, а частота среза была выбрана равной 2кГц. В результате отработки данного скрипта получим следующий рисунок:

15.FilterInAction

На данном графике красным цветом показан начальный сигнал (полезный сигнал+помеха), зелёным цветом показан отфильтрованный сигнал (полезный сигнал). Как видим, фильтрация удалась!

И в завершение… Надеюсь, у меня получилось изложить тему доступным Вам, читателю, языком? В следующей статье мы рассмотрим амплитудно-частотную и фазо-частотную характеристики разработанного нами с Вами низкочастотного RC фильтра. Ждите новые посты касательно цифровой обработки сигналов в контексте их применения на практике. Я же буду ждать Ваши отзывы внизу страницы! Вы даже не представляете как приятно получать такого рода отдачу:) Всего хорошего!)

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

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

P.P.P.S. И ещё. Дабы не возникало трудностей на этапе получения конечной формулы при подстановке \alpha=\frac{T}{T+2\tau}, будет правильным показать данный переход:

-(\frac{T}{T+2\tau} - \frac{2\tau}{T+2\tau}) = \frac{2\tau}{T+2\tau} - \frac{T}{T+2\tau} = \frac{2\tau-T}{T+2\tau} = \frac{(2\tau+T)-2T}{T+2\tau} = 1-2\alpha