close

Вход

Забыли?

вход по аккаунту

?

Gilmutdinov

код для вставкиСкачать
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное автономное
образовательное учреждение высшего образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
ВВЕДЕНИЕ В ЦИФРОВУЮ ОБРАБОТКУ
СИГНАЛОВ
Учебное пособие
Санкт-Петербург
2016
УДК 004.627
ББК 30
Г47
Рецензенты:
кандидат технических наук, доцент А. Б. Сергиенко;
кандидат технических наук, доцент Е. А. Бакин
Утверждено
редакционно-издательским советом университета в
качестве учебного пособия
Гильмутдинов, М. Р.
Г47
Введение в цифровую обработку сигналов: учеб. пособие / М. Р. Гильмутдинов, Н. Д. Егоров, А. И. Веселов. – СПб.: ГУАП, 2016. – 95 с.
ISBN 978-5-8088-1153-9
Рассматриваются теоретические и практические аспекты современных
методов формирования и преобразования цифровых сигналов. Рассмотрение
теоретических вопросов сопровождается описанием практических работ для
самостоятельного исследования. Предназначено для студентов, обучающихся по
направлениям «Информационная безопасность» и «Инфокоммуникационные
технологии и системы связи», а также для самостоятельной работы аспирантов.
УДК 004.627
ББК 30
ISBN 978-5-8088-1153-9
© Гильмутдинов М. Р., Егоров Н. Д.,
Веселов А. И., 2016
© Санкт-Петербургский государственный
университет аэрокосмического
приборостроения, 2016
ПРЕДИСЛОВИЕ
Цифровая обработка сигналов является фундаментальной дисциплиной,
на основе которой базируется большинство современных методов представления, анализа, передачи и хранения информации. В связи с этим понимание
основных принципов цифровой обработки сигналов является важным требованием для будущих специалистов в данных областях.
Материал, изложенный в настоящем пособии, является вводным для последующего углубленного изучения приложений цифровой обработки сигналов. Целями данного пособия являются:
1. изучение способов формирования цифровых сигналов с помощью методов дискретизации и квантования,
2. исследование методов анализа сигналов с помощью спектральных преобразований,
3. анализ базовых подходов в обработке сигналов с помощью линейных систем.
Основными источниками при составлении данного учебного пособия являлись учебное пособие А.Б. Сергиенко «Цифровая обработка сигналов»,
книги Ю. Сато «Без паники! Цифровая обработка сигналов», И.Е. Бочаровой «Compression for Multimedia» и материалы сайта dsplib.ru, созданного и
поддерживаемого С. Бахуриным.
3
I. ДИСКРЕТИЗАЦИЯ АНАЛОГОВОГО
СИГНАЛА. ФУРЬЕ-АНАЛИЗ
СИГНАЛОВ
Цели исследования: изучение свойств дискретизации аналоговых сигналов; изучение методов Фурье-анализа аналоговых сигналов.
1. Теоретические основы
1.1. Основные понятия
Сигналом (от латинского signum – знак) обычно называют намеренное
изменение состояния (возмущение) какой-либо физической величины, целью которого является передача (или сообщение) информации. Сигналы, как
правило, используют для передачи информации в какой-либо физической
среде (металлический проводник, воздух) или в вакууме. Наиболее распространенными средствами возмущения является импульс или волна, например, электромагнитная или акустическая. Таким образом, для передачи информации в пространстве и/или во времени, требуется формирование импульса с заданными параметрами или волны, параметры которой изменяются по определенному закону.
Для формирования (регистрации) информации указанным способом
требуется специальное устройство, называемое источником, а для получения (воспроизведения) — устройство, называемое приемником. Работа обоих устройств должна быть согласованной, т.е. приемник должен знать заранее или определять в процессе передачи, какой вид информации, по какому
правилу передается и как эту информацию интерпретировать. Полная схема
от формирования сигнала до его воспроизведения изображена на рис. 1.1,
однако в данном учебном пособии подробно будут рассматриваться методы,
относящиеся к процессам формирования цифрового сигнала из аналогового и обратно, а также методы преобразования (фильтрации) аналогового и
цифрового сигналов.
4
u(t)
Регистратор
Сигнал’
Проигрыват u'(t)
ель
АналогоЦифровой
Преобразователь
un
Кодер
Источника
Кодер
Канала
Модулятор
ЦифроАналоговый
Преобразователь
u'n
Декодер
Источника
Декодер
Канала
Демодулятор
Канал
Сигнал
Рис. 1.1. Типовая схема формирования и передачи сигналов
Одним из основополагающих понятий в теории сигналов является понятие спектра. Любой сигнал можно представить в виде суммы некоторых составляющих с заданными (известными заранее) параметрами. Совокупность
составляющих заданной формы, образующих в сумме требуемый сигнал, называют спектром, а сами составляющие сигнала –– спектральными составляющими [1]. В зависимости от типа сигнала, его спектр может состоять из
конечного или бесконечного числа составляющих. Наиболее популярными в
технике являются спектры на базе гармонических составляющих, например,
тригонометрических рядов Фурье.
В данном исследовании будут рассматриваться сигналы, являющиеся
непрерывными или дискретными функциями от времени. В следующих подразделах будет приведена классификация сигналов по различным параметрам, представлен математический аппарат, используемый для описания сигналов и спектров, а также приведено описание ряда Фурье и интеграла Фурье, применяемых для спектрального представления непрерывных периодических и непериодических сигналов.
1.2. Классификация сигналов
1. Аналоговые, дискретные и цифровые.
Аналоговый сигнал u(t) является непрерывной функцией от непрерывного аргумента. Все регистрируемые сигналы, существующие в природе,
относятся к этому классу. Наиболее распространенным аргументом является время t, однако в задачах обработки визуальной информации, объекты описываются непрерывной пространственной функцией интенсивности регистрируемого светового потока. Дискретные сигналы u[n], как
5
функции от дискретного аргумента, определены только в фиксированные
моменты времени. Частным случаем дискретного сигнала, используемого в цифровой обработке является цифровой сигнал un , который является дискретным не только по аргументу, но и по возможным значениям.
Примеры сигналов данных классов изображены на рисунке 1.2. В данном исследовании будет рассмотрена задача дискретизации аналогового
сигнала по времени. Наиболее распространенный на практике способ перехода от аналогового сигнала к дискретному и обратно был сформулирован в работах Котельникова [5] и Шеннона [6]. Способам квантования,
используемым в задаче перехода от дискретного представления сигнала к
цифровому, посвящено отдельное исследование.
а)
Рис. 1.2.
б)
в)
Примеры сигналов различных видов: а–– аналоговый; б–– дискретный; в–– цифро-
вой
2. Детерминированные и случайные сигналы.
Детерминированными [2] называют сигналы, вид и параметры которых
известны заранее. Это позволяет рассчитать значение сигнала при любом
значении аргумента.
6
Случайный сигнал описывается случайной функцией того же аргумента.
Конкретное значение этой случайной функции при заданном значении аргумента принимается с заданной вероятностью.
3. Периодические сигналы.
Сигнал u(t) является периодическим, если выполняется следующее равенство:
u(t + nT ) = u(t), при ∀t ∈ R, ∀n ∈ Z.
(1.1)
Минимальную величину T , при которой выполняется данное равенство,
называют периодом сигнала.
4. Гармонические сигналы. Гармонические сигналы задаются с помощью
функции вида:
u(t) = A cos(ωt + ϕ),
(1.2)
где A является амплитудой, ω –– угловой частотой, а ϕ –– начальной фазой
сигнала. В радиотехнических системах гармонический сигнал используют в качестве тестовых данных. Гармонические функции также используются в качестве базисных, например, для рядов Фурье.
Среди множества математических функций следует выделить дельтафункцию Дирака, которая широко используется в теоретических исследованиях сигналов:
{
∞, t = 0
δ(t) =
(1.3)
0, t ̸= 0.
Данная функция не является физически реализуемой. Ее интеграл равен единице:
∫∞
δ(t)dt = 1.
(1.4)
−∞
Важным свойством дельта-функции является фильтрующее свойство [3]:
∫∞
u(t)δ(t − t0 )dt = u(t0 ).
(1.5)
−∞
1.3. Энергия и мощность сигнала
Наиболее распространенными на практике количественными характеристиками, используемыми для измерения и сравнения сигналов, являются
7
энергия и мощность. Данные величины имеют интерпретацию, основанную
на использовании этих понятий в физике. Рассмотрим электрическую схему,
изображенную на рис. 1.3.
U
R
u(t)
Рис. 1.3. Пример электрической схемы
При постоянном напряжении u(t) = U источника, мощность, выделяемая на сопротивлении R, рассчитывается по формуле:
U2
.
(1.6)
R
Количество энергии, выделяемой на сопротивлении в виде тепла, за время T
составляет:
U2
(1.7)
E=
T.
R
В случае переменного напряжения, вычисляется мгновенная мощность,
выделяемая на сопротивлении в момент времени t:
P =
u2 (t)
.
(1.8)
R
Энергия, выделяемая за время T , в этом случае вычисляется с помощью интеграла:
∫T
∫T
1
E = p(t)dt =
u2 (t)dt,
(1.9)
R
p(t) =
0
0
а средняя мощность за то же время составляет:
P =
E
1
=
T
RT
∫T
u2 (t)dt.
0
8
(1.10)
Поскольку при анализе сигналов значение сопротивления R не важно, для
упрощения формулы (1.10) оно принимается за единицу.
Энергию, как правило, используют в качестве характеристики для сигналов конечной длительности. Энергия периодических сигналов бесконечна,
поэтому для периодических сигналов используют усреднение энергии на одном периоде сигнала:
T /2
∫
1
P =
u2 (t)dt.
(1.11)
T
−T /2
1.4. Способы математического описания сигналов
Будем полагать, что значения сигнала являются элементами поля действительных R или комплексных C чисел. Если нет уточнений, то поумолчанию речь идет о поле действительных чисел. Также, если нет уточнений, то речь будет идти о сигнале, который представляет собой непрерывную функцию от времени. Любой сигнал, как непрерывную во времени
функцию u(t) на произвольном интервале [a; b], можно представить набором из N его значений, дискретных по времени. Для описания этого набора
используется N -мерный вектор ⃗u = (u1 , u2 , . . . , uN ). Всевозможные значения N -мерных векторов образуют N -мерное векторное пространство [2],
для которого определены операции сложения (коммутативная группа с существованием нулевого вектора) и умножения вектора на скаляр (с существованием единичного скаляра). Векторное пространство также называют
линейным. Пространство, включающее в себя всевозможные векторы длины
N < ∞ над полем действительных чисел, называют эвклидовым и обозначают RN , если определить для него операцию скалярного произведения векторов [2] (⃗u, ⃗v ) между парой любых векторов ⃗u и ⃗v из этого пространства:
(⃗u, ⃗v ) =
N
∑
ui vi ; ui , vi ∈ R.
(1.12)
i=1
Если пространство определено над полем комплексных чисел, то такое пространство CN будет являться частным случаем гильбертова пространства.
Элементы второго вектора в скалярном произведении (1.12), определенном для пространства над полем комплексных чисел, берутся комплексно9
сопряженными:
(⃗u, ⃗v ) =
N
∑
ui vi∗ ; ui , vi ∈ C.
(1.13)
i=1
Символом (∗ ) обозначается комплексное сопряжение.
При устремлении N к бесконечности векторное пространство становится пространством функций [2]. Вектор в этом случае совпадает с исходным
сигналом:
⃗u −−−−→ u(t)
N →∞
Скалярное произведение двух функций u(t) и v(t) в общем случае определяется формулой:
∫b
(u(t), v(t)) =
γ(t)u(t)v ∗ (t)dt,
(1.14)
a
где γ(t) –– весовая функция. Для простоты ее выбирают в виде константы
1/(b − a), что позволяет вынести γ(t) за знак интеграла. Такой нормирующий множитель полезен тем, что позволяет избежать зависимости от пределов интегрирования, поскольку при увеличении интервала [a; b] значение
скалярного произведения двух функций может неограниченно возрастать.
1
(u(t), v(t)) =
b−a
∫b
u(t)v ∗ (t)dt.
(1.15)
a
Скалярное произведение обладает следующими свойствами, которые
также являются определяющими для эвклидовых пространств (вещественных чисел). Здесь приводятся без доказательства:
– коммутативность: (⃗u, ⃗v ) = (⃗v , ⃗u);
– дистрибутивность: (⃗u + ⃗z, ⃗v ) = (⃗u, ⃗v ) + (⃗z, ⃗v );
– сочетательное свойство: α(⃗u, ⃗v ) = (α⃗u, ⃗v ) = (⃗u, α⃗v ), где α – скаляр, принадлежащий тому же полю;
– неотрицательность: (⃗u, ⃗u) ⩾ 0.
Для описания свойств векторов и функций, а также их отношений друг к
другу в исследовании будут использоваться следующие функционалы и операции.
10
1. Норма. Нормой L2 вектора ∥⃗u∥2 является функционал, заданный в векторном пространстве размерности N и определяющий длину ⃗u в этом пространстве:
√
√
∥⃗u∥2 = (⃗u, ⃗u) = u21 + u22 + · · · + u2N .
(1.16)
Векторное пространство с заданной на нем нормой называют нормированным пространством. В дальнейшем будем рассматривать только L2
норму, поэтому нижний индекс в обозначении использоваться не будет:
∥⃗u∥ = ∥⃗u∥2 .
Если ∥⃗u∥ = 1, то такой вектор называют единичным вектором или ортом.
Любой вектор можно сделать единичным, поделив его на его же норму:
⃗u
.
∥⃗u∥
В пространстве функций для вычисления нормы функции на интервале
[a; b] сумма заменяется интегралом:
v
u
∫b
u
u 1
∥u(t)∥ = t
u2 (t)dt.
(1.17)
b−a
a
2. Расстояние. Расстояние d(⃗u, ⃗v ) между парой векторов ⃗u и ⃗v в N -мерном
пространстве определяется как норма вектора [2], являющегося результатом их разности:
√
d(⃗u, ⃗v ) = ∥⃗u − ⃗v ∥ = (u1 − v1 )2 + · · · + (uN − vN )2 .
(1.18)
Подобное определение расстояния позволяет считать его метрикой рассматриваемого пространства, а само пространство метрическим, поскольку оно удовлетворяет трем необходимым для этого условиям:
2.1 d(⃗u, ⃗v ) = 0 тогда и только тогда, когда ⃗u = ⃗v ;
2.2 d(⃗u, ⃗v ) = d(⃗v , ⃗u);
2.3 d(⃗u, ⃗v ) ⩽ d(⃗u, ⃗z) + d(⃗z, ⃗v ) –– так называемое правило треугольника.
11
3. Коэффициент корреляции. Скалярное произведение двух векторов ⃗u и
⃗v также можно определить как произведение норм векторов ⃗u и ⃗v и угла
между ними ∠(⃗u, ⃗v ):
(⃗u, ⃗v ) = ∥⃗u∥∥⃗v ∥ cos (∠(⃗u, ⃗v )) .
(1.19)
Совпадение формул (1.12) и (1.19) можно легко доказать с помощью теоремы косинусов:
∥⃗u − ⃗v ∥2 = ∥⃗u∥2 + ∥⃗v ∥2 − 2∥⃗u∥∥⃗v ∥ cos(∠(⃗u, ⃗v )).
(1.20)
Использование подобного определения скалярного произведения важно
при анализе схожести сигналов. О степени схожести или различия двух
сигналов, представленных в векторной форме можно судить по значению
косинуса угла между векторами, описывающими эти сигналы. Значение
косинуса лежит в пределах [−1; 1]. В случае, когда векторы, описывающие
оба сигнала максимально разнонаправлены, т.е. находятся под углом 180◦
друг к другу, значение косинуса равно −1. Если векторы имеют одно и то
же направление, то косинус угла между ними равен 1, что указывает на
максимальное с точностью до отношения норм совпадение между векторами. Косинус ортогональных векторов равен нулю. Таким образом, косинус угла между двумя векторами, описывающими сигнал, можно интерпретировать как коэффициент корреляции ru,v между двумя сигналами,
описываемыми векторами ⃗u и ⃗v [2]:
N
∑
ru,v = cos(∠(⃗u, ⃗v )) =
ui vi
(⃗u, ⃗v )
= √ i=1 √
,
∥⃗u∥∥⃗v ∥
N
N
∑
∑
u2i
vi2
i=1
(1.21)
i=1
что не противоречит определению коэффициента корреляции, используемого в теории вероятностей и математической статистике для случайных величин с нулевым математическим ожиданием или выборок с нулевым средним. Аналогичное определение коэффициента корреляции можно привести для функций:
∫b
ru,v = √
∫b
a
12
u(t)v(t)dt
√
.
b
∑
u2 (t)dt
v 2 (t)dt
a
a
(1.22)
4. Функция взаимной корреляции. В англоязычной литературе также
используется термин кросс-корреляции (cross-correlation). С помощью
функции взаимной корреляции можно оценить степень схожести двух векторов, один из которых сдвинут относительно другого на τ элементов:
R⃗u,⃗v (τ ) =
N
−τ
∑
1
⃗uk⃗vk+τ .
N −τ
(1.23)
k=1
Всего в суммировании может участвовать N − τ элементов векторов. Однако, если считать векторы периодическими последовательностями, то
суммирование можно выполнить по всем N элементам:
R⃗u,⃗v (τ ) =
N −1
1 ∑
⃗uk+1⃗v((k+τ )modN )+1 .
N
(1.24)
k=0
В пространстве функций аналогичное определение функции взаимной
корреляции выглядит следующим образом:
1
Ru,v (τ ) =
b−a
∫b
u(t)v(t + τ )dt.
(1.25)
a
Следует отметить, что интервал суммирования/интегрирования должен
быть достаточно большим, для того, чтобы получить точное представление о степени похожести двух функций.
5. Функция автокорреляции. Функция Ru,u (τ ) демонстрирует степень похожести исходного сигнала u(t) и его сдвинутых на время τ копий. Основные свойства функции автокорреляции:
– Ru,u (τ ) = Ru,u (−τ ) для вещественного сигнала u(t);
– max{|Ru,u (τ )|} ⩽ ∥u∥2 ;
τ
– arg max{Ru,u (τ )} = 0.
τ
6. Свертка. Свертка двух функций u(t) и v(t) определяется формулой:
∫∞
(u ∗ v)(t) =
∫∞
u(τ )v(t − τ )dτ =
−∞
v(τ )u(t − τ )dτ.
(1.26)
−∞
13
В дискретном случае для двух бесконечных последовательностей свертка
описывается формулой:
(u ∗ v)(n) =
∞
∑
u[i]v[n − i]
(1.27)
i=−∞
7. Линейная комбинация. Линейная
⃗v (1) , . . . , ⃗v (N ) определяется как:
⃗u =
N
∑
комбинация
αi⃗v (i) ,
из
N
векторов
(1.28)
i=1
где αi являются скалярами. Результат линейной комбинации ⃗u принадлежит тому же линейному пространству, к которому принадлежат векторы
⃗v (i) .
8. Базис. Любые N произвольных линейно-независимых векторов
{⃗v (1) , . . . , ⃗v (N ) } образуют базис N -мерного векторного линейного
пространства, т.е. любой вектор из этого пространства можно представить в виде линейной комбинации базисных векторов ⃗v (i) . В N -мерном
векторном пространстве можно выбрать бесконечное количество базисов.
Если все базисные векторы взаимно ортогональны и имеют единичную
норму, то базис является ортонормированным.Убедиться в том, что набор векторов образует ортонормированный базис легко. Для этого нужно
найти все возможные попарные скалярные произведения векторов. Базис
будет являться ортонормированным, если будет выполнено следующее
условие:
{
0, если i ̸= j,
(i) (j)
(1.29)
(⃗v , ⃗v ) =
1, если i = j.
Чтобы определить коэффициент αi , используемый для базисного вектора ⃗v (i) в линейной комбинации для произвольного вектора ⃗u, достаточно
найти скалярное произведение вектора ⃗u и соответствующего вектора из
ортонормированного базиса ⃗v (i) :


N
N (
)
∑
∑
(⃗v (i) , ⃗u) = ⃗v (i) ,
αj ⃗v (j)  =
⃗v (i) , αj ⃗v (j) = αi .
(1.30)
j=1
14
j=1
9. Правило Парсеваля. Для разложения вектора по ортонормированному
базису справедливо правило Парсеваля: квадрат нормы вектора равен сумме квадратов коэффициентов разложения по ортонормированному базису:
N
N
∑
∑
(⃗u, ⃗u) =
|αi |2 =
(⃗u, ⃗v (i) )2 .
(1.31)
i=1
i=1
10. Формулы Эйлера. Формулы Эйлера используются в Фурье анализе для
установления взаимосвязи тригонометрических функций синуса и косинуса с комплексной экспонентой на комплексной плоскости.
e±jθ = cos(θ) ± j sin(θ),
(1.32)
cos(θ) =
ejθ + e−jθ
,
2
(1.33)
sin(θ) =
ejθ − e−jθ
,
2j
(1.34)
где знаком j обозначается мнимая единица.
1.5. Фурье анализ
1.5.1. Общие сведения
Как в теоретическом анализе, так и на практике в цифровой обработке сигналов широко используются ряды и интегралы Фурье. Разложение в
ряд Фурье, где базисными функциями являются синусы, косинусы или комплексные экспоненты, применяется для анализа непрерывных периодических сигналов. Для непрерывных сигналов, имеющих конечную длительность (непериодических), в Фурье анализе на основе интегралов Фурье применяется (интегральное/непрерывное) преобразование Фурье. Однако наибольшее распространение на практике получило дискретное преобразование Фурье, применяемое для дискретных периодических сигналов. Широкое практическое применение дискретного преобразования Фурье обусловлено существованием быстрых алгоритмов, позволяющих эффективно реализовывать это преобразование с помощью аппаратных средств сигнальных
процессоров и программируемых логических матриц.
15
Для существования разложения периодической непрерывной функции
u(t) в ряд Фурье необходимо выполнение условий Дирихле [3] на периоде
функции:
– функция не должна иметь разрывов второго рода;
– число разрывов первого рода должно быть конечным;
– число экстремумов функции должно быть конечным.
1.5.2. Синусно-косинусная форма ряда Фурье
Ортонормированный базис ряда Фурье в синусно-косинусной форме,
представляет собой систему функций [2]:
√
√
√
√
(1.35)
{1, 2 cos(t), 2 cos(2t), . . . , 2 sin(t), 2 sin(2t), . . . }
Если функция u(t) имеет период 2π, то ее разложение в ряд Фурье, представленный в синусно-косинусной форме, на интервале [−π, π] имеет вид:
∞
u(t) =
a0 ∑
+
(ak cos(kt) + bk sin(kt)) ,
2
(1.36)
k=1
где коэффициенты ak , bk для k-ой гармоники вычисляются по следующим
формулам:
1
ak =
π
bk =
1
π
∫π
u(t) cos(kt)dt,
−π
∫π
u(t) sin(kt)dt.
(1.37)
−π
Разложение (1.36) соответствует формуле линейной комбинации (1.28),
где в качестве базисных используются функции из (1.35). Способ вычисления коэффициентов ak или bk соответствует формуле (1.30), которая использует скалярное произведение разложимой в ряд Фурье функции u(t) и косинуса или синуса соответствующей частоты k.
При обобщении (1.36) на произвольный период T , в разложение добавляется шкалирующий коэффициент:
ω1 = 2π/T,
16
(1.38)
как отношение значений старого периода к новому. ω1 также называют круговой частотой. В общем виде разложение функции u(t) с произвольным
периодом T выглядит следующим образом [2, 3]:
)
(
))
(
∞ (
a0 ∑
2π
2π
u(t) =
,
+
kt + bk sin
kt
ak cos
2
T
T
(1.39)
k=1
2
ak =
T
(
T /2
∫
u(t) cos
)
2π
kt dt,
T
−T /2
2
bk =
T
(
T /2
∫
u(t) sin
)
2π
kt dt.
T
(1.40)
−T /2
1.5.3. Разложение в комплексный ряд Фурье
В качестве ортонормированного базиса для комплексной формы ряда
Фурье используют систему комплексных экспонент:
{. . . , e−j2t , e−j1t , 1, ej1t , ej2t , . . .} = {ejkt },
(k = 0, ±1, ±2, . . .). (1.41)
Осуществить переход от синусно-косинусной формы к комплексной можно осуществить, воспользовавшись формулами Эйлера (1.32)-(1.34). Разложение в комплексный ряд Фурье для функции с произвольным периодом T
осуществляется с использованием формулы круговой частоты (1.38):
f (t) =
∞
∑
Ċk ejω1 kt ,
(1.42)
f (t)e−jω1 kt dt.
(1.43)
k=−∞
1
Ċk =
T
T /2
∫
−T /2
Обратите внимание, что при вычислении скалярного произведения в формуле (1.43) используется комплексно-сопряженная функция комплексной экспоненты
( jω1 kt )∗
e
= e−jω1 kt .
17
Комплексный коэффициент Ċk часто представляют в форме амплитуды
|Ċk | и фазы ϕk = ∠Ċk :
Ċk = |Ċk |ejϕk t .
(1.44)
Амплитуда и фаза позволяют связать Ċk с коэффициентами синуснокосинусной формы:
√
a2k + b2k
|Ċk | =
,
(1.45)
2
( )
bk
∠Ċk = Z + arctan
,
(1.46)
ak
где
Z=







0, если ak > 0
π, если ak < 0, bk > 0
,
−π, если ak < 0, bk < 0
которые легко получить из представления на комплексной плоскости, воспользовавшись теоремой Пифагора.
Множество коэффициентов |Ċk | называют амплитудным спектром, а
множество ∠Ċk –– фазовым спектром сигнала u(t). Множество квадратов
амплитуд |Ċk |2 принято называть спектром мощности.
1.5.4. Некоторые свойства рядов Фурье
1. Коэффициенты ряда Фурье при разложении интегрируемых функций
убывают с увеличением k:
ak , bk −−−−→ 0.
k→∞
(1.47)
Данное свойство следует из теоремы Римана (лемме Римана — Лебега)
[4], согласно которой для абсолютно интегрируемой функции u(t) на конечном или бесконечном интервале (a, b) справедливо:
∫b
∫b
u(t) cos(ωt)dt = lim
lim
ω→∞
u(t) sin(ωt)dt = 0.
ω→∞
a
(1.48)
a
2. При разложении в ряд Фурье четной функции u(−t) = u(t), коэффициенты при синусах bk равны нулю. Если пределы интегрирования
18
[−T /2; T /2] симметричны относительно нуля, то значение коэффициента ak равно удвоенному интегралу в пределах [0; T /2]. Ċk является чисто
действительной величиной.
3. При разложении в ряд Фурье нечетной функции u(−t) = u(t) коэффициенты при косинусах ak равны нулю. Ċk является чисто мнимой величиной.
4. При разложении вещественной функции в ряд Фурье в комплексной форме коэффициенты разложения при частотах −k и k являются комплексно
сопряженными числами:
∗
Ċk = Ċ−k
.
(1.49)
Из формулы (1.49) следует, что амплитудный спектр вещественной функции является симметричным относительно нуля:
|Ċk | = |Ċ−k |,
(1.50)
а фазовый спектр на частотах −k и k имеет противоположные знаки:
∠Ċk = −∠Ċ−k .
(1.51)
1.6. Интегральное преобразование Фурье
1.6.1. Переход от рядов Фурье к интегральному преобразованию Фурье
Чтобы перейти от Фурье анализа непрерывного периодического сигнала к случаю непериодического сигнала достаточно рассмотреть ситуацию,
когда период сигнала устремлен к бесконечности. Эффект от увеличения периода легко продемонстрировать на примере разложения в ряд Фурье последовательности прямоугольных импульсов длительности τ с периодом T ,
которая определяется следующим образом:
{
A, если |t| ⩽ τ /2;
u(t) =
0, если τ /2 < |t| ⩽ T.
Поскольку функция u(t) является четной, коэффициенты bk будут равны нулю, а коэффициенты ak вычисляются по формуле:
(
)
(
)
2A
πkτ
2Aτ
πkτ
ak =
sin
=
sinc
.
(1.52)
πk
T
T
T
19
Рассмотрим особенности спектра амплитуд |Ċ|, найденного для u(t) и
изображенного на рис. 1.4. Высота главного лепестка sinc-функции равна
(Aτ )/T . Частоты гармоник, для которых вычисляются значения ak , кратны 2π/T , поэтому расстояние между соседними значениями спектра также
равно этому значению при использовании в качестве оси абсцисс круговой
частоты ω = 2πk/T . Ширина главного лепестка определяется аргументом
sinc-функции. Поскольку sinc-функция принимает нулевое значение, когда
значение аргумента, кратно π, т.е. когда k кратно T /τ . Таким образом, по угловой частоте ω ширина каждого лепестка рассматриваемого амплитудного
спектра равна 2π/τ и не зависит от периода.
а)
б)
в)
Рис. 1.4.
г)
Пример изменения амплитудного спектра для последовательности прямоугольных импульсов при удвоении их периода следования: а–– последовательность прямоугольных импульсов с периодом T1 ; б–– амплитудный спектр последовательности прямоугольных импульсов с периодом T1 ; в–– последовательность прямоугольных импульсов с периодом T2 = 2T1 ;
г–– амплитудный спектр последовательности прямоугольных импульсов с периодом T2 = 2T1
20
Из формулы (1.52) видно, что увеличение периода сигнала в M раз приводит к уменьшению амплитуды коэффициентов спектра во столько же раз.
При этом количество гармоник, входящих в “лепесток” sinc-функции и определяющих его ширину увеличится в M раз, поскольку ширина лепестка не
зависит от . Другими словами формальное расстояние между гармониками,
определяемое по формуле (2π/T ) также сократится в M раз. Сигнал перестает быть периодическим при устремлении T к бесконечности. При этом
последовательность гармоник превращается в функцию от непрерывного аргумента ω, который заменяет kω1 .
Для формального описания прямого преобразования Фурье функции
u(t) часто используют обозначение F{u(t)}, а для обратного F −1 {U̇ (ω)}.
Формулы прямого и обратного преобразования Фурье будут выглядеть соответственно:
∫∞
U̇ (ω) = F{u(t)} =
u(t)e−jωt dt,
(1.53)
−∞
u(t) = F
−1
1
{U̇ (ω)} =
2π
∫∞
U̇ (ω)ejωt dω.
(1.54)
−∞
Функцию U̇ (ω) называют спектральной функцией или спектральной
плотностью. Модуль спектральной функции |U̇ (ω)| называют амплитудным спектром, а аргумент ϕ(U̇ (ω)) –– фазовым спектром.
Для существования преобразования Фурье помимо условий Дирихле
функция u(t) должна быть абсолютно интегрируема, т.е.
∫∞
|u(t)|dt < ∞.
−∞
Если использовать в преобразовании Фурье не круговую частоту ω, а частоту f = ω/ (2π), то использовать нормирующий коэффициент не нужно:
∫∞
U̇ (f ) =
u(t)e−j2πf t dt,
(1.55)
U̇ (f )ej2πf t df.
(1.56)
−∞
∫∞
u(t) =
−∞
21
1.6.2. Некоторые свойства интегрального преобразования Фурье
1. Если преобразование Фурье выполняется от вещественной функции u(t),
то:
U̇ (−ω) = U̇ ∗ (ω).
(1.57)
Отсюда следует, что для вещественной функции амплитудный спектр является четной функцией:
|U̇ (−ω)| = |U̇ (ω)|,
а фазовый спектр –– нечетной:
ϕ(U̇ (−ω)) = −ϕ(U̇ (ω)).
2. Спектральная функция вещественной четной функции является вещественной четной функцией.
3. Спектральная функция вещественной нечетной функции является чисто
мнимой нечетной функцией.
4. Преобразование Фурье является линейным:
F{α1 u1 (t) + α2 u2 (t)} = α1 F{u1 (t)} + α2 F{u2 (t)}.
(1.58)
5. Сдвиг сигнала u(t) на величину τ приводит к сдвигу фазы спектральной
функции:
F{u(t − τ )} = F{u(t)}e−jωτ .
(1.59)
Амплитудный спектр при этом остается без изменений.
6. При масштабировании временной оси спектр сигнала u(at) выглядит следующим образом:
∫∞
F{u(at)} =
u(at)e
−jωt
−∞
1
dt =
a
z=at
∫∞
u(z)e−j a z dz =
ω
−∞
1 (ω)
U̇
.
|a|
a
(1.60)
7. Преобразование Фурье свертки двух функций равно произведению их
спектров:
F{(u ∗ g)(t)} = U̇ (f )Ġ(f ).
(1.61)
8. Правило Парсеваля.
∫∞
∫∞
2
u (t)dt =
−∞
22
1
|U̇ (f )| df =
2π
∫∞
|U̇ (ω)|2 dω.
2
−∞
−∞
(1.62)
1.7. Постановка задачи аналого-цифрового и цифро-аналогового
преобразования
Сигналы, формируемые детекторами и датчиками, а также сигналы, используемые устройствами воспроизведения имеют непрерывную природу.
В качестве наглядного примера таких устройств можно привести акустический микрофон и динамик. В состав обоих устройств входит мембрана. Колебания воздуха, формируемые при передаче звука, передаются на мембрану.
Колебания мембраны тем или иным способом преобразуются в напряжение,
изменение которого во времени можно описать функцией от непрерывного аргумента u(t). В свою очередь для того, чтобы воспроизвести звук, зафиксированный с помощью мембраны микрофона, необходимо, используя
напряжение u(t), организовать аналогичные колебания воздуха, что выполняется с помощью мембраны динамика. Хранить полностью изменения амплитуды u(t) во времени можно на аналоговом носителе информации, например магнитной пленке. Однако такой способ является весьма требовательным к физическому объему носителя. Современные средства хранения
информации используют цифровое представление un , формируемое на базе
u(t), которое позволяет значительно уменьшить физические габариты носителя. Однако при этом происходит частичная потеря информации об исходной функции напряжения u(t). Процесс формирования цифрового сигнала
un по аналоговому сигналу u(t) называют аналого-цифровым преобразованием (АЦП), а процедуру восстановления аналогового сигнала u′ (t) по un ––
цифро-аналоговым преобразованием (ЦАП).
АЦП можно структурно разделить на две процедуры (3.1):
1. дискретизация (дискретизация по времени),
2. квантование (дискретизация по уровню).
Простейшей вариант ЦАП (широтно-импульсный модулятор) структурно
состоит из следующих процедур:
1. формирование по цифровым данным первичного кусочно-константного
аналогового сигнала, амплитуда которого постоянна на интервале
[nT, (n + 1)T ) и равна аппроксимирующему значению, соответствующему un ,
2. пропускание первичного аналогового сигнала через восстанавливающий
фильтр, который на практике реализуется с помощью низкочастотного
23
фильтра с частотой среза обратно пропорциональной периоду дискретизации (подробнее будет рассмотрено ниже в подразделах 1.8 и 1.9).
В данном разделе будут подробно изучаться процедура дискретизации
аналогового сигнала по времени, а также процедуры восстановления аналогового сигнала в устройстве ЦАП.
1.8. Дискретизация аналогового сигнала. Теорема Котельникова
В данном разделе будут рассмотрены теоретические аспекты наиболее
популярного способа дискретизации по времени –– равномерная дискретизация. При переходе от аналогового представления сигнала u(t), имеющего
спектр U (f ), к дискретным отсчетам u[n] при равномерной дискретизации
необходимо определить период дискретизации Td или частоту дискретизации fd , которые обратно пропорциональны друг другу:
fd =
1
.
Td
(1.63)
Для класса сигналов с ограниченным спектром существуют такие значения
периода (частоты) дискретизации, при которых возможно точное восстановление исходного сигнала u(t) по его дискретным отсчетам u[n]. Взаимосвязь
значения периода Td (частоты fd ) дискретизации со значением максимальной ненулевой частоты в спектре исходного сигнала устанавливает теорема, сформулированная независимо в работах Котельникова [5] в 1933 году
и Шеннона [6] в 1949 году. В русскоязычной литературе она получила распространение под названием теоремы Котельникова. В англоязычной литературе ее называют теоремой Шеннона или теоремой отсчетов (sampling
theorem).
Теорема. Пусть аналоговый сигнал u(t) имеет ограниченный частотой
fc спектр, т.е. U (f ) = 0 при |f | > fc . Тогда восстановление исходного аналогового сигнала возможно, если равномерная дискретизация выполняется
с частотой fd , минимум вдвое превышающую частоту fc :
fd > 2fc ,
(1.64)
1
.
2fc
(1.65)
или
Td <
24
В цифровой обработке сигналов также используют понятие частоты Найквиста fN , которая равна половине частоты дискретизации fd . Т.о. теорему
Котельникова можно также сформулировать следующим образом: исходный
сигнал u(t) можно восстановить со сколь угодно высокой точностью, если
его максимальная частота меньше частоты Найквиста.
Для доказательства теоремы рассмотрим спектр дискретизуемого сигнала U̇ (f ), который можно получить, применив прямое преобразование Фурье
к сигналу u(t) [7]. Если |f | < fc , то пределы интегрирования в обратном преобразовании Фурье можно сократить до [−fc ; fc ]. В то же время, поскольку
спектр сигнала U̇ (f ) ограничен значением fc , его можно продолжить периодически с периодом 2fc , получив функцию Ũ (f ) (рис.1.5). Непрерывную
периодическую функцию Ũ (f ) можно разложить в ряд Фурье:
Ũ (f ) =
∞
∑
Ċk e
j2πf k
2fc
,
(1.66)
k=−∞
комплексные коэффициенты которого вычисляются по формуле:
1
Ċk =
2fc
∫∞
Ũ (f )e
−j2πf k
2fc
df.
(1.67)
−∞
.
|U(f)|
-fd
-fd/2 -fc
fc fd/2
fd
f
Рис. 1.5. Амплитудный спектр дискретизированного сигнала
Теперь определим значение спектра функции u(kTd ), значения которой
отличны от нуля только в точках kTd временной оси, и положим Td = 2f1c .
25
Тогда
(
u(kTd ) = u
k
2fc
)
∫∞
=
U̇ (f )e
j2πf k
2fc
df.
(1.68)
−∞
Сравнивая формулы (1.67) и (1.68), нетрудно увидеть соотношение
(
)
1
−k
Ċk =
u
.
2fc
2fc
(1.69)
Таким образом, зная значения функции u(t) в точках u(kTd ), можно определить значения коэффициентов Ċk , воспользовавшись формулой (1.69). С
помощью коэффициентов Ċk можно определить функцию Ũ (f ). Зная Ũ (f ),
можно определить U̇ (f ) и, наконец, восстановить u(t):
∫fc
U̇ (f )ej2πf t df =
u(t) =
∫fc
−fc
[
∞
∑
=
Ċk e
]
k
j 2πf
2fc
ej2πf t df =
(1.70)
k=−∞
−fc
∫fc
∞
∑
=
Ċk
k=−∞
ej2πf ( 2fc +t) df.
k
−fc
Теперь воспользуемся формулой Эйлера (1.32) и обратим внимание, что
интервал интегрирования [−fc ; fc ] симметричен относительно начала координат. Поскольку функция синуса является нечетной, то значение интеграла
будет равно нулю. Для четной функции косинуса результат интегрирования
на интервале [−fc ; fc ] равен удвоенному интегралу в пределах неотрицательной части [0; fc ]. Таким образом выражение (1.70) можно записать так
(
)
k
u(t) = 2
Ċk cos 2πf
+ t df =
2fc
k=−∞
0
(
)

∞
sin 2πfc 2fkc + t
∑
(
) .
=
Ċk 
π 2fkc + t
k=−∞
∞
∑
26
∫fc
(1.71)
С учетом формулы (1.67) и, положив Td =
ному выражению:
u(t) =
1
2fc ,
перейдем к окончатель-
]
sin 2πfc (kTd + t)
=
π(kTd + t)
k=−∞
[
]
∞
∑
sin π (kTTdd+t)
=
=
Ċk
π(kT + t)
k=−∞
[
]
∞
d)
∑
sin π (t−nT
Td
=
u(nTd )
,
d)
π (t−nT
n=−∞
Td
∞
∑
[
Ċk
(1.72)
где n = −k. Таким образом доказано, что для точного восстановления сигнала, описываемого непрерывной функцией u(t), спектр которого ограничен
частотой fc , достаточно формировать дискретные отсчеты u[n] с периодом
Td < 2f1c , меньшим значения обратно пропорционального удвоенной максимальной частоте исходного сигнала.
Функцию sin(x)
= sinc (x) называют sinc-функцией, а ряд, полученный
x
в формуле (1.72), называют интерполяционным рядом Котельникова:
( (
))
t
sinc π
−n
,
(1.73)
Td
поскольку с помощью него по узловым точкам u(nTd ) можно восстанавливать промежуточные значения функции u(t). В западной литературе
формулу (1.72) называют интерполяционной формулой Уитекера-Шеннона
(Whittaker–Shannon).
1.9. Восстановление аналогового сигнала по дискретным отсчетам
Выражение (1.73) по сути представляет собой дискретную свертку сигнала u(nTd ) и идеального низкочастотного фильтра с импульсным откликом h(t) = sinc(π Ttd ). Спектр импульсного отклика будет представлять собой
импульс единичной высоты, ограниченный интервалом
[ прямоугольный
]
π
π
− Td ; Td . Импульсный отклик фильтра и амплитуда его Фурье-образа изображены на рис. 1.6. Данный фильтр невозможно реализовать на практике,
поскольку фильтр имеет бесконечно протяженный импульсный отклик, поэтому вместо него используют аппроксимации, имеющие близкую к прямо27
а)
б)
Рис. 1.6. Примеры характеристик идеального низкочастотного фильтра: а–– импульсный
отклик, имеющий форму sinc(π Tt ); б–– амплитудный спектр
d
угольной амплитудно-частотную характеристику (подробно этот вопрос будет рассмотрен в следующих разделах). Процедуру восстановления аналогового сигнала по последовательности дискретных отсчетов можно описать
в виде следующих формальных шагов.
1. Формирование на основе последовательности u[n] сигнала un (t), значение которого постоянно на интервале [nTd , (n + 1)Td ) и равно аппроксимирующему значению для цифрового сигнала un .
2. Фильтрация сигнала un (t) с помощью низкочастотного фильтра с частотой среза π/Td .
2. Порядок самостоятельного выполнения исследования
В ходе выполнения исследования необходимо выполнить следующие задания.
1. Выполнить аналитический расчет спектров для двух классов сигналов
(периодических и непериодических). В отчете представить графики исходных сигналов и полученных спектров, а также описать их характерные
особенности.
1.1 Периодические сигналы с периодом T :
1.1.1. Последовательность прямоугольных импульсов:
{
A, |t| ≤ τ /2
u(t) =
,
(1.74)
0, |t| > τ /2
28
где τ < T .
1.1.2. Пилообразный сигнал:
u(t) =
2A
1
1
(t − kT ), (k − )T < t ≤ (k + )T.
T
2
2
(1.75)
1.1.3. Последовательность треугольных импульсов:
u(t) = A(1 − 4
|t − kT |
1
1
), (k − )T ≤ t < (k + )T. (1.76)
T
2
2
1.1.4. Косинус:
u(t) = A cos(2πf0 t + θ).
(1.77)
u(t) = A sin(2πf0 t + θ).
(1.78)
1.1.5. Синус:
1.2 Непериодические сигналы:
1.2.1. Односторонний экспоненциальный импульс:
{
Ae−at , t ≥ 0
u(t) =
.
0, t < 0
1.2.2. Прямоугольный импульс:
u(t) =
{
A, |t| ≤ τ /2
0, |t| > τ /2
.
1.2.3. Несимметричный треугольный импульс:

 At, 0≤t≤T
T
u(t) =
.
 0, t < 0, t > T
1.2.4. Задержанный прямоугольный импульс:
{
A, 0 ≤ t ≤ τ
.
u(t) =
0, t < 0, t > τ
1.2.5. Симметричный треугольный импульс:
(
)

 A 1 − |t| , |t| ≤ T
T
u(t) =
.

0, |t| > T
(1.79)
(1.80)
(1.81)
(1.82)
(1.83)
29
1.2.6. Гауссовский импульс:
u(t) = Ae−t
2
/a2
.
(1.84)
2. Проанализировать свойства двух функций u1 (t) = sin(2πf1 t) и u2 (t) =
sin(2πf2 t) а интервале [−T /2; T /2]. Написать программу, которая позволит:
2.1 Вычислить все значения функций u1 (t) и u2 (t) на заданном интервале с шагом 10−3 .
2.2 Вычислить приближенное значение скалярного произведения двух
функций (u1 (t), u2 (t)), воспользовавшись аппроксимацией формулы (1.15).
2.3 Вычислить нормы обеих функций.
2.4 Определить, являются ли исходные функции ортогональными друг
к другу.
2.5 Как нужно изменить исходные функции, что они могли являться
элементами ортонормированного базиса? Выполните данную модификацию и продемонстрируйте результат.
2.6 Останутся ли исследуемые функции элементами ортонормированного базиса, если:
2.6.1. частоты f1 и f2 удвоятся;
2.6.2. интервал [−T /2; T /2] увеличится вдвое;
2.6.3. интервал [−T /2; T /2] уменьшится вдвое.
Объясните почему.
3. Исследовать процедуру дискретизации синусоидального сигнала u(t) с
частотой F Гц. Длительность наблюдения сигнала 10/F секунд. Написать программу, которая позволит:
3.1 сформировать выборки отсчетов u(i) [n] (результаты дискретизации)
(i)
исследуемого сигнала с частотами дискретизации fd равными:
3.1.1.
3.1.2.
3.1.3.
3.1.4.
3.1.5.
30
1.5F Гц,
1.75F Гц,
2F Гц,
3F Гц,
1000F Гц.
3.2 применить ряд Котельникова (1.73) для восстановления исходного
сигнала по его дискретным отсчетам с помощью формулы (1.72).
3.3 вывести графики исходного сигнала и восстановленных сигналов
u(i) (t), где i –– индекс частоты дискретизации.
Сделать выводы о точности восстановления исследуемого сигнала u(t)
при использовании различных частот дискретизации.
Объяснить, чем вызваны искажения на краях построенных графиков.
4. Реализовать процедуру передискретизации изображения с помощью интерполяционного ряда Котельникова. Формат исходного изображения ––
BMP24, разрешение исходного изображения W × H пикселей. Результатом передискретизации будет изображение размером nW × mH. Выполните передискретизацию с различными комбинациями значений m и n:
4.1 m > 1, n > 1,
4.2 m < 1, n < 1,
4.3 m > 1, n < 1,
4.4 m < 1, n > 1.
31
II. ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ
ФУРЬЕ
Цели исследования: изучение методов Фурье-анализа дискретных и
цифровых сигналов.
1. Теоретические основы
1.1. Дискретное преобразование Фурье
1.1.1. Переход от рядов Фурье к дискретному преобразованию Фурье
Дискретное преобразование Фурье является наиболее популярным математическим аппаратом в Фурье-анализе, поскольку оно применимо к цифровым сигналам и для данного преобразования существует множество вариантов быстрых реализаций, которые реализуются в сигнальных процессорах.
Дискретное преобразование Фурье, применяемое для дискретных периодических сигналов легко получить, используя разложение в ряд Фурье непрерывной периодической функции u(t), из которой рассматривается только N значений un , где n = 0, . . . , N − 1. Здесь намеренно используется обозначение цифрового сигнала, поскольку на практике преобразование Фурье выполняется именно для такого класса сигналов, хотя применение
преобразования Фурье для дискретных отсчетов также является справедливым. Итак, разложение в ряд Фурье сигнала u(t) осуществляется по формуле
(1.42). При вычислении разложения для фиксированных значений un с шаn
гом N
T будет выглядеть следующим образом:
un =
∞
∑
Ċk ej
2πn
N k
.
(2.1)
k=−∞
Легко выполнить дискретизацию исходного сигнала таким образом, чтобы
он остался периодическим, т.е. un+mN = un для любого целого m. Посколь2πn
2πn
ку ej N (k+mN ) = ej N k выражение (2.1) можно переписать следующим
32
образом, получив формулу обратного дискретного преобразования Фурье:
( ∞
)
N
−1
N
−1
∑
∑
∑
2πn
j 2πn
k
N
U̇k ej N k ,
un =
Ċk+mN e
=
(2.2)
k=0
m=−∞
k=0
где n принимает значения от 0 до N − 1. Формула прямого преобразования
Фурье подразумевает процедуру вычисления значений U̇k . Для этого найдем скалярное произведение векторного представления un и вектора на базе
2πn
комплексной экспоненты ej N m и применим к нему формулу (2.2):
)
( −1
N
−1
N
−1 N
∑
∑
∑
2πn
2πn
−j 2πn
m
j
k
un e N
=
U̇k e N
e−j N m .
(2.3)
n=0
n=0
k=0
Изменим порядок суммирования и объединим показатели экспонент, получив:
N
−1
N
−1
N
−1
∑
∑
∑
2πn
U̇k
e−j N (k−m) =
(2.4)
U̇k (N δkm ) = N Um .
n=0
k=0
k=0
Объединив формулы (2.2) и (2.4) и заменив в итоговом варианте формальную
переменную m на k, получим формулу прямого дискретного преобразования
Фурье:
N −1
2πn
1 ∑
U̇k =
un e−j N k .
(2.5)
N n=0
Нормирующий коэффициент 1/N , как правило, применяется при вычислении обратного ДПФ вместо прямого.
Поскольку прямое и обратное ДПФ можно интерпретировать в терминах
операций над векторами, это преобразование удобно представлять в матричной форме соответственно:
⃗ = ⃗uFH .
U
(2.6)
⃗ F.
⃗u = U
(2.7)
Матрица F имеет вид:

2π0
ej N 0
 j 2π0
 e N 1
F=
..


.
ej
2π0
N (N −1)
2π1
ej
ej N 0
2π1
ej N 1
..
.
···
···
..
.
2π1
N (N −1)
···
2π(N −1)
ej N 0
2π(N −1)
ej N 1
..
.
ej
2π(N −1)
(N −1)
N



,


а FH является эрмитово сопряженной с F матрицей.
33
1.1.2. Некоторые свойства дискретного преобразования Фурье
Свойства дискретного преобразования Фурье схожи со свойствами ранее
рассмотренного преобразования Фурье.
1. Линейность.
(1)
(2)
α1 uk + α2 uk ↔ α1 U̇n(1) + α2 U̇n(2) .
2. Задержка.
uk−τ ↔ U̇n e−j
2πτ n
N
.
3. Симметрия спектра.
Для вещественного дискретного сигнала u[k] выполняются следующие
свойства симметрии спектра:
U̇−k = U̇N −k = U̇k∗ .
Значение спектра при нулевой частоте называют постоянной составляющей и вычисляют по формуле:
U̇0 =
N
−1
∑
uk .
(2.8)
k=0
Если размер преобразования N четный, то частоту N /2 называют частотой Найквиста и вычисляют значение спектра на этой частоте как знакопеременную функцию исходного сигнала:
U̇N /2 =
N
−1
∑
uk (−1)k .
(2.9)
k=0
4. Равенство Парсеваля.
N
−1
∑
k=0
u2k =
N −1
1 ∑
|U̇n |2 .
N n=0
(2.10)
1.2. Спектр дискретного сигнала
В данном подразделе рассматривается особенности спектра дискретного
сигнала u[n], полученного путем дискретизации из аналогового сигнала u(t),
34
т.е. значения отсчетов u[n] являются значениями исходного сигнала, взятого
с периодом T :
u[n] = u(nT ).
(2.11)
Более детально со свойствами спектра дискретного сигнала можно ознакомиться в [3].
Воспользовавшись фильтрующим свойством дельта-функции, непрерывный сигнал, сформированный из дискретных отсчетов можно описать
следующей формулой:
ud (t) =
∞
∑
u[n]δ(nT − t) = u(t)
n=−∞
∞
∑
δ(nT − t).
(2.12)
n=−∞
Сумма дельта-функций является периодической последовательностью с периодом T , которая раскладывается в ряд Фурье с коэффициентом:
T /2
∫
1
Ċk =
T
δ(t)e−j T
2π
kt
dt =
1
.
T
−T /2
Таким образом, заменив сумму из (2.12) на ее разложение в ряд Фурье получим:
(
)
∞
∞
2π
1 ∑ j 2π kt
1 ∑
ud (t) = u(t)
=
(2.13)
e T
u(t)ej T kt .
T
T
k=−∞
k=−∞
Переходя к спектру
(
)
∞
1 ∑
2πk
˙
Ud (ω) =
U̇ ω −
,
T
T
(2.14)
k=−∞
получаем, что спектр дискретного сигнала является суммой бесконечного числа копий спектра исходного сигнала U̇ (ω), сдвинутых на величину
2πk/T . На этом примере также можно наглядно продемонстрировать теорему Котельникова. Исходный сигнал может быть восстановлен в процессе фильтрации сигнала ud (t) фильтром с частотой среза π/T , как показано
в подразделе 1.9 Однако, если период дискретизации аналогового сигнала
превышает величину 1/ (2fc ), то сдвинутые копии спектра исходного сигнала наложатся друг на друга, проявив так называемый эффект элайзинга
(от английского alising –– наложение). Восстановить исходный сигнал в этом
случае в точности не удастся.
35
1.3. Растекание спектра
При определении дискретного преобразования Фурье говорилось, что
дискретизация исходного сигнала выполняется таким образом, чтобы он
остался периодическим. При этом, если значения начальных и конечных отсчетов сильно различаются, то на краях периодически повторяемых сегментов будут возникать скачки. Наличие подобных скачков на краях сегментов
приводит к явлению, называемому растеканием спектра.
Данное явление наиболее наглядно проявляется при анализе спектра
гармонического сигнала. Предположим, что в качестве анализируемого сигнала выступает дискретный сигнал следующего вида:
u(t) = A sin(ωkTd + ϕ), k = 0, 1, · · · , N − 1.
(2.15)
Если анализируемая последовательность содержит целое количество периодов (т.е. N ωTd /(wπ)–– целое число), то все спектральные коэффициенты, кроме двух, являются нулевыми значениями:

ωTd
AN jϕ


2 e , n = 2π N,

(
)
AN −jϕ
d
(2.16)
U̇k =
, n = 1 − ωT
2 e
2π N,


 0, иначе.
Такой набор спектральных коэффициентов аналогичен спектру
непрерывного гармонического сигнала. В альтернативном случае, когда
N ωTd /(wπ) не является целым числом, спектр имеет большее количество
ненулевых значений. Пример амплитудного спектра для подобной последовательности представлен на рисунке 2.1. Подобное поведение может
объяснятся тем, что при периодическом повторении данного набора значений получаемая последовательность не будет являться набором значений
непрерывной синусоиды.
Для уменьшения эффекта растекания спектра дискретного преобразования Фурье используются весовые функции, также известные как окна. Перед
произведением дискретного преобразования Фурье обрабатываемый сигнал
умножают на вещественную неотрицательную весовую функцию, значения
которой убывают на краях сегмента. Тогда формула прямого преобразования
Фурье с использованием весовой функции имеет следующий вид:
U̇k =
36
N −1
2πn
1 ∑
un vn e−j N k .
N n=0
(2.17)
а)
б)
Рис. 2.1. Пример дискретного преобразования Фурье для нецелого числа периодов гармонческого сигнала: а–– графики аналогового и дискретизированного исходного сигнала вида sin();
б–– амплитудный спектр анализируемого сигнала
На практике наиболее часто применяются следующие весовые функции.
1. функция окна Ханна (Хэннинга):
(
(
))
2πn
v(n) = 0.5 1 − cos
.
N −1
(2.18)
2. функция окна Хэмминга:
(
v(n) = 0.53836 − 0.46164 cos
2πn
N −1
)
.
(2.19)
Примеры других окон можно найти в [3].
2. Порядок самостоятельного выполнения исследования
В ходе выполнения исследования дискретного преобразования Фурье
необходимо выполнить следующие задания.
37
1. Написать программу вычисления прямого и обратного дискретного преобразования Фурье в матричной форме.
2. Продемонстрировать с помощью написанной программы свойства линейности, сдвига сигнала во времени и равенство Парсеваля.
3. Произвести декодирование аудио-файла с записью тонального сигнала
(Dual-Tone Multi-Frequency (DTMF)) сигнала в формате WAV PCM 16 bit,
mono. Данный способ кодирования предполагает, что кодируемое значение представляется в виде пары различных частот (f1 /f2 ) в соответствии
с приведенной таблицей 2.1. Затем, сигнал представляется в виде отсчетов суммы двух синусоид соответствующих частот. Для декодирования
сигнала необходимо произвести прямое дискретное преобразование для
имеющегося набора отсчетов и определить частоты используемых для кодирования синусоид по полученному амплитудному спектру.
Таблица 2.1.: Таблица частот для кодирования DTMF сигналов
(f1 /f2 ) 1209 Гц
697 Гц
1
770 Гц
4
852 Гц
7
*
941 Гц
1336 Гц
2
5
8
0
1477 Гц
3
6
9
#
1633 Гц
A
B
C
D
4. Получить у преподавателя два файла с изображениями в формате BMP24.
Выполнить оценку смещения между двумя изображениями путем анализа
фазового спектра.
38
III. МЕТОДЫ КВАНТОВАНИЯ В
СИСТЕМАХ ОБРАБОТКИ СИГНАЛОВ
Цель исследования: изучение методов квантования, используемых в системах цифровой обработки сигналов.
1. Теоретические основы
1.1. Классификация методов квантования
В широком смысле под квантованием понимают правило отображения
элементов x из некоторого исходного множества X в элементы другого множества Y :
y = Q(x),
(3.1)
причем процедура отображения такова, что несколько элементов из X (для
непрерывного множества это подмножество значений) могут отображаться
в один и тот же элемент из Y . Подмножество из X, все элементы которого
отображаются в один и тот же элемент из Y , образуют квант. Множество
Y принято называть аппроксимирующим множеством. С точки зрения теории множеств, множество X мощнее Y . Как правило, множество Y является
конечным или счетным.
Конечный результат процедуры квантования зависит от конкретной прикладной задачи, в которой она применяется. Это может быть аппроксимирующий элемент y или номер кванта, в который попадает x. Результат представляется в виде индекса или кода, удобного для дальнейшего хранения или
передачи в цифровой форме.
Пример 1. Применение квантования в цифровой технике.
Квантование является одним из методов, применяемых для перевода аналогового сигнала в цифровую форму. Процедура аналого-цифрового преобразования схематично изображена на рис. 3.1. С помощью дискретизации аналоговый сигнал преобразуется в дискретный, после чего дискретные отсчеты преобразуются в значения цифрового сигнала с помощью процедуры квантования. Квантование также называют дискретизацией сигнала по
39
X(t)
X[n]
Дискретизация
X
Xn
Квантование
Xn
X(nT)
t
t
n
Рис. 3.1. Схема аналого-цифрового преобразования
уровню, имея в виду, что по завершению процедуры значения цифрового сигнала выбираются из конечного множества, мощность которого, как правило, равна степени числа 2.
В данном исследовании будут рассмотрены три способа квантования.
1. Скалярное равномерное.
2. Скалярное неравномерное.
3. Векторное.
Квантование называют скалярным, если элементами множеств X и Y
являются числовые (скалярные) величины. Если элементы обоих множеств
являются векторами, то квантование называют векторным. В этом случае
при отображении по формуле (3.1) длины всех векторов x̄ и ȳ должны совпадать.
1.2. Критерии оценки эффективности квантования
Поскольку по определению множество X мощнее аппроксимирующего
множества Y , квантование в общем случае является преобразованием с потерями. Эффективность квантования, как правило, оценивается по ошибке,
вносимой при выполнении процедуры. Однако в некоторых прикладных задачах могут также использоваться дополнительные критерии. Например, в
задаче сжатия источника с потерями, в которой также применяется квантование, для оценки эффективности совместно с ошибкой используют степень
достигаемого сжатия.
40
Введем в рассмотрение модель дискретного источника, выход которого будет подвергаться квантованию. Данные источника представляют собой независимые, одинаково распределенные случайные величины X (будем
обозначать случайную величину по аналогии с множеством квантуемых значений, поскольку эти понятия взаимосвязаны). Также будем полагать, что
случайная величина X является непрерывной и имеет плотность распределения f (x).
Для количественного измерения ошибки квантования будем использовать метрику эвклидова расстояния. Если квантованию подвергаются векторы x̄ длины n, то в качестве расстояния ∥x̄−ȳ∥ будем использовать евклидово
расстояние в n-мерном пространстве:
√
∥x̄ − ȳ∥ ≜ (x1 − y1 )2 + . . . + (xn − yn )2 ,
(3.2)
где xi и yi определяют i-ю позицию в соответствующем векторе. Для скалярных величин x и y в качестве расстояния ∥x − y∥ будем использовать
абсолютную разность двух этих величин.
С помощью формулы (3.2) определяется ошибка квантования только для
одного элемента исходного множества X. В качестве интегрального критерия эффективности обычно используют взвешенный интеграл квадратов
ошибок, который называют средним квадратом ошибки:
∫∞
ε2X
∫∞
∥x − Q(x)∥2 f (x)dx,
2
=
e (x)f (x)dx =
−∞
(3.3)
−∞
где f (x) задает плотность распределения случайной величины, принимающей значения из X.
Также в качестве критерия оценки ошибки, вносимой при скалярном
квантовании, используют формулы отношения сигнал-шум [8], [9]:


SN R = 10 lg 

M
∑
i=1
M
∑
(x −

x2i
Q(x))2

,

(3.4)
i=1
SN R = 10 lg
2
σX
.
2
εX
(3.5)
41
Заметим, что формула (3.5) имеет тот же смысл, что и (3.4), но в качестве
2
энергии сигнала используется σX
–– дисперсия случайной величины, которая подвергается квантованию (обратите внимание, что X должна быть центрированной случайной величиной, т.е. M [X] = 0), а в качестве энергии
шума –– ε2X , значение среднего квадрата ошибки, определяемой по формуле
(3.3). Формулу (3.4) следует трактовать, как логарифм отношения энергии
сигнала, представленной в виде суммы квадратов значений квантуемой выборки, к энергии ошибки, которая выражается как сумма квадратов ошибок
квантования.
1.3. Скалярное квантование
1.3.1. Общие сведения
Поскольку элементами множества X и Y для скалярного квантования
являются числа, для простоты будем далее полагать, что множество X есть
подмножество вещественных чисел. Тогда кванты будут представлять собой
интервалы вещественной оси, а определение процедуры квантования можно интерпретировать как задание способа разбиения числовой оси на эти
интервалы и выбора в каждом из них единственного аппроксимирующего
значения. Таким образом, для того чтобы определить процедуру скалярного
квантования, необходимо задать следующие величины.
1. Число квантов N . Кванты при этом нумеруются от 1 до N . Для упрощения
кодирования номера кванта с помощью равномерного кода, число квантов
выбирают как N = 2R , где R определяет число двоичных разрядов равномерного кода.
2. Границы квантов ti , i = 0, ..., N . j-й квант определяет полуоткрытый интервал [tj−1 ; tj ). В общем случае неважно, находится закрытая граница
’[’ справа или слева. В данной работе будем полагать, что такая граница
всегда будет располагаться слева. Исключение могут составлять первый
и последний кванты. Если множество X неограничено слева, то первый
квант (−∞, t1 ) будет открытым. Если множество X ограничено справа
значением tN , то последний квант [tN −1 , tN ] будет закрытым.
3. Аппроксимирующие значения yi ∈ [t(i−1) ; ti ), i = 1, ..., N .
42
Сама процедура скалярного квантования сводится к замене вещественного
числа x на аппроксимирующее значение кванта, в который оно попадает:
x ∈ [t(i−1) , ti ) 7−→ Q(x) = yi .
Заметьте, что yi также являются вещественным числом.
Формулу (3.3) среднего квадрата ошибки можно переписать следующим
образом:
∫ti
N
∑
2
εX =
(x − yi )2 f (x)dx.
(3.6)
i=1t
(i−1)
Процедуру отображения при скалярном квантовании часто изображают
графически в виде ступенчатой функции [10] в декартовой системе координат, как это показано на рис. 3.2.
Здесь на оси аргумента отложены элементы исходного множества X, а
по оси ординат –– значения аппроксимирующего множества Y . Если множество X является неограниченным, то t0 = −∞, а tN = ∞. Для ограниченного множества X границы крайних квантов и определяются предельными
элементами множества X.
В зависимости от положения нулевого значения различают два вида скалярного квантования:
1. нулевое значение расположено на границе одного из квантов (рис. 3.3а).
В англоязычной литературе этот вид квантования называют midrise. При
этом ноль не является аппроксимирующим значением;
2. нулевое значение расположено внутри одного из квантов, часто в его середине (рис. 3.3б). В англоязычной литературе этот вид квантования называют midtread. При этом ноль, как правило, является аппроксимирующим значением. Примером такого вида скалярного квантования является
арифметическое округление.
Если все кванты, кроме крайних для случая, когда множество X неограничено, имеют одинаковый размер, то квантование называют равномерным
(в англоязычной литературе uniform scalar quantization). В других случаях
квантование называется неравномерным.
Пример 2. Равномерное скалярное квантование случайной величины, распределенной по
равномерному закону на интервале [−3, 3].
Пусть равномерному квантованию подвергаются отсчеты xt ∼ U [−3, 3], распределенные по
равномерному закону на интервале [−3, 3], и число квантов равно 8. В случае равномерного
43
xMIN y
1
t0
y2
t1
y3
t2
0
y4
t3
xMAX
t4
X
y4
xMIN
t0
xMAX
t1
t2
y3
t3
t4
X
y2
y1
Рис. 3.2. Пример скалярного квантования N = 4: сверху рисунка приведен пример разбиения числовой оси; снизу –– ступенчатая функция
квантования ∆ = 43 , ti = −3 + i∆ и yi =
описанием схемы квантования.
i
1
2
3
4
5
6
7
8
t(i−1) +ti
.
2
[t(i−1) ; ti )
[−3; −2.25)
[−2.25; −1.5)
[−1.5; −0.75)
[−0.75; 0)
[0; 0.75)
[0.75; 1.5)
[1.5; 2.25)
[2.25; 3]
Ниже представлена таблица с полным
yi
−2.625
−2
−1.125
−0.375
0.375
1.125
2
2.625
код
000
001
010
011
100
101
110
111
Рассмотренный пример относится к равномерному скалярному квантованию вида midrise.
44
а)
б)
Рис. 3.3. Виды скалярного квантования: а–– midrise – ноль на границе кванта; б–– midtread
– ноль внутри кванта
1.3.2. Анализ эффективности равномерного скалярного квантования
Равномерное скалярное квантование является популярным в силу простоты его реализации. Все кванты за исключением крайних одинаковы и равны некоторой величине ∆. Для того чтобы определить равномерное скалярное квантование, необходимо задать только число квантов N . Границы квантов рассчитываются, исходя из области значений случайной величины X.
Аппроксимирующие значения выбираются в середине квантов.
Пример 3. Округление вещественных чисел.
Наиболее известным примером равномерного скалярного квантователя является округление вещественных чисел r ∈ ℜ до целых z ∈ Z. При округлении X ∼ ℜ, Y ∼ Z и ∆ = 1. Границы
квантов совпадают со значениями [z − 0.5, z + 0.5), где z ∈ Z. Особенностью данного примера является то, что при округлении элементы из непрерывного множества преобразуются в
элементы из счетного множества. Сама процедура квантования может быть представлена в виде
арифметического выражения:
y = ⌊x + 0.5⌋,
где ⌊.⌋ –– операция отбрасывания дробной части, реализуемая на языке С с помощью функции
floor().
Размер кванта ∆ выбирается в соответствии с количеством квантов N
и диапазоном значений, принимаемых случайной величиной X. Наиболее
простую зависимость ошибки квантования от количества квантов N можно
продемонстрировать для случая, когда данные распределены по равномер45
ному закону U [a, b]. Равномерное распределение задается двумя параметрами –– минимальным a и максимальным b значениями, которые может принимать случайная величина X.
a)
f(X)
X ~ U(a, b)
y1
a
t0
y2
t1
б)
1
b-a
y3
y4
0
t2
t3
b
t4
X
e(x) = x - Q(x)
Δ/2
a
t0
y1
y2
t1
y3
t2
y4
0
t3
b
t4
X
-Δ/2
Рис. 3.4. Ошибка квантования случайной величины X ∼ U [a, b]
В качестве меры вносимых при квантовании искажений будем использовать отношение «сигнал-шум», задаваемое формулой (3.5). Дисперсия равномерно распределенной случайной величины равна:
2
σX
=
(b − a)2
(2R ∆)2
=
12
12
(3.7)
Поскольку при квантовании равномерно распределенной случайной величины вероятности ошибок квантования внутри каждого кванта ∆i будут одинаковые, то, подставляя функцию плотности равномерного распределения в
46
формулу (3.6), получим средний квадрат ошибки вида
ε2X
N
∫ti
∑
=
i=1 t(i−1)
(x −
yi )2 N1∆ dx
=
1
N∆
N
∑
yi + ∆
2
∫
(x − yi )2 dx =
i=1 yi − ∆
2
=
1
∆
yi + ∆
2
∫
(3.8)
(x − yi ) dx.
2
yi − ∆
2
Ошибка внутри кванта с номером i имеет вид функции e(x) = x−yi , т.е. прямой пересекающей ось абсцисс в точке, равной аппроксимирующему значению yi этого кванта, как показано на рис.3.4б. Максимальное по модулю значение ошибки не превышает значения ∆
2 . Выполним подстановку x−yi = z,
в результате которой выражение (3.8) примет вид
∆
ε2X
∫2
1
=
∆
z 2 dz.
(3.9)
−∆
2
Таким образом, вычислив интеграл (3.8), получим
∆2
.
(3.10)
12
Подставляя формулы (3.7) и (3.10) в (3.5), получим окончательную функциональную зависимость отношения «сигнал-шум» от логарифма (по основанию два) числа квантов R, которая имеет линейный характер:
(
)
SN R(R) = 10 lg (2R )2 = 20R lg 2 ≃ 6.02R dB.
(3.11)
ε2X =
Для случайных величин, распределенных не по равномерному закону,
значение среднеквадратической ошибки будет отличаться от приведенной в
формуле (3.10). Покажем, что при большом числе квантов для произвольного
распределения ошибку равномерного квантование также можно оценивать
с помощью формулы (3.11). Для этого перепишем формулу (3.6) в форме
условного математического ожидания:
ε2X =
N
∑
[
] {
}
M e2 (x) | x ∈ [t(i−1) ; ti ) Pr x ∈ [t(i−1) ; ti ) .
(3.12)
i=1
При этом
N
∑
{
}
Pr x ∈ [t(i−1) ; ti ) = 1
(3.13)
i=1
47
и условное математическое ожидание квадрата ошибки квантования случайной величины X при условии, что ее значение попадает в i-й квант, определяемый на интервале [t(i−1) ; ti ), равно:
∫ti
[
]
M e2 (x) | x ∈ [t(i−1) ; ti ) =
( {
})
(x − yi )2 f x | x ∈ [t(i−1) ; ti ) dx.
t(i−1)
(3.14)
( {
})
Условная плотность вероятности f x | x ∈ [t(i−1) ; ti ) , представленная в
формуле (3.14), определяется отношением плотности вероятности f (x) к вероятности попадания непрерывной случайной величины X в этот интервал.
Таким образом
[
]
M e2 (x) | x ∈ [t(i−1) ; ti ) =
∫ti
(x − yi )2
t(i−1)
f (x)
∫ti
dx.
(3.15)
f (y)dy
t(i−1)
Для равномерного распределения f (x) = N1∆ при любом допустимом
∫ti
значении x, а
f (y)dy = N1 . Подставив соответствующие значения в
t(i−1)
(3.15) и (3.12), придем к результату, полученному в формуле (3.8).
Для распределения, отличного от равномерного, значение f (x) на любом
из интервалов не будет постоянным. Для того, чтобы оценить ε2X для произвольного распределения, построим на интервале [t(i−1) ; ti ) прямоугольник
со сторонами ∆ и h (рис. 3.5), таким образом, чтобы выполнялось условие
∫ti
f (y)dy = hi ∆.
t(i−1)
При этом на интервале [t(i−1) ; ti ) для распределения, отличного от равномерного, f (x) может не совпадать с hi . Однако при большом числе квантов
интервал ∆ становится малым и можно считать, что значение f (x) на всем
интервале постоянно, т.е. f (x) ≈ hi . В итоге выражение (3.12) сводится к
[
]
N
{
}
∫ti
∑
2
2 1
εX ≈
(x − yi ) ∆ dx Pr x ∈ [t(i−1) ; ti ) =
i=1 t(i−1)
(3.16)
N
{
} ∆2
∑
∆2
= 12
Pr x ∈ [t(i−1) ; ti ) = 12 .
i=1
48
Рис. 3.5. Приближенное вычисление вероятности попадания в квант [t(i−1) ; ti ) при равномерном скалярном квантовании
Таким образом из совпадения (3.10) и (3.16) следует, что формулу (3.11)
также можно использовать для приближенной оценки ошибки равномерного
квантования значений случайной величины с произвольным распределением. Точность этой оценки необходимо будет проанализировать в задании 14
в ходе выполнения данного исследования.
1.3.3. Неравномерное скалярное квантование. Алгоритм Макса-Ллойда
Равномерное квантование не является эффективным для случая, когда
плотность распределения случайной величины отлична от равномерного закона. При минимизации среднеквадратической ошибки ε2X размеры квантов
уже не будут одинаковыми и оптимальные аппроксимирующие значения не
обязательно будут располагаться в середине кванта. На рис. 3.6 этот факт поясняется на качественном уровне. Для минимизации ошибок размеры квантов в области наиболее частых значений должны быть меньше, т.е. скалярное
квантование становится неравномерным.
Для формирования параметров оптимального неравномерного квантования сформулируем следующую оптимизационную задачу. Пусть задана
функция плотности вероятности f (x) непрерывной случайной величины X
и количество квантов N . Требуется выбрать значения границ квантов {ti }
и аппроксимирующих значений {yi }, минимизирующие целевую функцию,
49
0.14
yi
0.12
ti
f(x)
0.1
0.08
0.06
0.04
0.02
0
-10
-8
-6
-4
-2
0
2
4
6
8
10
x
Рис. 3.6. Пример неравномерного скалярного квантования
которой является средний квадрат ошибки квантования (3.6). Результатом решения задачи минимизации целевой функции (3.6) является одновременное
выполнение двух следующих условий:
ti =
y(i+1) + yi
;
2
∫ti
yi =
xf (x) dx
t(i−1)
∫ti
(3.17)
.
(3.18)
f (x) dx
t(i−1)
На практике функция плотности вероятности f (x) неизвестна, и для
формирования параметров оптимального неравномерного квантования используют выборку случайной величины X. Итеративный алгоритм, который
на основе выборки находит приближенное значение {ti } и {yi } в соответствии с формулами (3.17) и (3.18) был сформулирован Ллойдом в 1956 г. и
опубликован Максом [11] в 1960 г. (работа Ллойда была опубликована в 1982
[12]). Алгоритм состоит из следующих шагов.
Входные данные:
– выборка {x1 , x2 , ..., xM };
50
– число квантов N , причем M ≫ N .
Вспомогательные переменные:
– k - номер текущей итерации алгоритма;
(k)
– ti
–
–
–– пороговые значения квантователя на итерации k;
(k)
yi –– аппроксимирующие значения квантователя на итерации k;
(k)
|Qi (x)| –– количество элементов из исходной выборки {x}, которые
по-
падают в i-й квант на итерации k;
– Tε2 –– предустановленный порог для ошибки квантования;
– ∆ε2 –– предустановленный порог для разности ошибок квантования на текущей (k) и предыдущей (k − 1) итерациях;
– Tk –– предустановленный порог для максимального числа итераций.
Инициализация:
– k = 1.
(1)
– Определение начальных значений ti . Начальные значения можно задать
произвольным образом, но так, чтобы выполнялось два условия:
(1)
(1)
t(i−1) < ti ;
(1)
|Qi (x)| ̸= 0,
(3.19)
∀i = 1, · · · , N.
(3.20)
Последнее условие означает, что в каждый квант должно попасть как минимум одно значение из исходной выборки {xj }. Наиболее простой спо(1)
соб –– использовать для инициализации {ti } значения границ равномерного квантователя:
(1)
t0 = min{xj };
j
(1)
tN
(1)
ti
(1)
= max{xj };
= t0 + i
j
(1)
tN
(1)
− t0
.
N
После этого надо убедиться, что условие (3.20) выполняется . В противном
(1)
случае необходимо изменить ранее сформированное множество {ti }.
51
(k)
(k)
Шаг 1. Для каждого кванта i, определяемого интервалом [t(i−1) ; ti )
вычисляется среднее значение среди элементов выборки, попавших в этот
интервал. Вычисленное среднее присваивается аппроксимирующему значе(k)
нию yi :
∑
1
(k)
yi = (k)
xj .
(3.21)
[
)
|Qi (x)|
(k)
(k)
xj ∈ t(i−1) ;ti
При работе с выборкой формула (3.21) соответствует формуле (3.18).
Шаг 2. Вычисление ошибки квантования ε2X (k) на текущей итерации k
(k)
в соответствии с текущими аппроксимирующими значениями yi и грани(k)
цами ti :
ε2X (k) =
M
1 ∑
2
(xj − Q(xj )) ,
M j=1
(k)
Q(xj ) = ys(k) : xj ∈ [ts−1 , t(k)
s ).
(3.22)
(3.23)
При работе с выборкой формула (3.22) соответствует (3.3). Напомним,
что Q(xj ) определяет правило выбора аппроксимирующего значения по
кванту, в который попадает xj .
Шаг 3. Анализ условий останова алгоритма:
2
1. εX (k) < Tε2 ;
2. ε2X (k) − ε2X (k − 1) < ∆ε2 ;
3. k == Tk .
Если одно из вышеперечисленных условий выполняется, то алгоритм завершает свою работу.
(k)
Шаг 4. Установка значений границ ti в соответствии с формулой (3.17)
Шаг 5. k = k + 1. Переход к шагу 1.
Строго говоря, алгоритм никогда не достигает оптимального значения,
так как одно из условий (3.17) или (3.18) остается невыполненным после
остановки. В оригинальном описании алгоритма [12] проверка на выполнение условия останова (шаги 2 и 3) выполняется после каждой модификации:
(k)
(k)
аппроксимирующих значений yi (шаг 1) и граничных значений ti (шаг
4). Изложенный вариант с одной проверкой останова обусловлен практическими соображениями ускорения процедуры. К тому же такой вариант позволяет не хранить границы квантов {ti }, так как номер кванта может быть
определен по ближайшему к x аппроксимирующему значению yi .
52
1.4. Векторное квантование. Алгоритм Линде-Бузо-Грея
Алгоритм построения оптимального векторного квантователя является
обобщением скалярного случая. Поэтому его иногда еще называют обобщенным алгоритмом Макса-Ллойда. Аппроксимирующие векторы ȳ и квантуемые векторы x̄ можно интерпретировать, как точки в n-мерном пространстве, где n –– длина векторов. Чтобы определить схему векторного квантования, необходимо сформировать множество из N аппроксимирующих векторов {ȳ}, которое называют кодовой книгой. Процедура квантования ȳ = Q(x̄)
заключается в выборе для x̄ наиболее близкого по заранее определенной
метрике вектора ȳ из кодовой книги. Границы между квантами определены
условно, исходя из выбранной метрики. В качестве метрики будем вновь использовать евклидово расстояние. Геометрической интерпретацией векторного квантования может служить диаграмма Вороного, пример которой приведен на рис. 3.7.
Рис. 3.7. Пример диаграммы Вороного
Пример 4. Палитра BMP файла как кодовая книга.
Рассмотрим в качестве примера кодовой книги палитру в формате изображений BMP8. Для
представления цвета пикселя также используется стандартная тройка базовых цветов (R, G, B),
однако тройки выбираются не произвольно, как в формате BMP24, а из набора 256 вариантов,
определенных в палитре, которая располагается между заголовком и данными. Поскольку выбор осуществляется из 28 вариантов, для представления пикселя необходим один байт. Более
подробное описание работы с палитрой в формате BMP8 представлено в подразделе 1.5
53
Обобщение метода Ллойда-Макса в англоязычной литературе также известно, как алгоритм Линде-Бузо-Грея [13], [14]. Алгоритм состоит из следующих шагов.
Входные данные:
– выборка векторов {x̄1 , x̄2 , ..., x̄M }; векторы имеют длину n;
– размер кодовой книги N , причем M ≫ N .
Вспомогательные переменные:
– k - номер текущей итерации алгоритма;
(k)
– ȳi , i = 1, · · · , N –– аппроксимирующие векторы кодовой книги на
итерации k;
– Si , i = 1, · · · , N –– набор счетчиков, фиксирующих количество векто(k)
ров исходной выборки, для которых вектор ȳi будет являться аппрок(k)
симирующим на текущей итерации. По смыслу Si аналогично |Qi (x)| в
алгоритме Ллойда-Макса.
– v̄i , i = 1, · · · , N –– набор из векторов-счетчиков для накопления суммы
квантуемых векторов, отнесенных к i-му кванту;
– Tε2 –– предустановленный порог для ошибки квантования;
– ∆ε2 –– предустановленный порог для разности ошибок квантования на текущей (k) и предыдущей (k − 1) итерациях;
– Tk –– предустановленный порог для максимального числа итераций.
Инициализация:
– k = 1.
– Формируем кодовую книгу ȳ (1) из N векторов исходной выборки произвольным образом. Формирование кодовой книги на основе первых N
векторов исходной выборки не рекомендуется в силу высокой корреляции между соседними векторами в выборке. Это особенно характерно для
изображений.
Шаг 1. Вычисление ошибки квантования ε2X (k) на текущей итерации k
с использованием кодовой книги {ȳ (k) }:
ε2X (k) =
M
1 ∑
∥x̄j − Q(x̄j )∥2 ,
M j=1
(k)
Q(x̄j ) = ȳs(k) : s = argmin{∥x̄j − ȳi ∥2 }.
i
При работе с выборкой формула (3.24) соответствует (3.3).
Шаг 2. Анализ условий останова алгоритма:
54
(3.24)
(3.25)
1. ε2X (k) < Tε2 ;
2. ε2X (k) − ε2X (k − 1) < ∆ε2 ;
3. k == Tk .
Если одно из вышеперечисленных условий выполняется, то алгоритм завершает свою работу.
Шаг 3. Обнуление всех счетчиков Si и векторов-счетчиков v̄i :
v̄i = 0̄,
i = 1, · · · , N ;
Si = 0,
i = 1, · · · , N.
Шаг 4. Для каждого вектора выборки x̄j , j = 1, · · · , N выполняются
следующие действия.
– Определяется наиболее близкий по евклидову расстоянию аппроксими(k)
рующий вектор ȳi .
(k)
t = argmin{∥x¯j − ȳi ∥2 },
i
(k)
где ∥x¯j − ȳi ∥ вычисляется по формуле (3.2).
– Модификация вектора v̄t :
v̄t = v̄t + x̄j .
Сложение векторов выполняется по правилам линейной алгебры (поэлементно).
– Наращивание счетчика St :
St = St + 1.
Шаг 5. Формирование новой кодовой книги, которая будет использоваться на следующей итерации k + 1.
(k+1)
ȳi
=
v̄i
,
Si
i = 1, · · · , N.
Деление вектора v̄i на скаляр осуществляется по правилам линейной алгебры (поэлементно).
Шаг 6. k = k + 1. Переход к шагу 1.
55
Границы квантов можно формально получить по построенной кодовой
книге, определив множество векторов {t̄i,j }, которые равноудалены от аппроксимирующих векторов в ȳi и ȳj . Для этого воспользуемся следующим
равенством:
2
2
∥t̄i,j − ȳi ∥ = ∥t̄i,j − ȳj ∥ .
Совокупность векторов {t̄i,j } для всех возможных значений пар i и j формирует так называемые многогранники Вороного.
1.5. Использование палитры в формате BMP
С помощью формата BMP24 можно представить до 224 ∼ 16 млн. различных цветовых оттенков одного пикселя. Однако далеко не все из них используются в одном изображении. Для сокращения размера файла ранее в
формате BMP активно использовался способ ссылочного представления значения пикселя вместо непосредственного, как для 24-битной версии. Для
организации такого представления используется специальная таблица, которую принято называть палитрой (palette).
В файле BMP палитра хранится после заголовочной структуры
BITMAPINFOHEADER. Палитра представляет собой логическую таблицу,
строки которой последовательно записаны в файле с помощью структур типа RGBQUAD. Структура RGBQUAD объявляется следующим образом:
typedef struct tagRGBQUAD {
BYTE rgbBlue;
BYTE rgbGreen;
BYTE rgbRed;
BYTE rgbReserved;
} RGBQUAD;
В первых трех полях хранится комбинация значений синей (rgbBlue), зеленой ( rgbGreen) и красной (rgbRed) компонент. Значение поля rgbReserved не
используется. Таким образом, строка таблицы описывает цветовой оттенок,
используя 24 бита. Значение пикселя в разделе данных интерпретируется как
ссылка на строку в таблице палитры.
Количество экземпляров структуры RGBQUAD, представленных
в разделе палитры, определяется в переменной biClrUsed структуры
BITMAPINFOHEADER. Очевидно, что при затрате k бит на представление
одного пикселя, размер палитры следует выбирать 2k для обеспечения
56
максимально возможного числа используемых цветовых оттенков в этом
файле. Чтобы обеспечить такой способ представления данных, следует
установить значение поля biClrUsed в ноль. В этом случае размер палитры
будет определяться полем biBitCount структуры BITMAPINFOHEADER,
которое в случае использования палитры может принимать значения 1, 2, 4
или 8.
2. Порядок самостоятельного выполнения исследования
В ходе выполнения исследования необходимо выполнить следующие задания.
1. Проверка точности оценки SNR, получаемой по формуле (3.11):
1.1 сформировать выборку объемом M = 10000 равномерно распределенной случайной величины X ∼ U [−a, a]. Значение параметра a
определяет преподаватель.
1.2 применить к выборке {x1 , ..., xM } равномерное квантование с числом квантов 2R , где R = 1, ..., 7.
1.3 По результатам эксперимента построить график SNR(R), воспользовавшись формулой (3.4) и сравнить его с теоретическим, полученным по формуле (3.11). Для вычисления SNR по теоретическим
значениям, используя дисперсию исследуемого закона распределения, следует пользоваться формулой (3.5).
1.4 выполнить аналогичные действия для выборки случайной величины, распределенной по нормальному закону N (0, σ 2 ), где σ = a3 .
Примечания по заданию:
– интервал квантования установить [−a; a] (как и для случайной
величины X),
– t0 = −∞, tN −1 = +∞ (N –– число квантов).
1.5 провести сравнение построенных теоретического и экспериментальных графиков, сделать выводы о возможности использования
теоретического расчета для экспериментальных выборок данных;
2. Анализ скалярного квантования:
2.1 получить у преподавателя три изображения в формате RGB24, различающиеся по яркости:
2.1.1. с преобладанием светлых тонов («засвеченное» изображение),
57
2.1.2. сбалансированное по яркости изображение,
2.1.3. затемненное изображение.
2.2 для каждого получить яркостную составляющую Y , воспользовавшись формулой [8], [9]:
Y
=
0.299R + 0.587G + 0.114B;
(3.26)
все дальнейшие действия проводить только с Y компонентой;
2.3 для каждого тестового изображения построить гистограмму [15] яркостной компоненты Y ;
2.4 выполнить равномерное скалярное квантование исходных данных
на 2R уровней, где R = 1, ..., 7 для выборки, сформированной на
основе яркостной Y . Построить график SNR(R) по результатам выполненного квантования.
2.5 для каждого значения R вычислить PSNR, воспользовавшись формулой [8], [9]:
P SN R = 10 lg
W H(2L − 1)2
,
H ∑
W
∑
2
(Yi,j − Ŷi,j )
(3.27)
i=1 j=1
где Ŷi,j –– аппроксимирующее значение, вычисленное для яркости
пикселя на позиции (i, j), а также оценить степень сжатия данных
после квантования, воспользовавшись формулой оценки энтропии
[10]:
∑
Ĥ(X) = −
p̂(x) log2 p̂(x).
(3.28)
x
2.6 по вычисленным значениям для равномерного квантования построить графики PSNR(R) и PSNR(H);
2.7 реализовать алгоритм Макса-Ллойда для оптимального неравномерного скалярного квантования на 2R уровней. Применить алгоритм к каждому изображению;
2.8 выполнить неравномерное скалярное квантование исходных данных на 2R уровней, где R = 1, ..., 7;
2.9 для каждого значения R вычислить PSNR, воспользовавшись формулой (3.27) и оценить энтропию H по формуле (3.28);
58
2.10 по вычисленным значениям для неравномерного квантования построить графики PSNR(R) и PSNR(H).
2.11 построить график ступенчатой функции Q(x) по параметрам оптимального неравномерного квантования для каждого изображения.
Число квантов определяет преподаватель;
2.12 провести сравнительный анализ графиков, полученных с помощью
неравномерного и равномерного скалярного квантования.
3. Создание палитры изображения с помощью методов векторного квантования:
3.1 используя метод Линде-Бузо-Грея, составить палитру по данным
компонент R, G и B исходного изображения в формате BMP24;
3.2 сформировать выходное изображение в формате BMP8 с палитрой;
3.3 оценить ошибку восстановления входного изображения. В качестве
критерия эффективности необходимо использовать PSNR;
3.4 проанализировать влияние округления на результат формирования
палитры (точность аппроксимации, число итераций). На шаге 5 алгоритма Линде-Бузо-Грея значения векторов кодовой книги, соответствующей палитре, перестают быть целочисленными. Округление элементов палитры можно осуществлять в процессе шага 5 (что
приведет к увеличению ошибки квантования) или по окончании работы алгоритма. В задании требуется проанализировать насколько
возрастет ошибка квантования для первого варианта.
59
IV. АНАЛОГОВЫЕ ЛИНЕЙНЫЕ
СИСТЕМЫ
Цель исследования: изучение способов описания аналоговых линейных
систем, а также методов построения реализуемых аналоговых фильтров.
1. Теоретические основы
1.1. Основные определения
В данном разделе рассматриваются линейные стационарные системы,
которые осуществляют преобразование детерминированных сигналов. Сигнал uin (t), который подвергается преобразованию, называют входным. Сигнал uout (t), являющийся результатом преобразования называют реакцией системы или выходным сигналом. Процедуру преобразования сигнала будем
обозначать в операторной форме
uout (t) = S {uin (t)} .
Систему называют линейной, если для нее выполняется правило суперпозиции –– реакция на линейную комбинацию N входных сигналов является линейной комбинацией N реакций на отдельные входные сигналы, причем весовые коэффициенты a(i) линейных комбинаций на входе и на выходе системы совпадают:
{N
}
N
∑
∑
(i) (i)
(i) (i)
a uout (t) = S
a uin (t) ,
(4.1)
i=1
где
(i)
uout (t)
i=1
{
}
(i)
= S uin (t) . Равенство (4.1) должно выполняться для любых
коэффициентов a(i) .
Система называется стационарной, если реакция системы на один и тот
же входной сигнал не зависит от момента времени, когда этот сигнал был
подан на вход системы. Таким образом, если входной сигнал и его смещенная
на произвольное время τ копия совпадают, то:
S {u(t ± τ )} = S {u(t)} .
60
1.2. Основные характеристики линейных систем
1.2.1. Импульсная характеристика системы
Основной характеристикой, используемой для анализа линейных стационарных систем является импульсная характеристика [3] h(t) –– реакция системы на входной сигнал в виде дельта-функции (для дискретных систем ––
единичный импульс):
h(t) = S {δ(t)} .
В соответствии с фильтрующим свойством дельта-функции:
∫∞
uin (τ )δ(t − τ )dτ.
uin (t) =
(4.2)
−∞
При прохождении через линейную систему, входной сигнал uin (t) преобразуется в выходной сигнал, а дельта-функция δ(t − τ ) в импульсную характеристику h(t − τ ). Функция uin (τ ) останется без изменений по скольку не
зависит от t. Тогда тождество (4.2) может быть представлено в следующем
виде [16], [3]:
∫∞
uout (t) =
uin (τ )h(t − τ )dτ.
(4.3)
−∞
Выражение (4.3) в [16] называют интегралом Дюамеля.
Нужно отметить, что для физической реализуемости подобной системы,
она должна обладать свойством причинности - то есть на выходе системы отличный от нуля отклик возникает только после того как был подан сигнал на
ее вход. Данное свойство выражается в том, что импульсная характеристика
реализуемой системы должна быть равна нулю при t < 0.
Так же для линейных систем определено свойство устойчивости. Система называется устойчивой, если при отсутствующем (нулевом) входном
сигнале выходной сигнал затухает (стремится к нулю) вне зависимости от
начальный условий:
lim uout (t) = 0 : uin (t) = 0.
t→∞
(4.4)
Данное свойство системы может быть выражено и через импульсную
характеристику:
lim h(t) = 0.
(4.5)
t→∞
61
Необходимое условие для устойчивости системы будет позже описано в подразделе 1.3.2.
1.2.2. Комплексный коэффициент передачи
Согласно выражению (4.3), выход линейной системы представляет собой свертку входного сигнала и импульсной характеристики системы [16].
Выполнив преобразование Фурье для данного равенства (4.3), получим:
U̇out (ω) = U̇in (ω)K̇(ω).
(4.6)
Фурье образ K̇(ω) импульсной характеристики называют комплексным коэффициентом передачи. Значение данного комплексного коэффициента
K̇(ω) показывает, как изменится при прохождении через линейную систему синусоида, заданная набором параметров (A, ω, ϕ):
uo ut(t) = A|K̇(ω)| sin(ωt + ϕ + ∠K̇(ω)).
(4.7)
Из формулы (4.7) видно, что модуль |K̇(ω)|, называемый амплитудночастотной характеристикой (АЧХ), показывает изменение амплитуды
входного сигнала при прохождении через данную систему. А аргумент
∠K̇(ω), называемый фазо-частотной характеристикой (ФЧХ), показывает изменение фазы входного сигнала на соответствующей частоте ω на угол
∠K̇(ω).
1.3. Способы описания линейных систем
1.3.1. Дифференциальное уравнение
Дифференциальное уравнение является наиболее общим и универсальным подходом для анализа взаимосвязи входного и выходного сигналов линейной системы, что позволяет исследовать ее характеристики. Все последующие изложенные методы так или иначе опираются на него.
n
∑
∑ dj
di
ai i y + a0 y(t) =
bj j x + b0 x(t).
dt
dt
i=1
j=1
m
(4.8)
Для выполнения условия физической реализуемости необходимо, чтобы
m ⩽ n. Для случая линейной системы коэффициенты ai и bj являются вещественными константами. Число n принято называть порядком системы.
62
1.3.2. Операторный метод. Функция передачи
Как было показано ранее (4.6), для описания линейной системы можно использовать спектральные методы. В частности рассмотренное ранее
преобразование Фурье, представляющее исследуемый сигнал в виде суммы
неограниченно большого числа элементарных слагаемых, каждое из которых периодически изменяется во времени по закону exp(jωt).
Обобщением данного представления будет использование в качестве
элементарных слагаемых комплексных экспоненциальных сигналов вида
exp(s), где s - комплексное число (s = σ + jω), называемое комплексной
частотой [16]. Вещественный сигнал может быть описан суммой пары экспоненциальных сигналов с комплексной частотой:
∗
est + es
s(t) =
2
t
=
ejω + e−jω
eσ+jω + eσ−jω
= eσt
= eσt cos(ωt). (4.9)
2
2
Из данного представления следует, что форма получаемого вещественного
сигнала зависит от выбора вещественной и мнимой части комплексной частоты:
– при σ = 0 и ω ̸= 0 получаются гармонические колебания вида cos(ωt),
изображенное на рис. 4.1в;
– при σ ̸= 0 и ω = 0 получаются возрастающие (рис. 4.1а) или убывающие
(рис. 4.1б) во времени экспоненты eσt ;
– при σ ̸= 0 и ω ̸= 0 получаемый вещественный сигнал имеет сложную
форму, в виде возрастающего или затухающего гармонического колебания
с частотой ω и огибающей eσt (пример затухающего колебания на рис.
4.1г).
Введенное понятие комплексной частоты часто используется для описания разнообразных линейных систем, поскольку экспоненциальные сигналы
вида (4.9) схожи с откликами таких систем. Преобразование при использовании комплексной частоты называется преобразованием Лапласа и осуществляется согласно следующей формуле:
∫∞
U (s) = L {u(t)} =
u(t)e−st dt.
(4.10)
0
Также как и преобразование Фурье, преобразование Лапласа является
линейным, поэтому для него справедливы все свойства последнего. Изобра63
а)
б)
в)
г)
Рис. 4.1. Формы получаемого вещественного сигнала при использовании комплексной частоты различных видов: а–– при σ > 0 и ω = 0; б–– при σ < 0 и ω = 0; в–– при σ = 0 и ω ̸= 0;
г–– при σ < 0 и ω ̸= 0
жение по Лапласу для производной функции f (x) k-го порядка при нулевых
начальных условиях выглядит следующим образом:
{ (k)
}
d
L
f (x) = sk F (s).
(4.11)
dtk
Применив оператор Лапласа к выражению (4.8), получим:


( n
)
m
∑
∑
ai si + a0 Y (s) = 
bj sj + b0  X(s).
i=1
(4.12)
j=1
Отсюда получаем передаточную функцию линейной системы:
m
∑
Y (s)
j=1
H(s) =
= ∑
n
X(s)
i=1
64
bj sj + b0
.
ai si + a0
(4.13)
Чтобы перейти от передаточной функции к комплексному коэффициенту
передачи, достаточно заменить аргумент первой на jω:
m
∑
bj (jw)j + b0
j=1
K̇(w) = H(jw) = ∑
n
.
(4.14)
ai (jw)i + a0
i=1
Передаточная функция линейной системы представляется в виде отношения двух полиномов. По скольку каждый полином может быть разложен
на множители, передаточная функция может быть представлена в виде отношения произведения линейных множителей многочлена числителя и произведения линейных множителей многочлена знаменателя:
∏m
j=1 (s − zj )
H(s) = k ∏n
,
(4.15)
i=1 (s − pi )
где k = bm /an называют коэффициентом усиления. zj - множество нулей
функции передачи, а pi - множество полюсов функции передачи. В районе
точек нулей АЧХ системы будет иметь относительно небольшие значения,
при этом в самих нулях |H(zi )| = 0. В районе точек полюсов АЧХ системы
будет иметь относительно большие значения и в самих полюсах |H(pi )| →
∞.
Передаточная функция линейной системы (4.14) в случае отсутствия
кратных полюсов может быть представлена как сумма простых дробей:
H(s) =
n−1
∑
i=0
Ci
+ k,
(s − pi )
(4.16)
где pi - i-ый полюс; Ci - i-ый вычет, k - целая часть функции передачи, отличная от нуля только в случае равенства степеней полиномов нулей и полюсов.
В случае наличия кратных полюсов, каждый l-кратный полюс pi дает l
слагаемых следующего вида:
l
∑
j=1
Ci,j
.
(s − pi )j
(4.17)
Форма записи передаточной функции в виде вычетов и полюсов позволяет получить формулу для импульсной характеристики системы путем
65
выполнения обратного преобразования Лапласа для каждого слагаемого. Результат обратного преобразования Лапласа для l-кратного полюса pi будет
представлен следующей формулой [17]:


l
l
∑
Ci,j  ∑ Ci,j tj−1 pi t
L−1
=
e .
(4.18)

(s − pi )j  j=1 (j − 1)!
j=1
Введенное определение импульсного отклика системы через функцию
от вычетов и полюсов позволяет определить множество допустимых значений для полюсов и вычетов, при которых система будет считаться устойчивой, т.е. для нее будет выполняться равенство lim h(t) = 0. Исходя из форt→∞
мулы (4.18), каждое слагаемое будет стремится к нулю если вещественная
часть соответствующего полюса будет отрицательной:
Re(pi ) < 0.
(4.19)
Отсюда можно сформулировать общее условие устойчивости для линейных систем [3]: линейная система является устойчивой тогда и только тогда, когда все полюсы ее функции передачи лежат в левой комплексной полуплоскости.
1.3.3. Фазовая и групповая задержки
Фазовая задержка как функция от круговой частоты ω рассчитывается
по следующей формуле:
τ (ω) = −
∠K̇(ω)
с.
ω
(4.20)
Групповая задержка вычисляется как производная ФЧХ:
τ (ω) = −
d∠K̇(ω)
.
dω
(4.21)
Напомним, что ФЧХ ∠K̇(ω) определяет фазовый сдвиг, который вносит система, описываемая комплексным коэффициентом передачи K̇(ω), в
спектр входного сигнала на частоте ω. При этом сдвиг фазы будет соответствовать временному сдвигу, рассчитываемому по формуле фазовой задержке (4.20). Причем на разных частотах фазовый сдвиг может оказаться раз66
ным. Насколько это окажется критичным для формы выходного сигнала можно определить по групповой задержке. Если ФЧХ системы линейна, то групповая задержка является константой и форма выходного сигнала будет повторять форму входного с некоторой задержкой (определяемой по все тому
же фазовому сдвигу). В случае нелинейной ФЧХ, разные частоты приобретают разные временные сдвиги, что будет приводить к искажению амплитуд
выходного сигнала во временной области. Поэтому при проектировании систем, для которых важно сохранить форму сигнала, необходимо учитывать,
что ФЧХ проектируемой системы должна быть линейной.
1.4. Аналоговые фильтры
1.4.1. Идеализированные фильтры
Основное предназначение фильтров - выделение или подавление определенных частот входного сигнала. Аналоговые линейные фильтры, как и
любые другие линейные системы, определяются с помощью передаточной
функции H(s).
Рассмотрим в качестве примера такого аналогового фильтра следующий
фильтр, пропускающий без искажений все частоты ниже частоты ωp , равной
1 рад/с и называемой частотой среза. Все частоты выше ωp полностью подавляются. В соответствии с определением данного фильтра его АЧХ будет
иметь следующий вид:
{
1, w < 1,
(4.22)
|K̇(w)| =
0, w > 1.
Приведенный в качестве примера фильтр (4.22) называется фильтромпрототипом, по скольку используется для построения фильтров других видов, и относится к виду фильтров нижних частот (ФНЧ). Остальные ФНЧ
отличаются от фильтра-прототипа тем, что их частота среза ωp произвольно задана, а их АЧХ (представлена на рис. 4.2) описывается формулой:
{
1, w < ωp
(4.23)
|K̇(w)| =
0, w > ωp .
Фильтры нижних частот, отличные от фильтра-прототипа, получаются с
помощью применения к последнему преобразования частоты [3]. Преобра67
Рис. 4.2. Амплитудно-частотная характеристика (АЧХ) фильтра нижних частот (ФНЧ)
зование частоты комплексного коэффициента передачи производит его масштабирование по оси аргумента следующим образом:
K̇new (ω) = K̇old (fω (ω)) ,
(4.24)
где под K̇old (ω) понимается комплексный коэффициент передачи исходного
фильтра, K̇new (ω) - комплексный коэффициент передачи преобразованного
фильтра, fω (ω) - функция преобразования частоты.
Данная процедура может быть выражена и через замену переменнойаргумента s в функции передачи H(s) исходного фильтра на некоторую подстановочную функцию fs (s):
Hnew (s) = Hold (fs (s)) .
(4.25)
Зная, что комплексный коэффициент передачи функции выражается из ее
функции передачи путем замены s = jω, можно выразить функцию преобразования частоты через подстановочную функцию следующим образом:
fω (ω) =
fs (jω)
.
j
(4.26)
При этом для выполнения равенств (4.25) и (4.26) подстановочная функция
передачи fs (s) должна удовлетворять двум основным требованиям:
1. fs (s) должна быть дробно-рациональной;
2. fs (s) при чисто мнимом результате должна давать чисто мнимый результат.
Для получения ФНЧ с произвольно заданным значением ωp из фильтрапрототипа с частотой среза 1 рад/c используется подстановочная функция
следующего вида:
s
fs (s) =
.
(4.27)
ωp
68
Функция преобразования частоты при такой подстановочной функции
fs (s) задается формулой:
ω
fω (ω) =
.
(4.28)
ωp
Фильтрами высоких частот (ФВЧ) называют такие фильтры, которые
пропускают без искажения все частоты выше заданной частоты среза ωp и
полностью подавляющие все остальные частоты.
Рис. 4.3. Амплитудно-частотная характеристика (АЧХ) фильтра высоких частот (ФВЧ)
Исходя из сравнения графиков АЧХ для ФВЧ (рис. 4.3) и ФНЧ (рис. 4.2)
фильтров видно, что построение ФВЧ фильтра с заданными параметрами
на основе ФНЧ фильтра-прототипа потребует инверсии частотной оси. Для
выполнения подобной процедуры можно воспользоваться следующей подстановочной функцией:
ωp
fs (s) =
.
(4.29)
s
Тогда функция преобразования частоты будет иметь следующий вид:
ωp
fω (ω) = − .
(4.30)
ω
Полосовыми фильтрами (ПФ) (рис. 4.4) называют такие фильтры, которые пропускают без искажений только частоты в заданном диапазоне
[ω1 ; ω2 ]. Остальные частоты, не входящие в данный диапазон, полностью подавляются. Данный тип фильтров может быть задан не только через указание
диапазона пропускаемых частот, но и через пару значений: средняя частота
(ωcp = (ω1 + ω2 )/2) и ширина полосы пропускания (△ω = ω2 − ω1 ).
Фильтры данного вида так же могут быть получены из ФНЧ фильтрапрототипа. Для этого нужно использовать следующую подстановочную
функцию:
s2 + ωc2
fs (s) =
,
(4.31)
△ωs
69
Рис. 4.4. Амплитудно-частотная характеристика (АЧХ) полосового фильтра (ПФ)
где ωc - частота, равная среднему геометрическому краев полосы пропускания:
√
ωc = ω1 ω2 .
(4.32)
Функция преобразования частоты будет иметь соответствующий вид:
fω (ω) =
ω 2 − ωc2
.
(△ω) ω
(4.33)
Режекторные фильтры (РФ) (рис. 4.5) (иногда именуются как заграждающие фильтры или полосно-задерживающие) полностью пропускают на
выход все частоты, кроме лежащих в заданном диапазоне [ω1 ; ω2 ]. Частоты,
лежащие в заданном диапазоне, полностью подавляются. Так же как и полосовой фильтр, РФ может быть задан через среднюю частоту (ωcp ) и ширину
полосы пропускания (△ω).
Рис. 4.5. Амплитудно-частотная характеристика (АЧХ) режекторного фильтра (РФ)
Поскольку РФ фильтр может быть получен путем инверсии ПФ, для получения РФ можно инвертировать подстановочную функцию ПФ (4.31):
fs (s) =
70
△ωs
.
+ ωc2
s2
(4.34)
Соответствующая функция преобразования частот :
fω (ω) =
(△ω) ω
.
ω 2 − ωc2
(4.35)
1.4.2. Аппроксимации идеализированных фильтров
Приведенные в предыдущем подразделе формы фильтров являются идеальными (имеют АЧХ прямоугольной формы) и не могут быть физически реализованы. На практике в цифровой обработке сигналов используются различные виды приближения (аппроксимации) фильтров с АЧХ прямоугольной формы. Реализуемые на практике фильтры характеризуются помимо
частот срезов ωp , частотами заграждения ωs , которые вместе формируют несколько различных диапазонов частот (данные диапазоны частот представлены на примере АЧХ для ФНЧ рис. 4.6 ):
1. полоса пропускания - это диапазон частот, в котором вносимое фильтром
ослабление не превышает заданной величины;
2. полоса подавления - это диапазон частот, в котором вносимое фильтром
ослабление не меньше заданной величины;
3. переходная полоса - это диапазон частот, содержащий все частоты, не вошедшие в первые две категории. То есть в частотах данного диапазона
ослабление сигнала больше, чем в полосе пропускания, но меньше чем в
полосе задерживания.
Степень вносимого фильтром ослабления в полосе пропускания определяется параметром, называемым допустимый уровень пульсации Gp [18].
Данный параметр вычисляется по следующей формуле:
1
,
Gp = √
1 + ε2p
(4.36)
где ε2p некоторый заданный коэффициент.
Параметр Gs определяет минимальное необходимое подавление [18] в
полосе подавления и вычисляется по формуле:
1
Gs = √
,
1 + ε2s
(4.37)
где ε2s некоторый заданный коэффициент.
71
Рис. 4.6. Пример АЧХ реализуемого на практике ФНЧ
Данные параметры принято выражать в децибелах:
Rp = −20 lg(Gp ) = 10 lg(1 + ε2p )dB;
Rs = −20 lg(Gs ) = 10 lg(1 + ε2s )dB.
(4.38)
Функция аппроксимации квадрата АЧХ фильтра-прототипа может быть
представлена в следующем виде [18], [3]:
|K̇(ω)|2 =
1
,
1 + ε2p FN2 (ω)
(4.39)
где FN (ω) – некоторая аппроксимирующая функция; N – порядок аппроксимирующей функции.
Аппроксимирующая функция может быть выражена из уравнения (4.39)
и условия необходимости обеспечения заданного уровня подавления Gs при
частоте ω ≤ ωp [18]:
FN (ωs /ωp ) = εs /εp .
(4.40)
Порядок аппроксимирующей функции фильтра задает максимальное количество нулей и полюсов передаточной функции H(s). Увеличение поряд72
ка фильтра позволяет уменьшить переходную полосу, но при этом негативно влияет на количество элементов (емкостей и индуктивностей), требуемое
при реализации данного фильтра. При вычислении порядка фильтра нужно
помнить, что для обеспечения заданной ширины переходной полосы результат вычисления следует округлять в большую сторону.
Таким образом, для аппроксимации необходимо задать аппроксимирующую функцию фильтра FN (ω) и ее порядок N . Ниже приведен список основных способов аппроксимации.
1. Аппроксимация по Баттерворту, при которой аппроксимирующая функция равна:
FN (ω) = ω N .
(4.41)
2. Аппроксимация по Чебышеву:
FN (ω) = TN (ω),
(4.42)
где TN (ω) - многочлен Чебышева порядка N .
3. Аппроксимация по Чебышеву второго рода:
FN (ω) =
1
.
TN (1/ω)
(4.43)
4. Аппроксимация по Кауэру (эллиптическая аппроксимация):
FN (ω) = RN (ω),
(4.44)
где RN (ω) - эллиптическая дробно-рациональная функция порядка N .
1.4.3. Расчет фильтра с заданными параметрами на основе
аппроксимации по Баттерворту
Построение фильтра с заданными параметрами при известной функции
аппроксимации начинается с определения ее порядка, при котором получаемый фильтр удовлетворяет предъявляемым к нему требованиям: допустимому уровню пульсации Gp в полосе пропускания и минимальному необходимому подавлению Gs в полосе подавления. При заданных параметрах
{ωs ; ωp ; εs ; εp } порядок аппроксимирующей функции Баттерворта (4.41) может быть выражен следующим образом [18]:
( )
( )N
ωs
εs
ωs
εs
FN
=
−→
=
(4.45)
ωp
εp
ωp
εp
73
Если прологарифмировать формулу (4.45), получем выражение для определения порядка фильтра N :
⌈
⌉
log(εs /εp )
N=
(4.46)
log(ωs /ωp )
Вычисленный для заданных входных параметров порядок фильтра может не
быть целым числом. В таком случае полученное значение N следует округлить в большую сторону до ближайшего целого.
После того, как был найден порядок фильтра N , полюса фильтра Баттерворта (нули, как видно из формул (4.39) и (4.41), у данного фильтра отсутствуют) могут быть определены исходя из формулы (4.39) путем замены
ω = s/j = −js следующим образом:
|H(s)|2 =
1
.
1 + ε2s (−1)N s2N
(4.47)
Для нахождения полюсов приравняем знаменатель к нулю:
1 + ε2p (−1)N s2N = 0 → (−1)N s2N = −
1
.
ε2p
(4.48)
Следует отдельно рассмотреть случаи четного и нечетного порядка
фильтра N . При четных N , (-1) в правой части равенства можно заменить на
комплексную экспоненту exp (j(2n + 1)π) , n = 1, 2, · · · , 2N − 1, 2N , тогда:
s2N =
1
exp (j(2n + 1)π) , n = 1, 2, · · · , 2N − 1, 2N.
ε2p
(4.49)
Прологарифмируем полученное равенство для того, чтобы получить
формулу для нахождения полюсов фильтра четного порядка N (пример найденных полюсов для N = 4 изображен на рис. 4.7а):
(
)
1
2N
ln(s ) = ln 2 exp (j(2n + 1)π) ;
εp
2N ln(s) = −2 ln(εp ) + j(2n + 1)π;
1
2n + 1
ln(s) = − ln(εp ) + j
π;
N
2N )
(
2n + 1
1
s = exp − ln(εp ) + j
π ;
N
2N
(
)
1
2n + 1
s = √ exp j
π ; n = 1, 2, · · · , 2N − 1, 2N.
N ε
2N
p
74
(4.50)
При нечетных N , (1) в правой части равенства (4.48) заменем на следующую комплексную экспоненту exp (j2nπ) , n = 1, 2, · · · , 2N −1, 2N , тогда:
s2N =
1
exp (j2nπ) , n = 1, 2, · · · , 2N − 1, 2N.
ε2p
(4.51)
Прологарифмируем полученное равенство для того, чтобы получить
формулу для нахождения полюсов фильтра нечетного порядка N (пример
найденных полюсов для N = 5 изображен на рис. 4.7б):
(
)
1
ln(s2N ) = ln 2 exp (j2nπ) ;
εp
2N ln(s) = −2 ln(εp ) + j2nπ;
1
2n
ln(s) = − ln(εp ) + j
π;
N
2N
(
)
n
1
s = exp − ln(εp ) + j π ;
N
N
( n )
1
s = √ exp j π ; n = 1, 2, · · · , 2N − 1, 2N.
N ε
N
p
(4.52)
Ранее, в подразделе 1.3.2, говорилось о том, что для того чтобы линейная система считалась устойчивой, полюсы данной системы должны находится в левой полуплоскости комплексной плоскости. Следовательно, из получаемых по формулам (для четных и нечетных N соответственно) полюсов нужно выбрать только те, что лежат в левой полуплоскости. Примеры
устойчивых фильтров различных порядков представлены на рис. 4.7в-г. Для
нахождения полюсов устойчивого фильтра можно использовать следующую
формулу [18]:
1
π
pn = √ exp(j( + θn )),
N ε
2
p
2n − 1
θn =
π, n = 1, 2, · · · , N − 1, N.
2N
(4.53)
√
где 1/ N εp - радиус окружности, на которой расположены получаемые полюсы; θn - угол, на который повернут n-ый полюс на окружности.
Предполагается (рис. 4.2), что фильтр-прототип полностью пропускает
частоты из полосы пропускания, следовательно |H(0)| → 1. В случае, если
75
а)
б)
в)
г)
Рис. 4.7. Примеры полюсов фильтров Баттерворта разного порядка: а–– неустойчивый
фильтр 4 порядка; б–– неустойчивый фильтр 5 порядка; в–– устойчивый фильтр 4 порядка; г––
устойчивый фильтр 5 порядка
|H(0)| ̸= 1 , коэффициент передачи фильтра следует нормировать следующим образом:
H(s)
Hnorm (s) =
.
(4.54)
|H(0)|
2. Порядок самостоятельного выполнения исследования
В ходе выполнения исследования необходимо выполнить следующие задания.
76
1. Написать программу для построения графиков АЧХ и ФЧХ аналоговой
линейной системы, введенной пользователем. В качестве входных параметров программа должна получить массивы нулей и полюсов, их размеры, диапазон частот (минимальная, максимальная, шаг), для которых
рассчитывается коэффициент передачи.
2. Определить все параметры (по вариантам), необходимые для построения
фильтра Баттерворта: порядок фильтра N , частоты среза ωp и заграждения
ωs , параметров пульсации в полосе пропускания εp и в полосе подавления
εs , если заданы:
2.1 ωp = 1, ωs , εp , εs ;
2.2 N , ωs , εp , εs ;
2.3 N , ωp = 1, εp , εs ;
2.4 N , ωp = 1, ωs , εs ;
2.5 N , ωp = 1, ωs , εp .
3. Определить полюсы фильтра Баттерворта порядка N . Построить графики
найденных полюсов.
4. Построить графики АЧХ и ФЧХ исследуемого фильтра с определенными
параметрами.
5. Повторить пункты 3 - 4 для порядков фильтра (N − 1) и (N + 1).
Сделать выводы о влиянии порядка фильтра на его фильтрующие свойства.
6. Для фильтра порядка n сформировать ФВЧ с частотой среза ω0 . Построить
графики АЧХ, ФЧХ.
7. Сформировать полосовой фильтр с частотами (ω1 , ω2 ). Построить графики АЧХ, ФЧХ.
8. Сформировать режекторный фильтр с частотами (ω1 , ω2 ).Построить графики АЧХ, ФЧХ.
77
V. ДИСКРЕТНЫЕ ЛИНЕЙНЫЕ
СИСТЕМЫ
Цель исследования: изучение способов описания и синтеза дискретных
фильтров.
1. Теоретические основы исследования
1.1. Особенности преобразования Лапласа для дискретного сигнала
В данной работе, в отличие от предыдущей, будут исследованы особенности работы линейных стационарных систем с цифровыми и дискретными
сигналами. Дискретной системой называется система, импульсная характеристика которой является дискретной, а коэффициенты передаточной функции рассчитаны точно (без ошибок округления) [18]. Цифровой системой
будем считать дискретную систему, коэффициенты передаточной характеристики которой рассчитаны не точно, а с ошибками округления, вызванными конечной разрядностью представления числа. По скольку на практике
разрядность представления используемых чисел ограничена, все рассчитанны фильтры будут относиться к цифровым. При этом, в случае работы с 64битными числами с плавающей точкой, ошибка округления будет незначительно мала, что позволяет говорить о ”почти дискретной” системе.
Ранее было показано, что влияние системы на входной сигнал может
быть выражено в виде отношения образов Лапласа входного и выходного
сигналов. Поэтому для начала рассмотрим подробнее представление образа
Лапласа для дискретного сигнала. Предположим, что имеется дискретизированный сигнал xd (t), вида [18], [3]:
∑
xd (t) =
x(t)δ(t − nTd ), n = 0, 1, . . . ,
(5.1)
n
где Td –– период следования дельта-функций.
Тогда преобразование Лапласа от такого дискретизированного сигнала
78
xd (t) будет равно:
∫∞
Xd (s) =
xd (t)e
−st
dt =
0
=
∞
∑∫
n
∫∞ ∑
x(t)δ(t − nTd )e−st dt =
n
0
x(t)δ(t − nTd )e−st dt =
∑
x(nTd )e−snTd dt,
(5.2)
n
0
где s = σ + jω –– комплексная частота. Для получения данного результата
использовалось фильтрующее свойство дельта-функции.
Главное отличие преобразования, выполняемого для дискретизированного сигнала от случая обработки аналогового сигнала заключается в периодичности получаемого образа (при этом период равен Ω = 2π/Td ) вдоль
мнимой оси. Пример различия образов Лапласа для аналогового и дискретного сигналов приведен на рис. 5.1. Приведенные на рисунке значения взяты
только для демонстрации главного отличия образов Лапласа и не относятся
к образу какого либо реализуемого сигнала. Само свойство периодичности
образа Лапласа дискретного сигнала может быть выражено следующим образом:
2π
Xd (σ + j(ω +
(5.3)
k)) = Xd (σ + jω), k = ±1, ±2, . . . ,
Td
Это свойство легко доказывается путем подстановки s = σ + j(ω + 2π
Td k)
в полученное ранее выражение для преобразования Лапласа от дискретного
сигнала:
∑
2π
2π
−(σ+j(ω+ T
k))nTd
d
Xd (σ + j(ω +
k)) =
x(nTd )e
=
Td
n
∑
2π
−j(ω+ T
k)nTd
d
x(nTd )e−σnTd e
=
=
n
=
∑
x(nTd )e−σnTd e−jωnTd e−j2πkn =
n
=
=
∑
∑
x(nTd )e−σnTd e−jωnTd =
n
x(nTd )e−(σ+jω)nTd = Xd (σ + jω).
(5.4)
n
79
а)
б)
Рис. 5.1. Примеры образов по Лапласу для аналогово и дискретного сигналов: а–– образ
Лапласа для аналогового сигнала; б–– образ Лапласа для дискретного сигнала
1.2. Z-преобразование
Использование и синтез цифровых и дискретных фильтров на практике
осложнено свойством периодичности их спектра (5.3) так как потребовало
бы обработки неограниченного количества значений. Для решения данной
проблемы вводится Z-преобразование [19], [20], [21], отображающее значения из комплексной s-плоскости в комплексную z-плоскость путем подстановки:
z = esTd .
(5.5)
При использовании данной подстановки преобразование Лапласа становится Z-преобразованием следующего вида:
∑
Xd (z) =
x(nTd )z −n .
(5.6)
n
80
Данное отображение на z-плоскость позволяет преобразовать все бесконечные периодические повторения нулей и полюсов передаточной характеристики в одну точку на плоскости z:
e(σ+jω)Td = e
2π
(σ+jω+ T
k)Td
d
= z.
(5.7)
При этом, в зависимости от вида значения комплексной частоты s,
отображение в z-плоскость будет происходить следующим образом (пример
отображения некоторых значений на z-плоскость представлен на рис. 5.2):
– если s = 0, то z = 1;
– если s = σ (т.е. чисто вещественно), то z = σTd - чисто вещественное и
неотрицательное значение (z ≥ 0). При этом если σ < 0 → z > 1 , а так
же при σ ≥ 0 → z ≥ 1;
– если s = jω (т.е. чисто мнимо), то z = ejωTd - точка, расположенная на
единичной окружности, повернутая на угол ωTd ;
– если s = σ + jω, то z = eσTd ejωTd - точка, расположенная на окружности
радиуса eσTd и повернутая на угол ωTd .
Рис. 5.2. Пример отображения значений из комплексной s плоскости в z-плоскость
81
Описанное преобразование тесно связано с преобразованиями Фурье и
Лапласа, поэтому его свойства схожи со свойствами последних. Незначительные различия в свойствах обусловлены спецификой обрабатываемых
дискретных сигналов. Выделим следующие свойства Z-преобразования.
– Линейность: Z-образ Y (z) линейной комбинации двух сигналов y(n) =
a1 x1 (n) + a2 x2 (n) равен линейной комбинации их образов:
∑
∑
∑
Y (z) =
y(n)z −n =
a1 x1 (n)z −n +
a2 x2 (n)z −n =
n
n
n
= a1 X1 (z) + a2 X2 (z).
(5.8)
– Задержка: Z-образ Y (z) дискретного сигнала x(n − m), задержанного на
m отсчетов, равен произведению z-образа исходного сигнала X(z) и коэффициента z −m , называемого оператором задержки:
∑
∑
Y (z) =
x(n−m)z −n = z −m
x(n−m)z −(n−m) = z −m X(z). (5.9)
n
n
– Z-преобразование свертки: Z-образ Y (z) свертки двух дискретных сиг∑
налов y(n) =
x1 (m)x2 (n − m) равен произведению их Z-образов:
m
Y (z) =
=
∑
m
∑
y(n)z −n =
n
x1 (m)
∑
(
∑ ∑
n
)
x1 (m)x2 (n − m) z −n =
m
x2 (n − m)z −n =
n
∑
x1 (m)X2 (z)z −m =
m
= X1 (z)X2 (z).
(5.10)
1.3. Способы описания дискретных фильтров
1.3.1. Разностное уравнение
Процедуру фильтрации дискретного сигнала можно обобщить вычислением взвешенной суммы от некоторого количества входных и предыдущих
выходных отсчетов:
y(k) = b0 x(k)+b1 x(k −1)+. . .+bm x(k −m)−a1 y(k −1)−. . .−an y(k −n),
(5.11)
где ai и bj –– некоторые вещественные весовые коэффициенты.
82
При этом, если сгруппировать слагаемые так, чтобы слева от знака равенства находились только выходные значения, а справа только выходные
значения, получим запись, называемую разностным уравнением:
y(k)+a1 y(k −1)+. . .+an y(k −n) = b0 x(k)+b1 x(k −1)+. . .+bm x(k −m).
(5.12)
Форма записи в виде разностного уравнения схожа с дифференциальным
уравнением для аналоговых систем (4.8), где операции дифференцирования
были заменены на операторы задержек. Остальные формы записи опираются
на него.
1.3.2. Функция передачи
Дискретные системы могут быть описаны множеством способов, схожих с теми, что использовались для аналоговых сигналов. Главным отличием описания является использование образов Z-преобразования вместо Фурье и Лапласа. Например, по аналогии с тем, как это было для аналоговых
сигналов, выходной сигнал является дискретной сверткой входного сигнала
и импульсной характеристики системы:
∑
y(k) =
x(m)h(k − m),
(5.13)
m
где h(·) –– импульсная характеристика дискретной системы.
Тогда, если применить Z-преобразование к данной формуле (5.13), можно выписать выражения для функции передачи дискретного фильтра H(z):
Y (z) = X(z)H(z), → H(z) =
∑
Y (z)
h(k)z −k .
=
X(z)
(5.14)
k
Для того, чтобы найти отношение образов z-преобразований выходного
и входного сигналов, применим Z-преобразование к разностному уравнению
(5.12):
Y (z)(1+a1 z −1 +a2 z −2 +. . .+an z −n ) = X(z)(b0 +b1 z −1 +b2 z −2 +. . .+bm z −m ).
(5.15)
83
Откуда H(z) = Y (z)/X(z) выражается, как:
M
∑
H(z) =
bm z
m=0
N
∑
1+
−m
.
an
(5.16)
z −n
n=1
Построение частотных характеристик (АЧХ и ФЧХ) для заданной передаточной функции происходит по аналогии с построением характеристик
для аналоговых фильтров за одним исключением - аргумент функции находится с помощью подстановки Z-преобразования:
K̇(ω) = Ḣ(ejωTd ).
(5.17)
Из данной формулы следует, что частотные характеристики дискретной системы, так же как и ее частотный спектр, является периодической функцией
с периодом, равным частоте дискретизации ωd = 2π/Td .
1.3.3. Нули и полюса
По аналогии со случаем обработки аналоговых сигналов, передаточная функция дискретной системы может быть задана как набор нулей и
полюсов. При этом, если раскрыть произведение нулей и полюсов, получим передаточную функцию в виде отношение двух многочленов –– дробнорациональную функцию[3]:
M
∑
H(z) =
bm z
m=0
N
∑
−m
an z −n
1+
∏M
= k ∏m=1
N
(1 − zm z −1 )
n=1 (1
− pn z −1 )
,
(5.18)
n=1
где под k = b0 понимается коэффициент усиления; {zi }–– множество нулей функции передачи; {pi }–– полюсы функции передачи. По аналогии с
свойствами нулей и полюсов для аналоговой системы, функция передачи
дискретной системы H(zi ) = 0 в точках, соответствующих ее нулям, и
H(pi ) → ∞ в точках полюсов.
84
1.3.4. Структуры дискретных фильтров
В соответствии с разностным уравнением дискретного фильтра (5.12),
процедура фильтрации происходит как взвешенное суммирование текущего и предыдущих входных значений. При этом, если во взвешенной сумме используются так же и предыдущие выходные значения процедуры, то
такой фильтр называется рекурсивным. В противном случа, когда предыдущие выходные значения не используются в фильтрации, фильтр называется
нерекурсивным. Мы будем рассматривать структурную схему рекурсивного
фильтра, так как нерекурсивный является его частным случаем.
Для реализации подобного фильтра потребуются элементы следующих
видов:
– операторы умножения (рис. 5.3а): позволяют вычислять произведения
слагаемых взвешенной суммы и соответствующих им коэффициентов (весов);
– сумматоры (рис. 5.3в): используются для сложения элементов из взвешенной суммы;
– линии задержки (рис. 5.3б): данные элементы аналогичны оператору задержки, применяются для получения предыдущих значений сигнала и
обозначаются, как z −n .
а)
б)
в)
Рис. 5.3. Примеры базовых элементов, применяемых для реализации дискретного фильтра:
а–– оператор умножения; б–– линия задержки; в–– сумматор
Реализация дискретного фильтра в соответствии с выведенной ранее
формулой (5.12) с помощью представленных ранее структурных элементов
показана на рис. 5.4 и называется прямой реализацией.
Как видно из приведенной схемы, вычисление выходного значения можно представить в виде двух операций: вычисления рекурсивной и нерекурсивной частей. Тогда, если результат получен после последовательного при85
X(k)
b0
+
+
Y(k)
z-1
z-1
-a1
b1
z-1
...
z-1
-a2
b2
z-1
...
z-1
bN
-aN
Рис. 5.4. Прямая форма фильтра
менения ряда линейных стационарных блоков вычислений, результат не изменится от порядка вычислений. Поменяем местами блоки вычисления рекурсивной и нерекурсивной частей. Из рисунка видно, что имеются две линии задержки, на вход которых подаются одинаковые значения. Если объединить две линии задержки, получим итоговую форму, называемую канонической (рис. 5.5).
Прямая и каноническая формы полностью эквивалентны с точки зрения
получаемого результата.Однако при реализации канонической формы требуется в два раза меньшее количество линий задержек.
1.4. Способы синтеза дискретных фильтров
1.4.1. Синтез дискретного фильтра по аналоговому прототипу
Одним из основных способов синтеза дискретного фильтра является
его синтез на основе аналогового прототипа. Популярность данного способа
объясняется достаточно хорошим развитием теории аппроксимации аналоговых фильтров с идеальной формой АЧХ. Сам синтез подразумевает первичных расчет функции передачи аналогового фильтра H(s) с последую86
X(k)
+
b0
+
Y(k)
z-1
b1
-a1
z-1
-a2
...
b2
z-1
-aN
bN
Рис. 5.5. Каноническая форма фильтра
щим ее преобразованием в функцию передачи дискретного фильтра H(z).
Методы вычисления H(s) уже были детально рассмотренны ранее в подразделе 1.4.2.
Нужно отметить, что при переходе от функции передачи аналогового
фильтра к дискретной, получаемый дискретный фильтр по своим параметрам не является идентичным аналоговому. Данный факт объясняется как минимум периодичностью частотных характеристик дискретного фильтра.
Для получения функции передачи дискретного фильтра из функции аналогового H(s) можно было бы воспользоваться следующей замены ее комплексного аргумента: z = exp(sTd ), где Td - период дискретизации сигнала. Но данная замена является крайне неудобной, так как получаемая H(z)
должна иметь дробно-рациональных характер, так же как и исходная H(s).
Если бы подстановочная функция имела дробно-рациональные характер,
процедура получения H(z) могла бы быть реализована с помощью известного алгоритма, подробно описанного в последующем подразделе 1.4.2. По
87
этому на практике вместо экспоненциального отображения в подстановке используется следующая дробно-рациональная функция, называемая билинейным Z-преобразованием:
s=
2 z−1
1 − z −1
= 2fd
.
Td z + 1
1 + z −1
(5.19)
Тогда z можно выразить через s следующим образом:
z=
2 + sTd
.
2 − sTd
(5.20)
Данная функция (5.19) является некоторой аппроксимацией Zпреобразования, получаемой из него с помощью применения к нему
разложения в ряд Тейлора при ограничении степени ряда равной единицы:
exp(x) = 1 + x +
⇒ z = exp(sTd ) =
∞
∑
x2
xn
x3
+
+ ··· =
,
2
6
n!
n=0
exp(sTd /2)
1 + sTd /2
≈
.
exp(−sTd /2)
1 − sTd /2
(5.21)
1.4.2. Алгоритм вычисления дробно-рациональной подстановки
В предыдущих разделах был рассмотрен способ синтеза дискретного фильтра на основе имеющейся функции передачи аналогового фильтра
H(s), задающейся набором коэффициентов {ai } и {bi }:
M
∑
H(s) =
bm pm
m=0
N
∑
.
(5.22)
an pn
1+
n=1
Тогда вычисление функции передачи дискретного фильтра H(z) происходит путем подстановки некоторой дробно-рациональной функции R(z) в
H(s):
M
∑
bm Rm (z)
m=0
H(z) = H(R(z)) =
,
(5.23)
N
∑
1+
an Rn (z)
n=1
88
где R(z) определено, как:
K
∑
B(z)
R(z) =
= k=0
Q
A(z)
∑
βk z k
.
αq
(5.24)
zq
q=0
Для алгоритма вычисляющего подстановку важно, чтобы степени многочленов числителя и знаменателя были одинаковы для исходной функции
H(s) и функции подстановки R(z), т.е. N = M и K = Q. На практике, даже
если эти равенства не выполняются, можно считать недостающие коэффициенты числителя или знаменателя равными нулям.
Подставим R(z) в H(s):
N
∑
H(R(z)) =
m=0
N
∑
n=1
(
bm
(
an
B(z)
A(z)
B(z)
A(z)
)m
)n =
2
N
2
N
B (z)
B (z)
b0 + b1 B(z)
A(z) + b2 A2 (z) + · · · + bN AN (z)
B (z)
B (z)
1 + a1 B(z)
A(z) + a2 A2 (z) + · · · + aN AN (z)
=
b0 AN (z) + b1 AN −1 (z)B(z) + b2 AN −2 (z)B 2 (z) + · · · + bN B N (z)
=
AN (z) + a1 AN −1 (z)B(z) + a2 AN −2 (z)B 2 (z) + · · · + aN B N (z)
N
N
K
∑
∑
bm AN −m (z)B m (z)
dm z m
m=0
m=0
=
= NK
.
N
∑
∑
N
−n
n
n
N
an A
(z)B (z)
cn z
A (z) +
=
n=1
n=0
(5.25)
Из полученной формулы видно, что искомые коэффициенты передаточной функции ⃗c и d⃗ могут быть получены путем вычисления следующих операций.
1. Вычисление коэффициентов многочленов A(z) и B(z), возведенных в
различные степени.
2. Вычисление коэффициентов всех N произведений многочленов вида
AN −n (z)B n (z), при n = 0, 1, · · · , N .
3. Домножение каждого полученного произведения на соответствующий
ему коэффициент числителя bi или знаменателя bj .
4. Приведение подобных слагаемых.
Алгоритм получает в качестве входных данных:
89
– вектор ⃗a = {a0 , a1 , · · · , aN }, содержащий коэффициенты знаменателя
H(s). Следует отметить, что коэффициент a0 = 1;
– вектор ⃗b = {b0 , b1 , · · · , bN }, содержащий коэффициенты числителя H(s);
– вектор α
⃗ = {α0 , α1 , · · · , αK }, содержащий коэффициенты знаменателя
R(z);
⃗ = {β0 , β1 , · · · , βK }, содержащий коэффициенты числителя
– вектор β
R(z);
– матрицу POWER_A, из N строк и (N − 1)(K − 1) + 1 столбцов для вычисления коэффициентов многочлена A(z) в различных степенях;
– матрицу POWER_B, из N строк и (N − 1)(K − 1) + 1 столбцов для вычисления коэффициентов многочлена B(z) в различных степенях.
Инициализация матриц: все элементы матриц, кроме элемента первой
строки и первого столбца, который инициализируется ’1’, инициализируются нулями:


1 0 0 ··· 0
0 0 0 · · · 0

P OW ER_A = P OW ER_B = 
(5.26)
. . . . . . . . . . . . . . . .
0 0 0 ··· 0
Шаг 1. Последовательно применим к каждой строке соответствующей
матрицы следующую операцию для вычисления коэффициентов многочленов, возведенных в различные степени:
⃗
P OW ER_B(i) = P OW ER_B(i − 1) ∗ β;
P OW ER_A(i) = P OW ER_A(i − 1) ∗ α
⃗,
(5.27)
где под оператором ∗ подразумевается линейная свертка.
Шаг 2. Введем еще одну матрицу M U LT _AB аналогичной размерности для хранения результатов произведений многочленов AN −n (z)B n (z) и
заполним ее, вычислив каждую ее строку по следующей формуле:
M U LT _AB(i) = P OW ER_A(N − i) ∗ P OW ER_B(i).
(5.28)
Шаг 3-4. Последние два шага могут быть выполнены умножением вектора с соответствующими коэффициентами на полученную ранее матрицу
коэффициентов произведений многочленов:
90
d⃗ = ⃗b ∗ M U LT _AB,
(5.29)
⃗c = ⃗a ∗ M U LT _AB,
(5.30)
где d⃗ и ⃗c –– вектора, содержащие искомые коэффициенты многочленов числителя и знаменателя передаточной функции дискретного фильтра H(z).
2. Порядок самостоятельного выполнения исследования
В ходе выполнения исследования необходимо выполнить следующие задания.
1. Использовать в качестве исходных данных параметры аналогового фильтра Баттерворта, реализованного в предыдущем исследовании. Рассчитать
передаточную характеристику H(s) аналогового фильтра.
2. Используя передаточную характеристику H(s) аналогового фильтра, рассчитать передаточную характеристику H(z) (т. е. набор коэффициентов
ai и bi ) цифрового фильтра Баттерворта, воспользовавшись билинейным
преобразованием (5.19).
3. Построить АЧХ и ФЧХ реализованного фильтра, используя подстановку
z = ejω .
4. Продемонстрировать работу фильтра при обработке гармонического сигнала, представленного в цифровой форме. Для этого написать программу фильтрации цифрового сигнала с помощью реализованного фильтра
Баттерворта. В программе использовать каноническую схему цифрового
фильтра с коэффициентами ai и bi , вычисленными ранее.
91
ЛИТЕРАТУРА
1. А.К.Лосев. Теория линейных электрических цепей. М.: Высшая школа,
1987. С. 527.
2. Ю.Сато. Без паники! Цифровая обработка сигналов. Додэка-XXI, 2010.
С. 176.
3. А.Б.Сергиенко. Цифровая обработка сигналов. СПб.: БХВ-Петербург,
2011. С. 768.
4. В.А.Зорич. Математический анализ, Часть II. М.: Наука, 1984. С. 640.
5. В.А.Котельников. О пропускной способности эфира и проволоки в электросвязи // Материалы к I Всесоюзному съезду по вопросам технической реконструкции дела связи и развития слаботочной промышленности. 1933.
6. Shannon C. E. Communication in the Presence of Noise // Proc. Institute of
Radio Engineers. 1949. Т. 37, № 1. С. 10–21.
7. E.Bocharova. Compression for Multimedia.
Press, 2010.
NY.: Cambrige University
8. Л.Шапиро Дж. Стокман. Компьютерное зрение. М.: БИНОМ. Лаборатория знаний, 2013. С. 753.
9. Р.Гонсалес Р. Вудс. Цифровая обработка изображений. М.: Техносфера,
2005. С. 1072.
10. У.Прэтт. Цифровая обработка изображений. Том 1.
С. 311.
М.: Мир, 1982.
11. J.Max. Quantizing for minimum distortion // IEEE Transactions on
Information Theory. 1960. Т. 6(1).
12. S.P.Lloyd. Least squares quantization in PCM // IEEE Transactions on
Information Theory. 1982. Т. 28(2).
13. Y.Linde A.Buzo R.M.Gray. An algorithm for vector quantizer design // IEEE
Transactions on Communications Technology. 1980. Т. COM-28(1).
14. Б.Д.Кудряшов. Теория информации. СПб.: Питер, 2009.
15. У.Прэтт. Цифровая обработка изображений. Том 2.
С. 477.
М.: Мир, 1982.
16. С.И.Баскаков. Радио/технические цепи и сигналы. М.: Высшая школа,
2000. С. 462.
92
17. S. Lipschutz M. Spiegel J. Liu. Mathematical Handbook of Formulas and
Tables. Schaum’s Outline Series (3rd ed.), 2009. С. 183.
18. Бахурин С. Теория и практика цифровой обработки сигналов.
http://www.dsplib.ru. Доступно на момент: 2016-11-16.
19. Л.Р.Рабинер Б. Голд. Теория и применение цифровой обработки сигналов. М.: Мир, 1978. С. 848.
20. А.В.Оппегейм Р.В. Шафер. Цифровая обработка сигналов. М.: Связь,
1979. С. 416.
21. Л.Р.Рабинер Р.В. Шафер. Цифровая обработка речевых сигналов. М.:
Радио и связь, 1981. С. 496.
93
СОДЕРЖАНИЕ
Предисловие
3
I. Дискретизация аналогового сигнала. Фурье-анализ сигналов
1
Теоретические основы . . . . . . . . . . . . . . . . . . . . . . .
1.1
Основные понятия . . . . . . . . . . . . . . . . . . . . . . .
1.2
Классификация сигналов . . . . . . . . . . . . . . . . . . .
1.3
Энергия и мощность сигнала . . . . . . . . . . . . . . . . .
1.4
Способы математического описания сигналов . . . . . . .
1.5
Фурье анализ . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6
Интегральное преобразование Фурье . . . . . . . . . . . .
1.7
Постановка задачи аналого-цифрового и цифроаналогового преобразования . . . . . . . . . . . . . . . . .
1.8
Дискретизация аналогового сигнала. Теорема Котельникова
1.9
Восстановление аналогового сигнала по дискретным отсчетам . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
Порядок самостоятельного выполнения исследования . . . . .
II. Дискретное преобразование Фурье
1
Теоретические основы . . . . . . . . . . . . . . . . . .
1.1
Дискретное преобразование Фурье . . . . . . . . .
1.2
Спектр дискретного сигнала . . . . . . . . . . . .
1.3
Растекание спектра . . . . . . . . . . . . . . . . . .
2
Порядок самостоятельного выполнения исследования
4
4
4
5
7
9
15
19
23
24
27
28
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
32
32
32
34
36
37
III. Методы квантования в системах обработки сигналов
1
Теоретические основы . . . . . . . . . . . . . . . . . . .
1.1
Классификация методов квантования . . . . . . . .
1.2
Критерии оценки эффективности квантования . . .
1.3
Скалярное квантование . . . . . . . . . . . . . . . .
1.4
Векторное квантование. Алгоритм Линде-Бузо-Грея
1.5
Использование палитры в формате BMP . . . . . . .
2
Порядок самостоятельного выполнения исследования .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
39
39
40
42
53
56
57
94
IV. Аналоговые линейные системы
1
Теоретические основы . . . . . . . . . . . . . . . . . .
1.1
Основные определения . . . . . . . . . . . . . . .
1.2
Основные характеристики линейных систем . . .
1.3
Способы описания линейных систем . . . . . . . .
1.4
Аналоговые фильтры . . . . . . . . . . . . . . . .
2
Порядок самостоятельного выполнения исследования
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
V. Дискретные линейные системы
1
Теоретические основы исследования . . . . . . . . . . . . . .
1.1
Особенности преобразования Лапласа для дискретного
сигнала . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2
Z-преобразование . . . . . . . . . . . . . . . . . . . . . .
1.3
Способы описания дискретных фильтров . . . . . . . . .
1.4
Способы синтеза дискретных фильтров . . . . . . . . . .
2
Порядок самостоятельного выполнения исследования . . . .
Литература
.
.
.
.
.
.
60
60
60
61
62
67
76
78
. 78
.
.
.
.
.
78
80
82
86
91
92
95
Учебное издание
Гильмутдинов Марат Равилевич,
Егоров Николай Дмитриевич,
Веселов Антон Игоревич
ВВЕДЕНИЕ В ЦИФРОВУЮ ОБРАБОТКУ СИГНАЛОВ
Учебное пособие
Публикуется в авторской редакции
Подготовка к печати А. Н. Колешко
Сдано в набор 15.12.2016. Подписано к печати 19.12.2016.
Формат 60 × 84 1/16. Усл. печ. л. 5, 5. Тираж 50 экз. Заказ №498.
Редакционно-издательский центр ГУАП 190000,
Санкт-Петербург, Б. Морская ул., 67
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 645 Кб
Теги
gilmutdinov
1/--страниц
Пожаловаться на содержимое документа