МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РЕСПУБЛИКИ КАЗАХСТАН

 Некоммерческое акционерное общество

Алматинский университет энергетики и связи

 

 

 

Л.Г. Богомолова, Г.С. Казиева

 

Цифровая обработка сигналов в телекоммуникационных системах.

 

Учебное пособие

 

 

Алматы 2011


УДК 004.383.3:621.394/397(075.8)

ББК32.88.5я73

К 14 Цифровая обработка сигналов в телекоммуникационных  системах:

Учебное пособие / Г.С. Казиева, Л.Г. Богомолова;

Алматы: АУЭС, 2011 – 81 с.

 

 

ISBN 978-601-7307-16-8 

 

                Учебное пособие излагаются как теоретические вопросы построения цифровых фильтров, так и практические примеры построения данных систем. 

    Учебное пособие предназначено для студентов  дневного и заочного факультетов специальностей специальности 5В071900 – Радиотехника, электроника и телекоммуникации.

 

Табл. 3, Ил.33, Библиогр. - 18 назв.

 

 

ББК32.88.5я73

 

 

РецензентЫ: 1. к.ф.м.н., доцент Каз НТУ, А.К. Шайкин

                            2.к.т. н., доцент АУЭС,  К.С. Чежимбаева 

 

Печатается по плану издания Министерства образования и науки Республики Казахстан на 2011 г.

 

 

ISBN 978-601-7307-16-8

 

 

©Некоммерческое АО «Алматинский университет энергетики и связи», 2011г.


Содержание

 

Введение

4

ГЛАВА 1. Основы цифровой обработки сигналов 

4

1.1 Дискретные и цифровые сигналы

4

1.2 Система цифровой обработки сигналов.

7

1.3 Преимущества и недостатки цифровых фильтров.

8

ГЛАВА 2. Дискретизация и восстановление непрерывных сигналов 

9

2.1 Выбор частоты дискретизации. Спекторы  дискретизованных сигналов

9

2.2  Восстановление непрерывного сигнала. Погрешности

дискретизации и восстановления сигнала  

10

ГЛАВА 3. Преобразование дискретных сигналов  

13

3.1 Дискретное преобразование Фурье 

13

3.2 Дискретная свертка 

16

3.3 z-преобразование. Основные свойства z-преобразования 

17

ГЛАВА 4. Характеристики цифровых фльтров  

23

4.1 Импульсная характеристика  

23

4.2 Алгоритм цифровой фильтрации 

24

4.3 Нерекурсивные и рекурсивные цифровые фильтры  

24

4.4 Системная (передаточная) функция 

26

4.5 Формы реализации цифровых фильтров  

27

ГЛАВА 5. Основы синтеза цифровых фильтров  

40

5.1 Метод инвариантной импульсной характеристики 

40

5.2 Метод билинейного z-преобразования 

42

5.3 Метод частотной выборки 

45

5.4 Метод взвешивания с помощью окна 

46

ГЛАВА 6. Деконволюция цифровых фильтров 

48

6.1. Понятие деконволюции 

48

6.2. Инверсия импульсного отклика фильтра 

51

6.3. Оптимальные фильтры деконволюции 

52

ГЛАВА 7. Аппрксимация сигналов и функций

58

7.1 Приближение сигналов рядами Тейлора 

59

7.2.  Интерполяция и экстраполяция сигналов 

60

7.3.  Сплайновая интерполяция  

64

7.4.  Спектральный метод  интерполяции   

65

ГЛАВА 8. Регрессия 

75

8.1.  Постановка задачи регрессии  

75

8.2.  Линейная  регрессия  

76

8.3.  Полиномиальная  регрессия  

77

Список литературы 

81

 

 ВВЕДЕНИЕ 

         Методы цифровой обработки сигналов в настоящее время получили повсеместное развитие и в значительной мере заменяют классические аналоговые методы.

Физические величины макромира, как основного объекта наших измерений и источника информационных сигналов, как правило, имеют непрерывную природу и отображаются непрерывными (аналоговыми) сигналами. Цифровая обработка сигналов (ЦОС или DSP - digital signal processing) оперирует исключительно с дискретными величинами, причем с квантованием как по координатам динамики своих изменений (во времени, в пространстве, и любым другим изменяемым аргументам), так и по амплитудным значениям физических величин. Математика дискретных преобразований зародилась в недрах аналоговой математики еще в 18 веке в рамках теории рядов и их применения для интерполяции и аппроксимации функций, однако ускоренное развитие она получила в 20 веке после появления первых вычислительных машин. В принципе, в своих основных положениях математический аппарат дискретных преобразований подобен преобразованиям аналоговых сигналов и систем. Однако дискретность данных требует учета этого фактора, а его игнорирование может приводить к существенным ошибкам. Кроме того, ряд методов дискретной математики не имеет аналогов в аналитической математике. Стимулом быстрого развития дискретной математики является и то, что стоимость цифровой обработки данных меньше аналоговой и продолжает снижаться, даже при очень сложных ее видах, а производительность вычислительных операций непрерывно возрастает. Немаловажным является также и то, что системы ЦОС отличаются высокой гибкостью. Их можно дополнять новыми программами и перепрограммировать на выполнение различных функций без изменения оборудования. В последние годы ЦОС оказывает постоянно возрастающее влияние на ключевые отрасли современной промышленности: телекоммуникации, средства информации, цифровое телевидение и пр. Следует ожидать, что в обозримом будущем интерес и к научным, и к прикладным вопросам цифровой обработки сигналов будет нарастать во всех отраслях науки и техники.

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

   

1 Глава. Основы цифрвой обработки сингалов

 

2.1. Дискретные и цифровые сигналы

 

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

         Большинство сигналов, с которыми приходится иметь дело, являются аналоговыми (непрерывными). Это связано с тем, что сигнал является отображением (аналогом) некоторого физического процесса, а почти все физические процессы непрерывны по своей природе.

         Аналоговыми называются сигналы, произвольные по величине и непрерывные по времени (см. рисунок 1.1, а).

Сигнал, дискретизированный только по времени, называют дискретным. Он еще не пригоден для обработки в цифровом устройстве. Дискретный сигнал представляет собой последовательность, элементы которой f(kT) в точности равны соответствующим значениям исходного аналогового сигнала f(t) (см. рисунок 1.1, б). Аналитически такой дискретный сигнал описывается выражением

,           (1.1)

где:

 - f(t) – исходный аналоговый сигнал;

-  - -функция (бесконечно короткий импульс бесконечно большой амплитуды, обладающий единичной площадью).

Следует заметить, что дискретный сигнал (1.1) – это абстракция: последовательность импульсов бесконечно большой амплитуды, а не  с амплитудами, равными отсчетам аналогового сигнала f(kT). Значения отсчетов аналогового сигнала определяют не амплитуды, а “веса” -функции. Отметим, что графически дискретный  сигнал  принято изображать в виде набора линий с высотами f(kT) (см. рисунок 1.1, б), хотя это и не соответствует математической записи (1.1).

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

Процесс квантования сводится к преобразованию непрерывных по уровню отсчетов аналогового сигнала f(kT) в квантованные значения fц(kT), амплитуды которых совпадают с ближайшими допустимыми уровнями квантования (см. рисунок 1.1, в). Таким образом, квантование представляет собой округление передаваемых значений с заданной точностью. Сигнал, дискретизированный и по времени, и по уровню (квантованный), называют цифровым.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 Рисунок 1.1 Дискретные сигналы


1.2. Система цифровой обработки сигналов. Дискретные и цифровые фильтры

 

         Рассмотрим, из каких элементов должна состоять система цифровой обработки сигналов.

         Практически система цифровой обработки аналоговых сигналов осуществляется в соответствии со структурной схемой, изображенной на рисунке 1.2.

 

 

 

 

 Рисунок   1.2 Структурная схема цифровой обработки сигналов.

        

Аналоговый сигнал x(t) поступает на вход дискретизатора (Д), на выходе которого появляются отсчеты x(kT). Реально, аналого-цифровой преобразователь (АЦП) выдает отсчеты xц(kT) в цифровой форме, т.е. в виде комбинаций двоичного кода, так что каждая комбинация этого кода является представлением соответствующего отсчета в двоичной системе счисления. Точность такого описания дискретных отсчетов сигнала x(t) зависит от числа разрядов в кодовых комбинациях и в принципе может быть сделана сколько угодно высокой.

         Последовательность отсчетов xц(kT) c выхода АЦП поступает на вход цифрового вычислительного устройства, которое по заданному алгоритму “пересчитывает” ее в последовательность выходных отсчетов y(kT). Алгоритмы обработки цифровых сигналов могут быть очень разнообразными как по характеру, так и по степени сложности. Цифровые устройства, производящие линейную обработку сигналов, называют цифровыми фильтрами (ЦФ). Разумеется, выходная последовательность отсчетов ЦФ y(kT) также представляет собой последовательность кодовых комбинаций двоичного кода. Для преобразования цифрового сигнала в аналоговый используют восстанавливающее устройство (штриховая линия на рисунке 1.2), состоящее из цифро-аналогового преобразователя (ЦАП) и выходного сглаживающего фильтра (СФ). ЦАП преобразует цифровой сигнал в импульсы прямоугольной формы, чьи амплитуды в точности равны записанным в двоичной системе счисления значениям y(kT). Импульсы прямоугольной формы подаются на СФ, и на выходе этого фильтра получается аналоговый сигнал y(t).

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

         Дискретные фильтры могут быть реализованы в виде некоторой аналоговой схемы, состоящей из линии задержки, сумматоров и перемножающих элементов.

Математический аппарат теории дискретных фильтров очень полезен

при теоретическом анализе прохождения сигналов через цифровые фильтры, так как позволяет, пренебрегая эффектами квантования, сформулировать и определить все основные параметры цифровых фильтров.

 

1.3. Преимущества и недостатки цифровых фильтров. Реализация цифровых фильтров

 

         Интенсивное развитие цифровой техники не означает, что должны быть полностью отброшены и забыты аналоговые устройства для обработки сигналов. Каждому из типов фильтров присущи свои преимущества и недостатки и в зависимости от конкретных условий следует применять тот или иной тип фильтра.

         По сравнению с аналоговыми фильтрами цифровые фильтры имеют следующие преимущества:

1) стабильность частотных и временных характеристик, малая их зависимость

2) от влажности, температуры, изменения питающих напряжений и т.д.;

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

3) отсутствие этапа предварительной настройки (в аналоговом фильтре нередко приходится подгонять номиналы элементов);

4) значительно более высокая точность обработки сигналов;

5) гибкая оперативная перестройка алгоритмов обработки сигналов, например, адаптивных (подстраивающихся) алгоритмов, изменяющихся при изменении параметров входного сигнала;

6) отсутствует задача согласования нагрузок;

7) уменьшенные масса и габаритные размеры, а также унифицированное и стандартизованное исполнение.

По сравнению с аналоговыми фильтрами цифровые фильтры имеют следующие недостатки:

1)  обязательные аналого-цифровые и цифро-аналоговые преобразователи, что удорожает и усложняет процесс цифровой обработки сигналов;

2) повышенное быстродействие процессора, реализующего ЦФ, аналого-цифровых и цифро-аналоговых преобразователей (особенно при высоких разрядностях входных переменных);

3) наличие специфических погрешностей, вызванных дискретизацией, квантованием сигнала и округлением значений обрабатываемого сигнала в процессе вычислений.

Сравнительный анализ позволяет сделать вывод, что если алгоритм обработки сигналов не сложен и не требует высокой точности, целесообразно применять аналоговые фильтры или дискретные фильтры. Они проще и дешевле, чем цифровые фильтры.

Цифровые фильтры могут быть реализованы в виде как специализированного вычислителя (аппаратурная реализация), так и в виде программ для универсальной ЭВМ (программная реализация). При аппаратурной реализации каждый функциональный узел (сумматор, элемент задержки, умножитель) цифрового фильтра выполняется, как правило, в виде отдельного модуля, например отдельной микросхемы. При программной реализации аналитический алгоритм работы цифрового фильтра представляет собой программу, записываемую в память универсальной ЭВМ.

Основным достоинством аппаратурной реализации цифровых фильтров является повышенное быстродействие; основным недостатком – недостаточная гибкость перестройки, связанная с заменой отдельных узлов цифрового фильтра, а также невозможность использовать данный специализированный вычислитель для других цифровых фильтров.

Основным достоинством программной реализации цифровых фильтров является достаточная гибкость перестройки, которая сводится лишь к вводу в память данной ЭВМ другой программы; основным недостатком – пониженное быстродействие, поскольку структура универсальной ЭВМ ориентирована на решение класса различных задач, а не только на реализацию данного цифрового фильтра.

 

                 2 Глава. Дискретизация и восстановление непрерывных сигналов

 

2.1. Выбор частоты дискретизации. Спектры дискретизированных сигналов

 

         Как указывалось, в результате дискретизации непрерывный сигнал заменяется дискретными отсчетными значениями, взятыми через определенные интервалы времени.

Очень важен правильный выбор интервала дискретизации. Если отсчеты сигнала брать слишком редко, то быстрые изменения, скачки дискретизируемого сигнала могут остаться незамеченными. Если отсчеты брать слишком часто, то это исключит ошибки, но для сохранения общего времени обработки сигнала неизменным приходится увеличивать скорость обработки, что не всегда возможно.

         Задача о выборе интервала дискретизации наиболее просто решается для сигналов с ограниченным спектром на основании теоремы Котельникова, или теоремы отсчетов.

         В соответствии с теоремой Котельникова, непрерывный сигнал f(t), в спектре которого не содержится частот выше fв, полностью описывается выборочными значениями f(kT), отсчитанными через интервалы времени  , где . Аналитически это выражается в виде ряда Котельникова

.         (2.1)

         Интервал времени  между соседними отсчетами называют интервалом Котельникова или интервалом Найквиста.

         В соответствии с теоремой Котельникова вместо непрерывного сигнала f(t) с ограниченным спектром можно передавать дискретную последовательность значений f(kT), причем интервал дискретизации  (, где - частота дискретизации), должен быть не более чем . Если отсчеты взять слишком редко, то это может привести к грубым ошибкам.

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

         Важной характеристикой дискретизированного колебания является его спектр. Наиболее просто определяется спектр дискретизированного сигнала , записанного в виде модулированной последовательности d - функций:

.       (2.2)

Если известен спектр исходного непрерывного сигнала f(t), то спектр дискретизированного сигнала определяется выражением:

,   (2.3)

где  - спектр исходного непрерывного сигнала f(t).

         Из (2.3) следует, что спектр дискретизированного сигнала представляет собой сумму всевозможных сдвигов спектра исходного непрерывного сигнала f(t) вверх и вниз по оси частот  на величины  ± n(2p/T).

 

2.2. Восстановление непрерывного сигнала. Погрешности дискретизации и восстановления сигнала

 

         Если спектр непрерывного исходного сигнала f(t) Sf(w) ограничен по ширине (см. рисунок 2.1, а), а интервал дискретизации T удовлетворяет условию теоремы Котельникова  и, следовательно, период повторения спектра дискретизированного сигнала  , то соседние части спектра, соответствующие Sf[w-n(2p/T)] и Sf[w-(n+1)(2p/T)], не перекрываются  (см. рисунок  2.1, б).

В этом случае легко указать способ восстановления непрерывного сигнала из дискретного, который состоит в том, что дискретный сигнал следует пропустить через идеальный фильтр нижних частот с полосой пропускания  (см. рисунок 2.1, в).

При этом из спектра дискретизированного сигнала  будет выделена средняя часть (см. рисунок 2.1, г), которая с точностью до постоянного множителя 1/T совпадает со спектром исходного непрерывного колебания f(t). Однако, если исходное непрерывное колебание таково, что его спектр с ростом частоты не обращается строго в нуль, то при любом выборе интервала дискретизации соседние составляющие спектра дискретизированного колебания будут частично перекрываться.

Если сигнал с таким спектром пропускать через идеальный фильтр нижних частот, то на выходе фильтра получится колебание, отличающееся от исходного непрерывного сигнала f(t).

Это отличие состоит не только в том, что “отрезана” часть спектра выше частоты , но также и в том, что на спектр этого колебания накладываются “хвосты” от соседних спектральных составляющих.

Наиболее простой и очевидный способ уменьшения ошибки дискретизации - это повышение частоты дискретизации. Увеличение частоты дискретизации раздвигает спектры и уменьшает их перекрытие. Однако, для получения достаточно малой ошибки частоту дискретизации приходится брать очень высокой, особенно, если спектр сигнала убывает медленно, что в ряде случаев бывает нежелательно.

Для уменьшения погрешности дискретизации можно перед дискретизацией пропустить сигнал через фильтр нижних частот с частотной характеристикой, близкой к прямоугольной. При этом спектр сигнала становится быстро убывающим, почти ограниченным, и дальнейшая дискретизация происходит практически без ошибок. Результирующая ошибка в этом случае определяется искажениями спектра при прохождении сигнала через фильтр нижних частот. Вследствие того, что на спектр сигнала в области частот  не накладываются “хвосты” от соседних составляющих, эта ошибка получается приблизительно в 2 раза меньше, чем при непосредственной дискретизации сигнала.

 

 
Пропускание сигнала через фильтр нижних частот перед дискретизацией является очень полезной мерой для снижения погрешности, если дискретизация сигнала производится при наличии широкополосного шума на входе. При прохождении через фильтр нижних частот дисперсия шума уменьшается и соответственно уменьшается ошибка дискретизации.

Еще одним источником ошибки является неидеальная фильтрация в процессе восстановления непрерывного сигнала из дискретного. Идеальная прямоугольная форма частотной характеристики фильтра нижних частот практически не может быть реализована; для сглаживания сигнала обычно применяют фильтры,    имеющие  монотонно    спадающую  характеристику

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

  

 

 

 

 

 

 

 Рисунок 1.1 - Спектры сигналов


Главный метод борьбы с этими погрешностями состоит в увеличении частоты дискретизации. Однако, увеличение частоты дискретизации приводит к усложнению и удорожанию устройства обработки сигналов. Поэтому в каждом конкретном случае приходится искать компромиссное решение, исходя из характера сигнала, требуемой точности его восстановления,

характеристик применяемого сглаживающего фильтра и других факторов. Все это приводит к тому, что в реальных устройствах частота дискретизации выбирается равной не 2fв , как следует из теоремы Котельникова, а в несколько раз большей. Обычно частота дискретизации (4...10)fв.

 

3 Глава. Преобразование дискретных сигналов

 

         3.1. Дискретное преобразование Фурье

 

         При теоретическом исследовании процессов цифровой фильтрации квантованием сигналов по уровню обычно пренебрегают, т.е. при таком анализе вместо цифровых сигналов рассматривают дискретные.

Таким образом,  дискретные неквантованные сигналы играют важную роль в теории цифровой фильтрации. Поэтому в дальнейшем будут рассмотрены только дискретные сигналы.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


Методы математического описания дискретных сигналов во многом аналогичны методам математического описания непрерывных сигналов. Обычному преобразованию Фурье соответствует дискретное преобразование Фурье, преобразованию Лапласа соответствует дискретное преобразование Лапласа и z-преобразование, свертке двух непрерывных сигналов соответствует дискретная свертка.

         Спектр дискретного сигнала  можно найти не только по спектру непрерывного сигнала, но и непосредственно по самим отсчетам f(nT), применяя прямое преобразование Фурье к сигналу Y(t):

.                                 (3.1)

         В свою очередь, согласно обратному преобразованию Фурье

.                       (3.2)

         Для периодических последовательностей с периодом NT вводится пара преобразований, называемая дискретным преобразованием Фурье (ДПФ):

, k=0,1,...,N-1;        (3.3)

, n=0,1,...,N-1,    (3.4)

где , основная частота периодической последовательности, причем соотношение (3.3) носит название прямого дискретного преобразования Фурье (ДПФ), а (3.4) – обратного дискретного преобразования Фурье (ОДПФ).

         При подстановке значения W в (3.3) и (3.4), в результате сокращения, длительность интервала дискретизации T выпадает из рассмотрения, а ДПФ и ОДПФ записывают в форме:

, k=0,1,...,N-1;           (3.5)

, n=0,1,...,N-1.          (3.6)

         Отметим, что последовательность , как и последовательность , периодична с периодом N отсчетов, так как , где m-целое число.

         Прямое ДПФ и ОДПФ могут быть использованы и для последовательности конечной длины N, определенной при n=0,1,2,...,N-1 и равной нулю при других n. Действительно, мы можем представить такую последовательность конечной длины N периодической последовательностью периода N, один период которой совпадает с данной последовательностью. В этом случае исходная последовательность имеет такое же ДПФ, что и периодическая последовательность, поскольку можно вычислить с использованием ДПФ один период периодической последовательности, а следовательно, и последовательность конечной длины. Следует только считать, что при n<0 и n>N-1 периодическая последовательность равна нулю.

         Несмотря на общность, есть существенное различие между обычным преобразованием Фурье и ДПФ. Заметим, что если сравнивать спектр конечного дискретного сигнала, определяемый формулой (3.1) (с учетом того, что f(nT)=0 при n<0 и n>N-1), и ДПФ этого сигнала (3.3), то ясно, что ДПФ представляет собой N отсчетов спектра, взятых на одном периоде  с интервалом дискретизации по частоте, равным .

         Пример 3.1. Пусть конечный дискретный сигнал  определяется в виде ; ; ; . Вычислить дискретное преобразование Фурье.

Для нашего случая (N=4), вычислив промежуточные величины , , , получим

;

;

;

.

         Рассмотрим один прием работы с ДПФ, который является одним из наиболее полезных. Пусть задан конечный дискретный сигнал , определенный при n=0,1,2,...,N-1. Отсчеты спектра - это ряд линий, “вписанных” в  в точках . Образуем новый дискретный сигнал , определенный при

n=0,1,2,...,L-1(L>N):

       (3.7)

Физически сигнал  не изменится, следовательно,  также останется прежним. Отсчеты спектра  будут другими, что очевидно из формулы ДПФ (3.3). Линии спектра будут “вписаны” в  в точках , k=0,1,2,...,L-1, т.е. более часто.

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

         Пример 3.2. Пусть конечный дискретный сигнал  определяется в виде ; . Вычислите дискретное преобразование Фурье для четырех отсчетов спектра.

         Образуем новый дискретный сигнал  в виде ; ; ; , который получен из  путем дополнения его недостающим числом нулевых отсчетов. Тогда, использовав ответ примера 3.1, находим ; ; ; .

 

3.2. Дискретная  свертка

 

         Если x(nT) и y(nT) – две периодические последовательности с периодами N, то можно ввести понятие о круговой (или периодической) свертке, которую получают как  значения периодической функции

       (3.8)

при k=0,1,2,...,N-1.

         Применение к обеим частям равенства (4.8) ДПФ приводит к соотношению:

,                         (3.9)

где ,  и  - ДПФ последовательностей f(nT), x(nT) и y(nT).

         Полученные формулы справедливы и для конечных последовательностей, если рассматривать x(nT) и y(nT) как эквивалентные им периодические последовательности с теми же ДПФ. В дальнейшем при работе с ДПФ будем сразу считать  заданные последовательности конечной длины периодически продолженными на всю временную ось.

         Как правило, для конечных последовательностей обычно нужна линейная, а не круговая свертка, поэтому в приведенные формулы следует внести уточнения.

         В общем случае нам может понадобиться свертка двух последовательностей неравной длины. Пусть x(nT) и y(nT) – две конечные последовательности с периодами и , т.е. x(nT) отлична от нуля при , а y(nT) – при . Линейной сверткой этих последовательностей называют последовательность f(nT), определяемую соотношением

.                    (3.10)

Нетрудно видеть, что длина  последовательности f(nT) равна сумме длин сворачиваемых последовательностей, точнее .

         Вычисление свертки непосредственно по формуле (3.10) при больших длинах и  очень громоздко. В том случае, когда и  превышают 50...100 отсчетов, значительно эффективнее путь, использующий ДПФ. Алгоритм вычисления круговой свертки может быть применен и для вычисления линейной свертки. Свертка периодических последовательностей периодична и имеет тот же период, что и сами последовательности.  Поскольку период свертки f(nT) равен  отсчетам, то для получения такого периода при круговой свертке необходимо, чтобы x(nT) и у(nT) содержали по отсчетов, что достигается дополнением каждой из двух последовательностей соответствующим числом нулевых отсчетов. Теперь можно вычислить в  точках ДПФ дополненных последовательностей,

перемножить их и выполнить ОДПФ произведения. В результате получаются

значения отсчетов свертки f(nT).

         Пример 3.3. Пусть конечные дискретные сигналы x(nT) и у(nT) с периодами и  определяются в виде x(0)=2; x(T)=1; y(0)=1; y(T)=-1; y(2T)=2. Рассмотрим вычисление линейной свертки для этих сигналов.

         Непосредственный подсчет по формуле (3.10) для нашего случая () приводит к следующему результату:

f(0)=x(0)y(0)=2×1=2;

f(T)=x(T)y(0) + x(0)y(T)=1×1+2(-1)=-1;

f(2T)=x(T)y(T) + x(0)y(2T)=1(-1)+2×2=3;

f(3T)=x(T)y(2T)=1×2=2.

         Одно из основных применений ДПФ – это вычисление спектров функций, заданных графически или таблично.

         Другое важное применение ДПФ – вычисление сигнала на  выходе фильтра с заданной частотной характеристикой. Если задан входной сигнал x(kT), то для него можно вычислить ДПФ . Если теперь умножим  на частотную характеристику фильтра, то получим ДПФ выходного сигнала: . После этого, с помощью ОДПФ, можно найти сигнал на выходе фильтра.

         Если входной сигнал имеет большую длительность, его обработку с помощью ДПФ можно производить по частям. Для этого берут первые N отсчетов входного сигнала, вычисляют их ДПФ и после умножения на частотную характеристику фильтра с помощью ОДПФ вычисляют первые N отсчетов выходного сигнала. После этого аналогичным путем обрабатывают следующие N отсчетов входного сигнала и т.д.

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

         Обработку сигналов с помощью ДПФ нельзя назвать цифровой фильтрацией в полном смысле слова. Обычно цифровые фильтры, работающие в реальном масштабе времени, производят обработку сигнала непрерывно по мере его поступления, а вычисление выходного сигнала с помощью ДПФ может быть произведено лишь после того, как станет известным полностью входной сигнал или хотя бы первая серия из N отсчетов. Поэтому при использовании ДПФ выходной сигнал может быть получен только с некоторым запаздыванием по отношением к входному сигналу. Однако в ряде практических применений такое запаздывание выходного сигнала не играет существенной роли, и тогда обработка сигналов с использованием ДПФ оказывается целесообразной.

 

         3.3. Z-преобразование. Основные свойства Z-преобразования

 

         Как указывалось, методы описания непрерывных и дискретных сигналов во многом аналогичны друг другу. Обычному преобразованию Фурье соответствует дискретное преобразование Фурье; по аналогии с обычным преобразованием Лапласа можно ввести дискретное преобразование Лапласа.

         Запишем преобразование Лапласа для непрерывных сигналов

.                   (3.11)

         Это соотношение справедливо для сигналов, тождественно равных нулю при t<0. Поскольку все реально существующие сигналы имеют начало, преобразование Лапласа можно применять практически для всех сигналов.

         Применим преобразование Лапласа к дискретному сигналу, записанному в виде последовательности d-функций

;

.      (3.12)

         Выражение (3.12) представляет собой дискретное преобразование Лапласа. Практическое использование дискретного преобразования Лапласа затруднено тем, что изображения по Лапласу F(p) дискретных сигналов являются всегда трансцендентными функциями аргумента p (,  ,  и т.п.), в то время как для аналоговых сигналов, имеющих практическое приложение, изображения по Лапласу – как правило, дробно-рациональные функции, работать с которыми значительно легче.

         Применительно к дискретным сигналам вводится специальное преобразование, так называемое z-преобразование, играющее ту же роль, что и преобразование Лапласа, но без отмеченных трудностей в практическом использовании.

         Переход от дискретного преобразования Лапласа к    z-преобразованию осуществляется простой заменой переменной :

.        (3.13)

         В записи (3.13) вида связь каждого отсчета сигнала f(nT) с соответствующим ему членом z-преобразования F(z). Если в сигнале есть отсчет, расположенный в точке nT по времени, то в изображении должен быть член степени  с коэффициентом, равным величине этого отсчета, и наоборот. z-преобразование можно применять и к абстрактным числовым последовательностям.   Поскольку z-преобразование – это степенной ряд переменной , то важно рассмотреть вопрос о его сходимости. В математическом обосновании z-преобразования этот вопрос изучен. Показано, что    z-изображение существует как для убывающих, так и для неограниченно возрастающих последовательностей (за исключением ряда специфических, не представляющих интереса). Действия сложения, умножения, деления z-изображений законы без проверки на сходимость.        В качестве примеров рассмотрим z-преобразование простейших сигналов. При этом всюду будем полагать, что сигнал f(nT) тождественно равен нулю при n<0.

 

1 Единичный импульс (см. рисунок 3.1)

 

                                        F(z)=1.

 

 

 

 

 

 

 


                            Рисунок 3.1 Единичный импульс

 

2 Дискретизированный единичный скачок (см. рисунок 3.2)

 

 

 


                                                                                     

 

 

Рисунок 3.2 Дискретизированный единичный скачок

f(nT)=1, n³0;

 

.

 

3 Линейно-нарастающий дискретный сигнал (см. рисунок 3.3)

 

 

 

 

 

 

 

     

Рисунок 3.3 Линейно-нарастающий дискретный сигнал

 

f(nT)=n, n³0;

 

.

4 Экспоненциально убывающий дискретный сигнал

 

;

.

5 Комплексная экспонента

 

;

.

 

 

6 Гармоническая функция

;

.

 

7 Степенная функция

 

;

.

         Обратное z-преобразование позволяет определить значения дискретного сигнала f(nT) по заданному z-изображению F(z). Формула для обратного     z-преобразования имеет следующий вид:

,       (3.14)

где интегрирование выполняется против часовой стрелки по контуру С, внутри которого располагаются все особые точки функции F(z).

Для вычисления интеграла воспользуемся теоремой о вычетах:

,     (3.15)

где  - подынтегральная функция;

 - нуль знаменателя (полюс) функции L(z);

 - вычет функции L(z) в особой точке (т.е. в полюсе ), который находится согласно соотношению

.   (3.16)

В частности, вычет в простом полюсе  равен

,   (3.17)

а вычет в полюсе второй кратности равен

.      (3.18)

Пример 3.4.  Найдите дискретный сигнал f(nT) по его изображению .

         Для нашего случая

,

получим

.

         Подынтегральная функция имеет один простой полюс  и ее вычет (3.17) равен

.

         Окончательно .

Пример 3.5. Найдите дискретный сигнал f(nT) по его изображению .

         Для нашего случая

,

получим

.

         Подынтегральная функция имеет полюс второй кратности  и ее вычет (3.18) равен

.

Окончательно f(nT)=n.

К основным свойствам z-преобразования относятся свойства, связанные со сверткой последовательностей и задержкой последовательностей.

1 Свертка последовательностей

Пусть последовательности x(nT) и у(nT) имеют соответственно                  z-изображения X(z) и Y(z). Свертка последовательностей x(nT) и y(nT)

имеет z-изображение F(z), равное произведению X(z) и Y(z):

         F(z)= X(z)Y(z).                                                                           (3.19)

         Выражение (4.19) аналогично формуле (3.9).

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

         2 Задержка последовательностей

         Если последовательность x(nT) имеет z-изображение X(z), то последовательность y(nT)=x(nT-mT) имеет z-изображение Y(z), равное

.       (3.20)

         Таким образом, задержка последовательности во времени на mT соответствует умножению z-изображения на .

         Пример 3.6. Рассмотрим вычисление дискретной свертки f(nT) с помощью z-преобразования для дискретных сигналов, заданных в примере 3.3. Пусть конечные дискретные сигналы x(nT) и y(nT) с периодами  и  определяются в виде x(0)=2; x(T)=1; y(0)=1; y(T)=-1; y(2T)=2.

         Найдем z-изображения X(z), Y(z) и F(z) сигналов x(nT), y(nT) и f(nT):

;

;

.

         Поскольку в общем виде , можно записать: f(0)=2; f(T)=-1; f(2T)=3; f(3T)=2, что, естественно, совпадает с результатами примера 3.3.

Пример 3.7. Пусть x(nT)=1, n³0; и пусть y(nT)=x(nT-2T), т.е.

 

Так как , то

 

4 ГЛАВА. Характеристики цифровых фильтров

 

         4.1 Импульсная характеристика

 

         На вход цифрового фильтра подается входной сигнал x(kT) в виде последовательности числовых значений, следующих с интервалом T. При поступлении каждого очередного значения сигнала x(kT) в цифровом фильтре производится расчет очередного значения выходного сигнала y(kT). Алгоритмы расчета могут быть самыми разнообразными. Сигнал на выходе цифрового фильтра y(kT) также представляет собой последовательность числовых значений, следующих с интервалом T.

         Если на вход цифрового фильтра подать простейший сигнал в виде единичного импульса

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

                           

4.2 Алгоритм цифровой фильтрации

 

         Подадим на вход фильтра произвольный сигнал x(kT), представляющий собой набор дискретных значений x(0), x(T),... .

         Под действием первого элемента x(0) на выходе фильтра формируется последовательность g(kT), умноженная на x(0); при действии x(T) – последовательность g(kT), умноженная на x(T) и сдвинутая вправо на величину T, и т.д. В результате на выходе получим последовательность y(nT), причем

y(0)=g(0)x(0);

y(T)=g(T)x(0)+g(0)x(T);

y(2T)=g(2T)x(0)+g(T)x(T)+g(0)x(2T);

---------------------------------------

.   (4.1)

 

         Таким образом, выходной сигнал определяется как дискретная свертка входного сигнала и импульсной характеристики. Формула (4.1) представляет собой алгоритм цифровой фильтрации.

 

         4.3 Нерекурсивные и рекурсивные цифровые фильтры

 

Если импульсная характеристика фильтра описывается последовательностью с конечным числом членов, то фильтр может быть реализован в виде схемы, изображенной на рисунке 4.1.

 

 

 

 

 

 

 

 

 

Рисунок 4.1 Структурная схема нерекурсивного фильтра

        

Здесь буквой T обозначены элементы задержки сигнала на время T (на одну ячейку памяти); g(0), g(T),..., g(MT) – элементы, умножающие сигнал на соответствующий коэффициент. Структурная схема, изображенная на рис. 5.1, не является электрической схемой цифрового фильтра; эта схема представляет собой графическое изображение алгоритма цифровой фильтрации и показывает последовательность арифметических операций, выполняемых при обработке сигнала. Отметим, что иногда структурная схема фильтра может соответствовать его электрической схеме, т.е. реализации фильтра на конкретных микросхемах. Можно принять, что операция алгебраического сложения выполняется с помощью сумматора, умножения – с помощью множительного устройства и для задержки используются элементы памяти.      В большей мере, чем схемной реализации, структурная схема соответствует программной реализации фильтра, т.е. реализации путем составления соответствующей программы на ЭВМ. Степень этого соответствия зависит от особенностей конкретной ЭВМ.           В цифровых фильтрах элементы, задерживающие сигнал на одну ячейку памяти, обычно отмечают символом , обозначающим задержку сигнала на языке z-преобразования; а дискретные значения импульсной характеристики g(0), g(T),..., g(MT) обозначают коэффициентами .

         Вернемся к схеме цифрового фильтра, изображенной на рис. 4.1. Такие фильтры, где для расчета используют лишь значения входного сигнала, называют нерекурсивными (трансверсальными). Алгоритм нерекурсивного фильтра легко записать, если известна импульсная характеристика фильтра. Для практической реализации алгоритма необходимо, чтобы импульсная характеристика содержала конечное число членов. Если импульсная характеристика содержит бесконечное число членов, но они быстро убывают по величине, то можно ограничиться конечным числом членов, отбросив те, значения которых малы. В случае, если элементы импульсной характеристики не убывают по величине, алгоритм нерекурсивного фильтра оказывается нереализуемым.

         Таким образом, алгоритм цифровой фильтрации (4.1) для нерекурсивного фильтра можно записать в виде

.        (4.2)

         Для реализации бесконечных импульсных характеристик используют рекурсивные цифровые фильтры. Алгоритм цифровой фильтрации для рекурсивного фильтра имеет вид

. (4.3)

         Алгоритму (4.3) соответствует схема, изображенная на рисунке 4.2.

По сравнению с нерекурсивным фильтром здесь добавилась цепь обратной связи, которая означает, что значения выходного сигнала используются в последующих расчетах.         Рекурсивные фильтры позволяют производить обработку сигнала с более высокой точностью, так как они позволяют более правильно реализовать импульсную характеристику без отбрасывания ее “хвоста”. Рекурсивные фильтры могут реализовать алгоритмы, вообще нереализуемые с помощью нерекурсивных фильтров. Следовательно, нет смысла применять нерекурсивные алгоритмы для создания цифровых фильтров с импульсной характеристикой большой протяженности. В этом случае целесообразнее использовать рекурсивные фильтры.

Область применения нерекурсивных алгоритмов – это реализация цифровых фильтров с импульсной характеристикой, содержащей небольшое число членов.

 

         4.4 Системная (передаточная) функция

 

         Определим параметр цифрового фильтра, аналогичный передаточной функции H(p) электрической цепи. Для этого применим z-преобразование к импульсной характеристике g(kT) цифрового фильтра:

 

 

 


Рисунок 4.2 Структурная схема рекурсивного фильтра

 

.       (4.4)

Функцию H(z) называют системной функцией фильтра.

         В соответствии с выражением (4.1) сигнал на выходе цифрового фильтра равен дискретной свертке входного сигнала и импульсной характеристике фильтра. Применяя к этому выражению теорему о z-преобразовании свертки, получим, что z-преобразование выходного сигнала Y(z) равно z-преобразованию входного сигнала X(z), умноженному на системную функцию фильтра:

Y(z)=X(z)H(z).            (4.5)

Следовательно, системная функция

                                               H(z)=Y(z)/X(z)     (4.6)

 играет роль передаточной функции цифрового фильтра.

         Системная функция рекурсивного фильтра, описываемого алгоритмом цифровой фильтрации (4.3), записывается в следующем виде:

.        (4.7)

         Соотношение (4.7) получается в результате применения z-преобразования к выражению (4.3) и определению H(z) согласно (4.6). Системная функция нерекурсивного фильтра, описываемого алгоритмом цифровой фильтрации (4.2),

,          (4.8)

или

,         (4.9)

т.к.  - дискретные значения импульсной характеристики.

 

         4.5 Формы реализации цифровых фильтров

 

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

,    (4.10)

где

,     (4.11)

.         (4.12)

         Такому представлению соответствует каноническая схема фильтра, представленная на рисунке 4.3 и являющаяся каскадным соединениям фильтров  и .

В отличие от нее схему на рисунке 4.2 называют прямой формой цифрового фильтра. Преимуществом канонической схемы фильтра является минимальное число элементов задержки, равное порядку фильтра. Благодаря этому удается сократить общий объем вычислений.

Кроме прямой и канонической форм реализации цифровых фильтров существуют и другие варианты.

 

4.6 Критерий устойчивости

 

         Фильтр называется устойчивым, если все полюсы H(z) располагаются внутри единичной окружности в плоскости z. Нули могут быть расположены в любой точке плоскости z. Это означает, что для устойчивости фильтра необходимо и достаточно выполнение условия

, k=1,2,...,N,            (4.13)

где - полюс (корень знаменателя) функции H(z). Очевидно, что нерекурсивный фильтр всегда устойчив.

 

 

 

 

 

 

 

 

 

 

 

Рисунок 4.3 Каноническая схема фильтра

                  

4.7 Частотные характеристики цифровых фильтров

 

         Как известно, для описания аналоговых фильтров часто используют частотную характеристику

.        (4.14)

         Для получения частотных характеристик цифровых фильтров надо в H(z) вместо  подставить .

         Для рекурсивного фильтра согласно (4.7)

,            (4.15)

для нерекурсивного фильтра согласно (4.8) и (4.9)

,         (4.16)

или

.      (4.17)

Нетрудно видеть, что H(jw) является периодической функцией с периодом , т.к. , n – любое целое число.

         В общем случае H(jw) – комплексная функция, которая может быть записана в виде

.       (4.18)

         Для вычисления амплитудно-частотной характеристики (АЧХ) A(w) цифрового фильтра нужно взять модуль H(jw)

,     (4.19)

а для вычисления фазочастотной характеристики (ФЧХ)   j(w) – аргумент H(jw)

.     (4.20)

Для рекурсивного фильтра справедливы соотношения:

;           (4.21)

;         (4.22)

а для нерекурсивного фильтра – соотношения:

;    (4.23)

.                     (4.24)

 

Пример 4.1 На вход нерекурсивного цифрового фильтра с импульсной характеристикой g(kT) подается сигнал x(kT). Пусть ; ; ; x(0)=0,4; x(T)=1; x(2T)=0,7. Определить системную функцию фильтра, алгоритм цифровой фильтрации, сигнал на выходе фильтра.  Изобразить схему фильтра.

         Системную функцию фильтра определяем по формуле (4.8)

.

         Алгоритм цифровой фильтрации определяется по формуле (4.2)

.

         Это выражение определяет одновременно сигнал на выходе фильтра.

         Вычислим y(nT):

y(0)=0,6x(0)+x(-T)+0,3x(-2T)=0,24;

y(T)=0,6x(T)+x(0)+0,3x(-T)=1;

y(2T)=0,6x(2T)+x(T)+0,3x(0)=1,54;

y(3T)=0,6x(3T)+x(2T)+0,3x(T)=1;

y(4T)=0,6x(4T)+x(3T)+0,3x(2T)=0,21;

y(nT)=0 при n³5.

         Схема фильтра определяется алгоритмом цифровой фильтрации и дана на рисунок 4.4.

Пример 4.2. Найти частотную характеристику, АЧХ и ФЧХ нерекурсивного фильтра с системной функцией

.

 

 

 

 

 

 

 

 

 

Рисунок 5.4 Нерекурсивный цифровой фильтр

 

               

         Частотная характеристика фильтра определяется по формуле (4.16)

.

         АЧХ и ФЧХ фильтра определяется по формулам (4.23) и (4.24)

,

.

         Пример 4.3. На вход рекурсивного цифрового фильтра с системной функцией

подается сигнал x(kT). Пусть x(0)=1; x(T)=2; x(2T)=4. Записать алгоритм цифровой фильтрации, найти сигнал на выходе фильтра. Определить, устойчив ли фильтр.

         По определению (4.6), системная функция

,

или

.

         Используя теорему о запаздывании, запишем выражение алгоритма цифровой фильтрации

y(nT)-2y(nT-T)+5y(nT-2T)=2x(nT)-x(nT-T)+3x(nT-2T),

или

y(nT)=2x(nT)-x(nT-T)+3x(nT-2T)+2y(nT-T)-5y(nT-2T).

Вычислим y(nT) (n=0,1,...,6):

y(0)=2x(0)-x(-T)+3x(-2T)+2y(-T)-5y(-2T)=2;

y(T)=2x(T)-x(0)+3x(-T)+2y(0)-5y(-T)=7;

y(2T)=2x(2T)-x(T)+3x(0)+2y(T)-5y(0)=13;

y(3T)=2x(3T)-x(2T)+3x(T)+2y(2T)-5y(T)=-4;

y(4T)=2x(4T)-x(3T)+3x(2T)+2y(3T)-5y(2T)=-61;

y(5T)=2x(5T)-x(4T)+3x(3T)+2y(4T)-5y(3T)=-102;

y(6T)=2x(6T)-x(5T)+3x(4T)+2y(5T)-5y(4T)=101.

Полюсы находим, приравнивая к нулю знаменатель системной функции. Для фильтра второго порядка

.

Учитывая, что ; , получим

,

следовательно,

; .

Фильтр устойчив при  и .

.

Так как  и , то фильтр неустойчив.

 

         4.7  Разработка рекурсивных цифровых фильтров

 

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

Для алгебраического преобразования непрерывной передаточной функции в многочлен по z используется билинейное преобразование, известное в теории комплексных переменных под названием дробно-линейного преобразования.

Этапы разработки рекурсивных фильтров включают:

1) Задание частотной характеристики или передаточной функции фильтра.

2) Аппроксимация и расчет коэффициентов b(n) и a(m) передаточной функции фильтра.  Этот этап может выполняться четырьмя методами:

-   Метод размещения нулей и полюсов на комплексной z-плоскости.

-   Метод инвариантного преобразования импульсной характеристики.

-   Согласованное z-преобразование.

-   Билинейное z-преобразование.

3) Выбор структуры реализации фильтра – параллельная или каскадная, блоками второго и/или первого порядка.

4) Программное или аппаратное обеспечение реализации фильтра.

Метод размещения нулей и полюсов применяется при разработке простых фильтров с ограниченным количеством нулей и полюсов, если параметры фильтра не обязательно задавать точно. Амплитудная характеристики системы может быть оценена по выражениям при перемещении точки ws по единичной окружности exp(-jwsDt)::

|H(w)| = Ui /Vj.                                        (4.7.1)

Каждой точке zs = exp(-jwsDt) может быть поставлен в соответствие вектор (zsni) на ni -нуль, модуль которого Ui = |(zsni)| отображает расстояние от zs до i-нуля, а равно и вектор (zspj) на pj-полюс с соответствующим расстоянием Vj = (zspj). Наибольшее влияние на изменение АЧХ по частоте оказывают нули и полюсы, расположенные ближе к единичной окружности. При расположении нуля непосредственно на окружности гармоника  ws в этой точке полностью обнуляется (коэффициент передачи фильтра равен нулю). И, наоборот, при перемещении  ws к полюсу, близкому к единичной окружности, происходит резкое нарастание коэффициента усиления системы.

В качестве иллюстрации метода выполним расчет фильтра со следующими параметрами:

- Полная режекция сигнала на частотах 0 и 250 Гц.

- Полоса пропускания с центром на fp = 125 Гц с шириной полосы по уровню 3 дб Dp =10 Гц.

- Частота дискретизации данных fD = 500 Гц.

При частоте дискретизации 500 Гц интервал временной дискретизации Dt = 1/fD, а частота Найквиста fN = 1/2Dt = 250 Гц. Соответственно, нули передаточной функции располагаются в точках n1 = exp(-j2p 0 Dt) = 1 и n2 = exp(-j2p 250 Dt) = -1.. Угол из начала координат z-плоскости на полюс p с учетом его сопряженности для получения действительных коэффициентов ±180o . fp/fN = ±90o. Значение радиуса r до полюса определяет ширину полосы пропускания и в первом приближении (при 0 < r < 1.1) оценивается по выражению:

r » 1 + (Dp/fD)p.  r » 1.063.

         Передаточная функция:

H(z) = (z-1)(z+1) / [(z-r exp(jp/2) (z-r exp(-jp/2)] = (z2 -1)/(z2 +r2) = Y(z)/X(z).

z2Y(z)+r2Y(z) = z2X(z)-X(z).  Y(z) = [z2X(z)-X(z)-z2Y(z)] / r2 .

Алгоритм фильтра:

y(k) = [x(k-2) – x(k) – y(k-2)] / r2.

При использовании символики z-1 полюс располагается внутри единичной окружности на том же радиусе со значением (при r < 1):

r » 1 - (Dp/fD)p.  r » 0.937.

         Передаточная функция и алгоритм фильтра:

H(z) = (z-1)(z+1) / [(z-r exp(jp/2) (z-r exp(-jp/2)] = (z2 -1)/(z2 +r2) = (1-z-2) / (1+r2 z-2).

y(k) = x(k) – x(k-2) – r2 y(k-2).

         Характеристики фильтров приведены на рисунке 4.7.1. Индексы h1, H1, y1 относятся к первому фильтру с полюсом за пределами единичной окружности, индексами h2, H2, y2 – внутри окружности (символика z-1). Импульсные отклики фильтров получены подачей на их входы импульса Кронекера, частотные характеристики вычислены по импульсным откликам. Значение r первого фильтра подобрано по АЧХ под равный коэффициент усиления гармоники fp со вторым фильтром, после чего коэффициенты фильтров нормированы по коэффициенту усиления к 1 на частоте fp

Как следует из рисунков, изменение многочлена по степеням z на 1/z хотя и изменило коэффициенты разностного уравнения фильтра, но практически не повлияло на его амплитудно-частотную характеристику. Однако, при этом произошло изменение области сходимости фильтра с соответствующим изменением фазовых углов направления на полюсы из всех точек единичной окружности, что отразилось на фазово-частотной характеристике и отсчетах импульсного отклика фильтра изменением фазы на p. Синтез систем непосредственно в z-области применяется, в основном, только для режекторных и селекторных фильтров и более детально рассматривается ниже.

Метод инвариантного преобразования импульсной характеристики применяется для получения из подходящей аналоговой передаточной функции H(s) с помощью преобразования Лапласа импульсной характеристики h(t), которая затем дискретизируется и подвергается z-преобразованию.

         Допустим, имеется простая аналоговая система с передаточной функцией:

H(s) = C / (s-p).

         Выполняем обратное преобразование Лапласа функции H(s) и дискретизируем результат преобразования с определенной постоянной времени Dt:

h(t) = TL-1[H(s)] = C exp(pt) → C exp(pnDt).

 

Выполняем z-преобразование и формируем передаточную функцию H(z):

H(z) =h(nDt) zn =C exp(pnDt) zn = C / (1-z exp(pDt)).

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

Согласованное z-преобразование применяется для преобразования аналоговых фильтров в эквивалентные цифровые непосредственным переводом всех полюсов и нулей с s-плоскости в z-плоскость:

 

9.2.1.GIF

 

Рисунок 4.7.1.- Инвариантное преобразование импульсной характеристики

 

(s-a) → z exp(aDt).

         Большинство полюсов и нулей являются комплексно сопряженными и реализуются фильтрами второго порядка:

(s-a)(s-a*) → 1 – 2z exp(Re(a)Dt) cos(Im(a)Dt) + z2 exp(Re(a)Dt).

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

Билинейное z-преобразование является основным методом получения коэффициентов рекурсивных БИХ-фильтров и использует следующую замену:

s = g (1-z) / (1+z),   g = 1 или 2/Dt,

при этом ось jw s-плоскости отображается в единичную окружность z-плоскости, правая половина s-плоскости – внутрь единичной окружности, а левая половина с полюсами устойчивых аналоговых фильтров – снаружи единичной окружности. Аналогичная замена при отрицательной символике z-1 с соответствующей сменой отображения:

s = g (z-1) / (z+1).

         Для фильтров верхних и нижних частот порядок фильтра H(z) равен порядку фильтра H(s). Для полосовых и заградительных фильтров порядок H(z) вдвое больше порядка H(s). Для сохранения частотных характеристик фильтра при нелинейном сжатии частотной шкалы аналогового фильтра (переход от ∞ к wN) предварительно выполняется деформация частотной шкалы аналогового фильтра. Более подробно эти вопросы рассмотрены ниже.

 

4.8. Режекторные и селекторные фильтры

Режекторный фильтр (фильтр-пробка) подавляет определенную частоту во входном сигнале. Он может быть спроектирован непосредственно по z-диаграмме.

Комплексная z-плоскость. Простейший режекторный фильтр типа НЦФ имеет один нуль на единичной окружности в z-плоскости в точке с частотой, которую необходимо подавить. Так, например, если из входного сигнала требуется исключить постоянную составляющую (нулевая частота), то импульсная реакция фильтра НЦФ имеет вид:

H(z) = 1-z.                                                  (4.8.1)

Нуль функции (4.8.1) равен n1=1. Как можно видеть на рисунке 4.8.1, коэффициент передачи сигнала H(w) на любой частоте wi от 0 до wN=p/Dt - частоты Найквиста, определяемый выражением (4.8.1), будет равен длине вектора Vn1, проведенного из нуля функции H(z) - точка n1 на действительной оси, до соответствующей частоты wi - точки z(wi) на единичной окружности. На частоте w = 0 длина этого вектора равна нулю. Амплитудно-частотная характеристика фильтра, приведенная пунктиром на рисунке 4.8.2 для передаточной функции (4.8.1), далека от идеальной для фильтр-пробки.

Режекторный фильтр постоянной составляющей сигнала.

Сконструируем простейший РЦФ, добавив к оператору (4.3.1) один полюс вне единичной окружности на малом расстоянии от нуля:

Hп(z) = G(1-z)/(1-az),  zp= 1/a.                                  (4.8.2)

Допустим, что полюс помещен в точке p1 = 1.01, при этом, а=0,99. Масштабный коэффициент G получим нормировкой H(z) к 1 на частоте Найквиста. Для приведенных условий G=0.995. Отсюда,  при t=1:

Hп(z) = 0,995(1-z)/(1-0.99z),

y(k)  = 0.995[x(k) –x(k-1)]+ 0.99y(k-1).

Отображение нуля n1 и полюса р1 на z-плоскости и АЧХ фильтра для исключения постоянной составляющей приведены на рисунке 4.8.1. Коэффициент передачи сигнала на произвольной частоте wi равен отношению длин векторов Vn1(z) и Vp1(z) соответственно из нуля и полюса до точки z(wi) на единичной окружности и близок к единице для всех частот, за исключением нулевой:

|Hп(z)| = G Vn1(z)/Vp1(z).

Режекторный фильтр произвольной частоты. При проектировании на подавление любой другой частоты частоты wv нули и полюсы располагаются на соответствующем радиусе z-плоскости. Радиальный угол направления на нуль и полюс определяются выражением:

jv = p·wv/wN.                                            (4.8.3)

Наличие двух знаков в выражении (4.8.3) отражает тот факт, что для получения вещественной функции фильтра нули и полюсы должны быть комплексно-сопряженными парами (их произведение дает вещественную функцию), т.е.:   Hv(z) = G(z-zn)(z-zn*)/[(z-zp)(z-zp*)].                           (4.8.4)

  

 

 

Рисунок 4.8.1. - Синтез фильтров.      Рисунок 4.8.2.- АЧХ  фильтра

         

Нули фильтра располагаются на единичной окружности:

zn = cos jv + j sin jv = Re zn + j Im zn.                          (4.8.5)

Полюсы - на полярном радиусе R:

zp = R cos jv + j R sin jv = Re zp + j Im zp.                       (4.8.6)

Пример положения нулей (n2 и n2*) и полюсов (р2 и р2*) приведен на рисунке 4.8.1. Подставляя (4.8.5-4.8.6) в (4.8.4), получаем:

Hv(z) =,                                   (4.8.7)

G = [1+(1+2Re zp)/R2] / (2+2Re zn).                               (4.8.8)

         При приведении уравнения (4.8.7) в типовую форму:

Hv(z) =,                                      (4.8.7)

b0 = 1,    b1 = -2·Re zn,    b2 = 1.                                  (4.8.9)

a1 = - (2·Re zp)/R2,     a2 = 1/R2.

         Соответственно, алгоритм вычислений:

y(k)  = G·[x(k) +b1·x(k-1)+x(k-2)] – a1·y(k-1) – a2·y(k-2).           (4.8.10)

         Для проверки вычисленного в примере фильтра на рисунке 4.8.4 приведен модельный входной сигнал, состоящий из суммы двух равных по амплитуде гармоник с частотой 50 и 53 Гц, и сигнал на выходе фильтра (смещен вверх). Справа на рисунке приведены спектры входного и выходного сигналов. Спектр выходного сигнала зарегистрирован после интервала установления реакции фильтра, который хорошо заметен на начальной части графика выходного сигнала. После установления сигнал на выходе фильтра практически полностью освобожден от гармоники 50 Гц.

При R → 1 ширина полосы подавления фильтра становится все более узкой, но при этом увеличивается длительность импульсной реакции фильтра и, соответственно, время установления фильтра при изменении спектра входного сигнала.

D05-02

 

Рисунок 4.8.4. - Спектры входного и выходного сигналов.

 

В первом приближении значимая часть импульсной реакции режекторных фильтров равна (4÷5)/(R-1). Отклик фильтра получен при подаче на вход РЦФ импульса Кронекера.

Селекторный фильтр. Если в уравнении (4.8.4) опустить нули, то получим селекторный фильтр, выделяющий сигналы одной частоты ωs – частоты селекции, с передаточной функцией:

Hs(z) = G/[(z-zp)(z-zp*)] = G1/(1+a1z+a2z2).                     (4.8.11)

         Характер передаточной функции (4.8.11) можно представить непосредственно по z-плоскости (см. рисунок 4.8.1). При расположении полюсов фильтра за пределами единичного круга (например, в точках р2 и р2*) значение коэффициента передачи фильтра на произвольной частоте ω на единичной окружности будет обратно пропорционально величине векторов из этих точек окружности на полюса фильтра.

При изменении ω от нуля до  ±π (движение по единичной окружности на z-плоскости по или против часовой стрелки) один из векторов (на полюс противоположной полуплоскости) изменяется в достаточно небольших пределах (не превышая значения 2), в то время как второй из векторов (на полюс в своей полуплоскости) будет сначала уменьшаться, достигает минимума при расположении ω на полярном радиусе полюса (на частоте селекции ωs), а затем снова начинает увеличиваться.

Соответственно, значение Hs(ω) максимально на частоте селекции  ±ωs и при R → 1 может быть очень высоким.

При необходимости фильтр может быть пронормирован к 1 на частоте селекции определением значения G1 по условию Hs(ω) = 1 при ω = ωs, т.е.:

G1 = 1+a1 z(ws)+a2 z(ws)2.

Фильтр в принципе не может иметь нулевого коэффициента передачи на других частотах главного диапазона. Если последнее является обязательным, то фильтр выполняется методом обращения режекторного фильтра Hv(z):

 

Hs(z) = 1-Hv(z).

Hs(z) = .      (4.8.12)

с0 = 1-G,   c1 = a1-Gb1,    c2 = a2-G.

 

D05-10

 

 

 

 

 

 

 

 

 

 

Рисунок 4.8.7.- График передаточной функции

  

Пример передаточной функции фильтра приведен на рисунке 4.8.7.           Пример применения фильтра для выделения гармонического сигнала на уровне шумов, мощность которых   больше мощности сигнала, приведен на рис. 4.8.8.

 

D05-04

 

 

 

 

 

 

 

 

 

Рисунок 4.8.8. - Фильтрация сигнала селекторным РЦФ.

 

 

5 ГЛАВА. Основы синтеза цифровых фильтров

 

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

         Задача синтеза цифровых фильтров сводится к нахождению коэффициентов (обозначаемых через ak и bk) его системной функции, обеспечивающих аппроксимацию заданных характеристик фильтра. Таким образом, задача синтеза фильтра сводится фактически к задаче аппроксимации. Задача реализации не представляет никакого труда, так как

по определенной системной функции фильтра  можно     сразу же изобразить

схему фильтра.

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

 

5.1. Метод инвариантной импульсной характеристики

 

         Наиболее просто задача синтеза цифрового фильтра решается в том случае, если известна импульсная характеристика фильтра-прототипа. Метод синтеза цифровых фильтров, основанный на использовании импульсной характеристики фильтра-прототипа, называют методом инвариантной импульсной характеристики. Его называют также методом стандартного z-преобразования. Согласно этому методу, для определения импульсной характеристики проектируемого цифрового фильтра необходимо подвергнуть дискретизации импульсную характеристику g(t) аналогового фильтра-прототипа. Значения импульсной характеристики цифрового фильтра должны быть равны значениям импульсной характеристики фильтра-прототипа в равномерно распределенные отсчетные моменты времени t=kT. Применяя к импульсной характеристике цифрового фильтра стандартное z-преобразование, можно найти системную функцию и составить алгоритм цифровой фильтрации.

         Пример 5.1. Рассмотрим синтез фильтра второго порядка, эквивалентного колебательному контуру, методом инвариантной импульсной характеристики.

         Импульсная характеристика колебательного контура имеет вид

,

где A – размерный множитель.

         Заменяя t на kT и опуская множитель A, запишем импульсную характеристику соответствующего цифрового фильтра

.

         Найдем системную функцию этого фильтра как z-преобразование от импульсной характеристики:

        

.

Таким образом,    .

         Алгоритм цифровой фильтрации данного фильтра запишем в виде

         .

         Схема цифрового фильтра, соответствующая этому алгоритму, приведена на рисунке 5.1.

Рассмотрим достоинства и недостатки метода инвариантной импульсной характеристики.

В результате дискретизации импульсной характеристики частотная характеристика цифрового фильтра представляет собой сумму сдвигов

частотной характеристики исходного аналогового фильтра.

 

 

 

 

 

 

 

 

 

 

  

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

         Отсюда следует, что недостатком этого метода является наличие эффекта наложения. Методы борьбы с этим недостатком состоят в увеличении частоты дискретизации и ограничении полосы пропускания аналогового фильтра-прототипа дополнительным фильтром. При использовании этих методов эффекты наложения могут оказаться пренебрежимо малыми и частотные характеристики аналогового и цифрового фильтров станут похожими друг на друга.

         Достоинства метода следующие: во-первых, точное совпадение импульсных характеристик цифрового фильтра и аналогового  фильтра- прототипа в отсчетные моменты времени t=kT; во-вторых, относительная простота расчета.

                

5.2. Метод билинейного z-преобразования

 

         Синтез цифровых фильтров по заданной частотной характеристике фильтра-прототипа более сложен. При этом используются различные методы, основным из которых является метод билинейного z-преобразования.

         Метод билинейного z-преобразования частотной характеристики заключается в пересчете известной передаточной функции H(p) аналогового фильтра-прототипа в системную (передаточную) функцию H(z) цифрового фильтра посредством билинейного z-преобразования.

         Билинейное z-преобразование комплексной переменной p в комплексную переменную z имеет следующий вид:

.            (5.2.1)

         Шкала соответствия между “аналоговыми” частотами  и “цифровыми” частотами , которую можно получить из (5.2.1), задается следующим соотношением:

,         (5.2.2)

т.е. вся ось “аналоговых” частот  преобразуется в совокупность интервалов “цифровых” частот вида , где n – любое целое число. Частотная характеристика аналогового фильтра-прототипа повторяется во всех частотных интервалах без эффекта наложения.

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

         Совокупность характерных частот цифрового фильтра известна. Использую соотношение (5.2.2), пересчитаем все характерные частоты цифрового фильтра в характерные частоты аналогового фильтра-прототипа. Затем по заданным требованиям выбирают вид аппроксимации идеальной характеристики и порядок фильтра-прототипа. Далее по справочникам определяют передаточную функцию H(p) аналогового фильтра-прототипа. Окончательно определяют системную функцию H(z) цифрового фильтра с помощью билинейного z-преобразования (5.2.1) и форму его реализации.

         Пример 5.2. Требуется рассчитать цифровой фильтр нижних частот, имеющий плоскую частотную характеристику в области нижних частот с затуханием 3 дБ на частоте среза fс = f1 = 1кГц и ослаблением не менее A = 20 дБ на частоте f2 = 2 кГц. Частота дискретизации fд = 10 кГц.

         Из теории аналоговых фильтров известно, что при соответствующем выборе параметров заданным требованиям могут удовлетворять фильтры Баттерворта, квадрат модуля частотной характеристики которых описывается выражением

.      (5.2.3)

         Прежде всего найдем частотную характеристику аналогового фильтра-прототипа. Между характерными частотами цифрового фильтра и фильтра-прототипа существует однозначная связь, которая устанавливается соотношением (5.2.3). В соответствии с техническим заданием частотная характеристика проектируемого цифрового фильтра имеет две характерные частоты: частоту среза , соответствующую затуханию 3 дБ, и частоту , на которой должно обеспечиваться затухание не менее 20 дБ. Частоте среза цифрового фильтра  соответствует частота  среза аналогового фильтра-прототипа:

,

где с – интервал дискретизации.

         Подставляя числовые значения, получим  1/c. Аналогично частоте  цифрового фильтра соответствует частота  аналогового фильтра:

 1/с.

         Таким образом, аналоговый фильтр-прототип должен иметь затухание 3дБ на частоте 1/с и затухание не менее 20 дБ на частоте 1/с. Отношение этих двух частот . Из выражения (5.2.3) определяем порядок фильтра-прототипа:

,

или

,

откуда

.

         Подставляя числовые значения, получим

.

Выбираем ближайшее большее целое значение n=3.

         Найдем передаточную функцию фильтра-прототипа. Операторный коэффициент передачи фильтра Баттерворта третьего порядка определяется выражением

 ,        (5.2.4)

где  1/с – частота среза аналогового фильтра-прототипа.

         Заменяя p на , после подстановки числовых значений и ряда упрощений получим выражение для системной функции цифрового фильтра

.                  (5.2.5)

Эта системная функция, представленная в виде произведения двух системных функций первого и второго порядков, может быть реализована, например, в виде канонической схемы фильтра (см. рисунок 5.2.1).

 

 

 


Метод, аналогичный рассмотренному, используется не только для синтеза фильтров нижних частот, но и для синтеза фильтров верхних частот, полосовых и режекторных. Для каждого типа фильтра получена формула билинейного z-преобразования, которая позволяет от аналогового фильтра-прототипа нижних частот перейти к синтезируемому цифровому фильтру.

Метод билинейного z-преобразования обладает следующими достоинствами:

- во-первых, физически реализуемый и устойчивый аналоговый фильтр отображается в физически реализуемый и устойчивый цифровой фильтр;

- во-вторых, отсутствуют проблемы, связанные с эффектом наложения, который характерен для метода инвариантного преобразования импульсной характеристики;                

-в-третьих, искажение частотного масштаба цифрового фильтра, преобразованного из аналогового фильтра-прототипа, можно учесть для широкого класса фильтров. Для этого достаточно пересчитать все характерные частоты аналогового фильтра-прототипа так, чтобы после выполнения билинейного z-преобразования (обратного пересчета) характерные частоты цифрового фильтра совпали бы с заданными.

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

Тем не менее метод билинейногоz-преобразования является самым распространенным аналитическим методом расчета цифровых фильтров.

 

         5.3. Метод частотной выборки

 

         До сих пор мы рассматривали методы синтеза цифровых фильтров с импульсной характеристикой неограниченной длины – так называемые фильтры с бесконечной импульсной характеристикой, или БИХ-фильтры. Однако наряду с БИХ-фильтрами широкое распространение получили цифровые фильтры с импульсной характеристикой конечной длины, или   КИХ-фильтры. Фильтры с конечной импульсной характеристикой реализуются, как правило, по нерекурсивной схеме и обладает рядом положительных качеств:          во-первых, легко создать КИХ-фильтры со строго линейной фазовой характеристикой. Действительно, для этого достаточно, чтобы импульсная характеристика цифрового фильтра g(kT) была симметричной. Этого нельзя достигнуть в БИХ-фильтрах, так как      из-за бесконечной протяженности импульсной характеристики g(kT) для k>0 и g(kT)=0 для k<0 эта характеристика принципиально не может быть строго симметричной; во-вторых, КИХ-фильтры, благодаря отсутствию обратных связей, всегда устойчивы.

         Фильтры с конечной импульсной характеристикой не имеют непосредственных аналогов среди пассивных электрических фильтров, поэтому методы их синтеза относятся к прямым методам. Одним из наиболее эффективных методов синтеза КИХ-фильтров является метод частотной выборки.       Метод частотной выборки формулируется полностью в частотной области и поэтому весьма удобен для проектирования фильтров по заданной частотной характеристике. Идея метода очень проста. Заданную частотную характеристику подвергают дискретизации как периодическую функцию с периодом, равным , разбивая интервал   на N равных частей. Интервал дискретизации по частоте  должен быть таким, чтобы передать характерные подробности частотной характеристики. К  образовавшейся последовательности отсчетов частотной характеристики  применяют обратное дискретное преобразование Фурье и находят импульсную характеристику. Значения g(kT) являются коэффициентами КИХ-фильтра.     Иногда получающаяся в процессе синтеза импульсная характеристика цифрового фильтра оказывается физически нереализуемой, так как отличны от нуля ее значения g(kT) при отрицательных kT. Фильтр можно сделать физически реализуемым, если его импульсную характеристику сместить вправо так, чтобы для моментов времени kT<0 она тождественно равнялась нулю.

         В качестве иллюстрации этого метода рассмотрим синтез фильтра нижних частот с идеальной прямоугольной частотной характеристикой. На рисунке 5.3.1, а изображена заданная частотная характеристика (пунктирная линия) и частотные отсчеты (точки). Частотная характеристика синтезированного фильтра (см. - сплошная линия на рисунке 5.3.1,а) точно пройдет через все заданные значения в точках отсчетов, но в промежутках между ними она может недопустимо отклоняться от идеальной характеристики. Это приведет к пульсации частотной характеристики в полосах пропускания и непропускания. При увеличении числа отсчетов увеличивается частота пульсаций частотной характеристики, а амплитуды пульсаций уменьшаются. Однако амплитуды последней пульсации в полосе пропускания и первой пульсации  в полосе непропускания остаются практически неизменными.

         Для уменьшения уровня пульсаций отсчетные значения частотной характеристики задают только в полосе пропускания фильтра, т.е. при  (см. рисунок 5.3.1, б), и полосе непропускания (), а в  переходной полосе () их подбирают так, чтобы обеспечить минимальный уровень пульсаций. Эта задача решается методами оптимизации на ЭВМ. Наличие отсчетов в переходной полосе приводит, естественно, к ее расширению. Этим приходится расплачиваться за уменьшение уровня пульсаций.

Еще один метод синтеза цифровых фильтров с конечной импульсной характеристикой основан на использовании так называемых временных окон.

        

5.4. Метод взвешивания с помощью окна

 

Рассмотрим идею этого метода на примере фильтра нижних частот.

Идеальный фильтр нижних частот с частотой среза  имеет отсчеты импульсной характеристики вида

,

где k – любое целое число.

         Наиболее простой подход к расчету КИХ-фильтров состоит в аппроксимации импульсной характеристики бесконечной длины идеального фильтра импульсной характеристикой конечной длины. Непосредственное усечение характеристики бесконечной длины приведет к пульсации частотной характеристики. Частотная характеристика будет иметь вид типичный для фильтров, синтезируемых методом частотной выборки (см. рисунок 5.4, а).

Уровень пульсаций частотной характеристики может быть уменьшен путем использования менее резкого усечения импульсной характеристики. Импульсную характеристику g(kT) умножают на весовую функцию (временное окно) W(kT), которая близка к единице в середине (kT=0) и плавно убывает к краям. В результате форма импульсной характеристики становится более плавной и, как следствие этого, улучшается вид частотной характеристики, которая также будет иметь вид типичный для фильтров, синтезируемых методом частотной выборки (см. рисунок 5.3, б).

 

 

 

 

 

 

 

 

 

 

 

Рисунок 5.4 - Частотная характеристика фильтра

        

         Метод умножения импульсной характеристики на функцию окна называется методом взвешивания с помощью окна.

         Существует несколько различных типов временных окон. Одним из наиболее распространенных является окно Хемминга, которое описывается функцией

, .

Окно Хемминга позволяет получить частотную характеристику с приемлемым уровнем пульсаций.

Достоинствами методов синтеза фильтров с конечной импульсной характеристикой (частотной выборки и временных окон) является их относительная простота, большой набор рассчитанных временных окон, а недостатками – негибкость (трудность определения частоты среза частотной характеристики), размывания частотной характеристики функцией окна или отсчетами в переходной полосе.

 

6 ГЛАВА.  Деконволюция цифрвых сигналов

 

Основное назначение деконволюции (deconvolution) – восстановление истинной формы сигнала, несущего информацию об исследуемом физическом или технологическом процессе, явлении природы и т.п., после его искажения при регистрации какой-либо линейной системой - измерительным трактом прибора (аппаратной или приборной функцией) или каналом связи. Естественно, что для восстановления необходимы сведения о характеристиках искажающей системы, и в первую очередь, об импульсном отклике системы или его частотной передаточной функции. Для выполнения деконволюции реализуются фильтры, частотные характеристики которого обратны частотной характеристике искажающих систем. Построение таких фильтров не всегда возможно. Так, невозможно в принципе восстановить в сигнале частоты, которые были полностью подавлены, а при восстановлении частотных составляющих, ослабленных до уровня шумов, одновременно происходит значительное усиление дисперсии шумов, в которых полезный сигнал может полностью затеряться.

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

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

 

6.1. Понятие деконволюции

 

При достаточно сложном физическом представлении во временной (координатной) области деконволюция проста для понимания в частотном представлении. Допустим, что в какой-либо регистрирующей системе происходит резонансное поглощение энергии и сдвиг по фазе определенного гармонического колебания в составе входного сигнала, например, преобразование гармоники sin w0t ® 0.5 sin (w0-p/4. Соответственно, для восстановления первоначальной формы сигнала операция деконволюции должна выполнить усиление данной гармоники в выходном сигнале в 2 раза и осуществить обратный сдвиг фазы на p/4. Для одной гармоники выполнение такой операции труда не представляет. Но практические задачи деконволюции значительно сложнее, так как требуется, как правило, восстановление полного спектра исходного сигнала, имеющего непрерывный характер. 

Определение деконволюции. Если для прямой свертки дискретного сигнала x(k) c импульсным откликом h(n) линейной системы (фильтра) имеем уравнение:

y(k) = h(n) x(k) ó H(z)X(z) = Y(z),   Y(z) =yk zkz = exp(-jj),

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

X(z) = Y(z)/H(z) = Y(z)H-1(z) ó y(k) h-1(n) = x(k),                 (6.1.1)

где индексом "-1" символически обозначена передаточная функция оператора обратного фильтра H-1(z) = 1/H(z), т.е. оператора деконволюции, инверсного прямому оператору свертки (импульсному отклику системы). Соответственно, при обратном z-преобразовании получаем оператор деконволюции:

H-1(z) = 1/H(z) ó h-1(n).                                        (6.1.2)

Очевидно, что если имеет место H(z)H-1(z) = 1, то при обратном z-преобразовании этого выражения мы должны иметь:

H(z)H-1(z) = 1 ó h(n) h-1(n) = do(n),                             (6.1.3)

где do(n) - импульс Кронекера. При этом последовательная свертка сигнала x(k) с оператором системы h(k) и оператором деконволюции h-1(k) будет эквивалентна свертке сигнала x(k) с импульсом Кронекера, что не должно изменить форму сигнала x(t).

При z = exp(-jj) все изложенное действительно и для спектрального представления операторов. Пример инверсии оператора через спектральное представление приведен на рисунке 6.1.1 (исходный оператор hn ® спектральная плотность H(ω) ® инверсная спектральная плотность H-1(ω) ® инверсный оператор h-1n на начальном интервале отсчетов).

Особенности деконволюции. Выражение (6.1.2) позволяет сделать некоторые выводы об особенностях выполнения деконволюции.

При ограниченной импульсной реакции h(n) инверсный оператор h-1(n) в общем случае не ограничен. Так, например, если импульсная реакция представлена нормированным диполем h(n) = {1,a} ó (1+az) = h(z), то имеем:

H-1(z) = 1/(1+az) = 1-az+a2z2-a3z3+ ....

h-1(n) = {1, -a, a2, -a3,....}.

Это действительно практически для любых операторов свертки, энергия которых на каких-либо ограниченных участках главного частотного диапазона близка к нулевой.

 

D07-01

 

Рисунок 6.1.1.- Инверсия оператора через спектральное представление

 

При инверсии спектральной функции таких операторов на этих участках возникают резкие спектральные пики, которые при обратном преобразовании Фурье дают медленно затухающие функции операторов. Пример такого явления приведен на рисунке 6.1.2.

Отсюда следует, что для точного выполнения деконволюции необходимо располагать бесконечно длинным инверсным оператором фильтра. Практически, деконволюция выполняется, если инверсный оператор достаточно быстро затухает и может быть ограничен. Но использование усеченных операторов приводит к появлению определенной погрешности деконволюции, величину которой следует контролировать.

 

D07-02

 

 

 

 

 

 

 

 

Рисунок 6.1.2. - Пример спектральных пиков

 

Заметим также, что передаточные функции систем, как правило, имеют низкочастотный характер. Инверсия операторов таких систем всегда связана с усилением высоких частот, что приводит к высоким коэффициентам усиления инверсными фильтрами дисперсии помех, что может приводить к потере полезного сигнала среди усиленных шумовых флюктуаций.

Кроме того, быстрое затухание оператора деконволюции является необходимым, но еще не достаточным условием реализации деконволюции.

Устойчивость фильтров деконволюции. Функция H(z) в выражении (5.2) имеет особые точки - нули функции, которые становятся полюсами функции H-1(z) = 1/H(z) и определяют устойчивость инверсного фильтра. Для того чтобы фильтр деконволюции был устойчивым, ряд 1/H(z) должен сходиться, т.е. полюса функции должны находиться вне единичного круга на z-плоскости (внутри круга при использовании символики z-1).     

Многочлен H(z) порядка N может быть разложен на N простых сомножителей - двучленов (диполей):

H(z) = (a-z)(b-z)(c-z)....,                                       (6.1.4)

где а, b, с,… - корни полинома. Обращение передаточной функции:

H-1(z) =                                      (6.1.5)

Если каждый из диполей функции (13.1.4) является минимально-фазовым диракоидом, т.е. корни диполей находится вне единичного круга на z-плоскости и модули нулевых членов диполей всегда больше следующих за ними первых членов (в данном случае: |а|>1, |b|>1, |с|>1), то функция H(z) также является минимально-фазовым диракоидом. Максимум энергии импульсного отклика минимально-фазового диракоида сосредоточен в его начальной части и последовательность отсчетов представляет собой затухающий ряд. При этом функция 1/H(z) также будет представлять собой передаточную функцию минимально-фазового оператора, обеспечивающую выполнение условия (6.1.3), а ее обратное преобразование - сходящийся ряд устойчивого фильтра. Так, например, фильтр, реализующий передаточную функцию (6.1.5), в самой общей форме может быть выполнен в виде включенных последовательно фильтров, каждый из которых имеет передаточную функцию следующего типа (для первого фильтра):

H1-1(z) = 1/(a-z) = (1+z/a+z2/a2+...)/a.

         Отсюда, равно как и непосредственно из выражения (6.1.5), следует также, что чем больше значения модулей корней фильтра (дальше от единичного круга полюс фильтра), тем быстрее затухает инверсный оператор.

Обращение недиракоидных функций. Если H(z) - реверсоид, т.е. корни составляющих его диполей находятся внутри и на единичном круге в z-плоскости, то устойчивое обращение H(z) является антиимпульсом (с отрицательными степенями z) и для его использования необходимо располагать "будущими" значениями входного сигнала.

Если диполи функции (6.1.4) представляют собой и диракоиды, и реверсоиды, то обращение будет центроидом и определяется полным рядом Лорана:

H-1(z) = ...+h-2z-2+h-1z-1+h0+ h1z1+h2z2+ ...,

т.е. оператор инверсного фильтра является двусторонним и не обязательно симметричным, как мы обычно рассматривали ранее двусторонние операторы.

 

6.2. Инверсия импульсного отклика фильтра

 

Вычисление коэффициентов инверсного фильтра по значениям каузального (одностороннего) оператора h(n) может быть проведено на основе выражения (6.2.3):

h-1(k)h(n-k) = do(n),                                       (6.2.1)

для чего достаточно развернуть его в систему n-уравнений при n = 0, 1, 2, … ,  k  ≤ n:

n = 0:   h-1(0)h(0) = 1,                                      h-1(0) = 1/h(0).

n = 1:   h-1(0)h(1)+h-1(1)h(0) = 0,                    h-1(1) = h-1(0)h(1) / h(0).

n = 2:h-1(0)h(2)+h-1(1)h(1)+h-1(2)h(0)=0,      h-1(2) = (h-1(0)h(2)+h-1(1)h(1))/h(0),   и т.д.

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

h-1(n) = (1/h0)h-1(k)h(n-k).                                (6.2.2)

Если фильтр деконволюции устойчив и ряд h-1(n) сходится, то появляется возможность ограничения количества членов ряда с определенной ошибкой восстановления исходного сигнала. Метрика приближения Е (квадратичная норма разности) определяется выражением:

Е2 =[do(n) - h(n) h-1(n)]2.                                    (6.2.3)

         Ошибка восстановления исходного сигнала появляется со сдвигом на длину оператора фильтра деконволюции и проявляется на интервале длины прямого оператора.

 

6.3. Оптимальные фильтры деконволюции

 

Можно рассчитать оптимальные фильтры деконволюции, метрика приближения которых меньше, чем у усеченных фильтров деконволюции. Для получения общего уравнения оптимальной деконволюции будем считать, что число коэффициентов оператора hn равно M+1, a число коэффициентов инверсного оператора hn-1 равно N+1.

Принцип оптимизации. Выходная функция приближения при использовании уравнения свертки (6.1.2) с ограничением числа членов оператора фильтра:

F = Е2 = [do(k)-sk]2.      sk =hn-1 hk-n.    (6.3.1)

Чтобы определить минимум функции, приравняем нулю  частные производные от Е по неизвестным коэффициентам фильтра:

dF/dhj-1 = -2hk-j [do(k) -hn-1 hk-n] = 0.                    (6.3.2)

hk-j hn-1 hk-n = hk-j do(k) = h-j.                       (6.3.3)

hn-1 hk-n hk-j  = hn-1 aj-n = h-j,    j = 0,1,2, ..., N,            (6.3.4)

где aj-n - функция автоковариации импульсной реакции h(n). Учитывая также, что hn = 0 при n<0 и aj = a-j (функция автоковариации является четной функцией), окончательное решение определяется следующей системой линейных уравнений:

  a0 h0-1 +      a1 h1-1 +        a2 h2-1 +        a3 h3-1 + ...    + aN hN-1          = h0                (6.3.5)

 a1 h0-1 +       a0 h1-1 +        a1 h2-1 +        a2 h3-1 + ...    + aN-1 hN-1       = 0

 a2 h0-1 +       a1 h1-1 +        a0 h2-1 +        a1 h3-1 + ...    + aN-2 hN-1       = 0

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

 aN h0-1 +      aN-1 h1-1 +     aN-2 h2-1 +     aN-3 h3-1 +...  +a0 hN-1             = 0

 

Как видно на рисунке 6.3.1, оптимизация инверсного оператора заключается в центрировании ошибок приближения и с распределением по интервалу суммарной длины прямого и инверсного оператора.

Уравнение оптимальной инверсии. Оптимальный инверсный фильтр может быть получен непосредственно с использованием z-образов импульсной реакции и автоковариационной функции прямого фильтра. Если для прямого фильтра мы имеем передаточную функцию H(z), то z-образ автоковариационной функции фильтра (как z-отображение спектральной плотности мощности) представляет собой произведение:

A(z) = H(z)H*(z),                                            (6.3.6)

где H*(z)- функция, комплексно сопряженная с H(z). Переносим H(z) в левую часть формулы и для функций диракоидного типа заменяем выражение 1/H(z) = H-1(z):

А(z)H-1(z) = H*(z).                                           (6.3.7)

Запишем последнее равенство в развернутом виде:

(a-Nz-N+ ... +a-1z-1+a0+a1z1+ ... +aNzN)(h0-1+h1-1z1+h2-1z2+ ... +hN-1zN) =

= h0*+h1*z-1+h2*z-2+ ... +hN*zN.                                    (6.3.8)

В выражении (6.3.8) сумма коэффициентов при одинаковых степенях z в левой части равенства должна быть равна коэффициенту при соответствующей степени z в правой части равенства, что позволяет составить следующую систему из N уравнений для коэффициентов при степенях z0, z1, z2, ... , zN:

a0 h0-1 +        a-1 h1-1 +       a-2 h2-1 +       a-3 h3-1 + ...   + a-N hN-1         = h0*             (13.3.9)

a1 h0-1 +        a0 h1-1 +        a-1 h2-1 +       a-2 h3-1 + ...   + a-N-1 hN-1     = 0

a2 h0-1 +        a1 h1-1 +        a0 h2-1 +        a5 h3-1 + ...    + a-N-2 hN-1     = 0         (6.3.9)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

aN h0-1 +       aN-1 h1-1 +     aN-2 h2-1 +     aN-3 h3-1 +...  +a0 h-N -1          = 0

 

В случае вещественных фильтров, когда ai = a-i и h0* = h0, уравнение (6.3.9) идентично уравнению (6.3.5).

Уравнение Левинсона. Практический способ расчета оптимальных инверсных фильтров по уравнению (6.3.9) предложен в 1947 году Н.Левинсоном.

Перепишем уравнение (6.3.9) в матричной форме:

                 (6.3.10)

Так как коэффициенты инверсного фильтра достаточно определить с точностью до произвольного масштабного множителя, приведем ho-1 к 1, a функцию автоковариации переведем в функцию коэффициентов корреляции делением обеих частей уравнения на ао. Обозначая Ai = ai/ao, Wi = hi-1/ho-1 и V ho*/(ho-1ao) =  hoho*/ao , получаем:

              (6.3.11)

где для значений W и V введен индекс j номеров предстоящих итераций по циклу вычисления коэффициентов фильтра.

При нулевой итерации (N=0, j=0) имеем только одно уравнение:

.                                           (6.3.12)

Благодаря проведенной нормировке решения уравнения (6.3.12) не требуется:

А0= 1, V0= 1, W00= 1.

Составим уравнение для двучленного фильтра (N=1, j=1):

                                       (6.3.13)

Перепишем уравнение (13.3.12) в прямой форме:

А0 W00 = V0.                                                (6.3.14)

Запишем вспомогательную систему, для чего к уравнению (6.3.14) добавим вторую строку с новой постоянной Ej:

A0 W00 + A1·0 = V0,

A1 W00 + A0 ·0 = E1.

В матричной форме:

 .                                     (6.3.15)

Реверсируем уравнение (13.3.15):

 .                                    (6.3.16)

Вычтем (13.3.16) из (13.3.15) с неопределенным множителем Rj:

 .                 (6.3.17)

Из верхней строки уравнения (7.3.16) можно получить значение Е1:

Е1= A1W00.                                                 (6.3.18)

Уравнение (6.3.13) можно  сделать  равнозначным  уравнению (6.3.17), если правую часть нижней строки уравнения (6.3.17) приравнять к правой части нижней строки уравнения (6.3.13):

E1 - R1V0 = 0,   R1 = E1/V0.                                      (6.3.19)

При этом из правых частей верхних строк уравнений (6.3.13, 6.3.17) будем иметь:

V1 = V0 - R1E1.                                               (6.3.20)

Приравнивая друг другу левые части уравнений (6.3.13, 6.3.17), получаем:

W01 = W00 - R1·0 = W00 = 1.

W11 = 0 - R1W00 = -R1W00.                                     (6.3.21)

Этим заканчивается первая итерация. Аналогично, для второй итерации:

                                (6.3.22)

                                (6.3.23)

                                (6.3.24)

         (6.3.25)

Из верхней строки уравнения (7.3.24):

Е2 = A1W11+A2W01.

Из правых частей нижней и верхней строк уравнений (6.3.22, 6.3.25):

R2 = E2/V1,

V2 = V1 - R2E2.

Новые коэффициенты из левых частей уравнений (6.3.22, 6.3.25):

W02 = W01 - R2  0= 1,

W12 = W11 - R2W11,

W22 = 0 - R2W01.

Анализ расчетов позволяет вывести следующие рекуррентные формулы:

Ej =AiWj-i,j ,   j = 1,2,...,M.                   (6.3.26)

Rj = Ej/Vj-1,

Vj = Vj-1 - RjEj,

Wi,j = Wi,j-1 - RjWj-1,j-1,  i = 0,1,.., j.

 

6.4. Рекурсивная  деконволюция

 

Уравнение фильтра рекурсивной деконволюции. Запишем уравнение (6.1.2) для инверсного фильтра в развернутой форме:

H-1(z) = 1/(h0+h1z+h2z2+ ...).                                    (6.4.1)

Так как для минимально-фазового оператора всегда выполняется условие h0  0, приведем (13.4.1) к виду:

H-1(z) = (1/h0)/(1+h1z/h0+h2z2/h0+...) = G/(1+q1z+q2z2+ ...),         (6.4.2)

где: G = 1/h0, q1 = h1/h0, q2 = h2/h0 и т.д. Но уравнение (6.4.2) есть не что иное, как уравнение передаточной функции рекурсивного фильтра, где цепь обратной связи образована коэффициентами нормированного оператора h(n). Алгоритм вычислений:

yk = G·xkqn·yk-n.

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

 

6.5. Фильтры сжатия сигналов

 

Рассмотренные выше методы относятся к расчету так называемых идеальных инверсных фильтров, т.е. фильтров полной инверсии системных операторов. Однако использование идеальных инверсных фильтров на практике не всегда возможно, т.к. регистрируемые данные обычно осложнены влиянием помех (шумов), а инверсные фильтры обычно имеют коэффициент усиления дисперсии шумов значительно больше 1. В этом случае задача точной деконволюции (восстановления истинной первоначальной формы сигнала), как правило, не ставится, а инверсные фильтры считаются оптимальными с точки зрения максимального приближения к форме полезного сигнала с определенным допустимым коэффициентом усиления дисперсии помех. Такие фильтры называются фильтрами неполной (частичной, ограниченной) деконволюции или фильтрами сжатия сигнала. При проектировании фильтров неполной деконволюции учитываются статистические характеристики помех во входном сигнале и их соотношение со статистическими характеристиками самого входного сигнала.

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

H-1(z) = H*(z)/[|H(z)|2+g2],                                     (6.5.1)

где g2 = k·sh2 - дисперсия шумов в единицах дисперсии оператора hn, sh2 – дисперсия значений оператора hn, (при условии суммы значений оператора, равной 1), k - отношение дисперсии шумов к дисперсии оператора hn. Коэффициент g2 играет роль регуляризирующего фактора при выполнении операции деконволюции информации.

На рисунке 6.5.1 пример формы оператора hn и спектральных функций (6.5.1) при разных значениях параметра g. При g = 0 выражение (6.5.1) обращается в идеальный инверсный фильтр 1/H(z). Во втором крайнем случае, при g2>>|H(z)|2, фильтр (6.5.1) переходит в фильтр, согласованный с сигналом по частотному спектру: H-1(z) = H*(z)/g2, который только максимизирует отношение сигнал/помеха.

На рисунке 6.5.2 приведена форма инверсных операторов, соответствующая их частотным характеристикам на рисунке 6.5.1(В), и результаты свертки инверсных операторов с прямым (для лучшего просмотра графики прямой оператор при свертке сдвинут вправо на 2 значения Dt).

При g=0 коэффициент усиления дисперсии шумов равен 11, при g= с2 > 0.3s2 равен 4.6. Однако снижение усиления дисперсии шумов сопровождается увеличением погрешности приближения, что можно видеть на рисунке 6.5.2 (В), при этом уменьшается амплитуда восстановления импульса Кронекера и появляются осцилляции после импульса. Но при наличии шумов и правильном выборе параметра g общее отношение амплитудных значений сигнал/ шум для оператора по (6.5.1) больше, чем для прямой инверсии по (6.5.1.2), что объясняется более существенным уменьшением коэффициента усиления дисперсии шумов при увеличении параметра g, чем увеличением погрешности приближения.

 

D07-12

 

 

 

 

 

 

 

 

 

 

 

 

 

Рисунок 6.5.1. - Форма инверсных операторов

 

Операторы оптимальных фильтров сжатия сигналов также могут вычисляться с учетом помех. Если сигнал s(k) и помеха статистически независимы, то функция автоковариации сигнала на входе фильтра:

ai = asi + bi,                                                   (6.5.2)

где asi и bi - функции автоковариации сигнала и помех. При помехе типа белого шума функция автоковариации помех представляет собой весовую дельта-функцию в точке 0:

bi = c2di,                                                     (6.5.3)

где с2- дисперсия помех. С учетом этого фактора расчет оптимальных инверсных фильтров может проводиться по вышеприведенным формулам (6.3.5, 6.3.9) с изменением значения коэффициента ао:

ao= ao + c2.                                                   (6.5.4)

На рисунке 6.5.3(А) приведены примеры операторов оптимальных инверсных фильтров, вычисленные по прямому оператору, приведенному на рисунке 6.5.1(А).

 

D07-13

 

 

 

 

 

 

 

 

 

 

Рисунок 6.5.2. - Операторы оптимальных фильтров

        

Значения коэффициента с2 заданы в долях дисперсии прямого оператора. Ввод коэффициента с2 в функцию автоковариации резко уменьшает значения коэффициентов инверсного оператора и, соответственно, уменьшает коэффициент усиления дисперсии помех.

 

D07-14

 

Рисунок 6.5.3. - Операторы оптимальных инверсных фильтров

 

Для приведенного примера, при исходном значении коэффициента усиления дисперсии шумов порядка 12 для с2=0, его значение уменьшается до 1.8 при с2=0.1s2 и становится меньше 1 при с2 > с2 > 0.3s2. Естественно, что общая погрешность приближения деконволюции при этом также существенно изменяется (см. рисунок 6.5.3(В)), но амплитуда значения сигнала на месте импульса Кронекера (там, где он должен быть) изменяется много меньше, чем коэффициент усиления дисперсии шумов, а, следовательно, отношение сигнал/шум при введении коэффициента с2, существенно увеличивается.

 

 

7. Аппроксимация сигналов и функций

 

Формулы сигналов, детально и точно описывающие определенные физические объекты и процессы, могут быть очень сложными и мало пригодными для практического использования, как при математическом анализе физических данных, так и в прикладных задачах, особенно при расчетах ожидаемых результатов измерений и при математическом моделировании физических процессов. Кроме того, практическая регистрация физических данных выполняется, как правило, с определенной погрешностью или с определенным уровнем шумов, которые по своим значениям могут быть выше теоретической погрешности прогнозирования сигналов при расчетах по сложным, хотя и очень точным формулам. Не имеет большого смысла и проектирование систем обработки и анализа сигналов по высокоточным формулам, если повышение точности расчетов не дает ощутимого эффекта в повышении точности обработки данных. Во всех этих условиях возникает задача аппроксимации – представления произвольных сложных функций f(x) простыми и удобными для практического использования функциями j(x) таким образом, чтобы отклонение j(x) от f(x) в области ее задания было наименьшим по определенному критерию приближения. Функции j(x) получили название функций аппроксимации.

Математика очень часто оперирует со специальными математическими функциями решения дифференциальных уравнений и интегралов, которые не имеют аналитических выражений и представляются табличными числовыми значениями  yi  для дискретных значений независимых переменных xi. Аналогичными таблицами {yi, xi} могут представляться и экспериментальные данные. Точки, в которых определены дискретные значения функций или данных, называются узловыми. Однако на практике могут понадобиться значения данных величин совсем в других точках, отличных от узловых, или с другим шагом дискретизации аргументов. Возникающая при этом задача вычисления значений функции в промежутках между узами называется задачей интерполяции, за пределами семейства узловых точек вперед или назад по переменным – задачей экстраполяции или прогнозирования. Решение этих задач также обычно выполняется с использованием аппроксимирующих функций.

Сглаживание статистических данных или аппроксимация данных с учетом их статистических параметров относится к задачам регрессии, и рассматриваются в следующей теме. Как правило, при регрессионном анализе усреднение данных производится методом наименьших квадратов (МНК).

Все вышеперечисленные задачи относятся к задачам приближения сигналов и функций и имеют многовековую историю, в процессе которой сформировались классические математические методы аппроксимации, интерполяции, экстраполяции и регрессии функций. В рамках настоящего курса мы не будем углубляться в строгую математическую теорию этих операций. Все современные математические системы (Mathcad, Matlab, Maple и пр.) имеют в своем составе универсальный аппарат выполнения таких операций, дающий пользователю возможность реализации любых практических задач по обработке данных без отвлечения на теоретические подробности их исполнения. Не имеют также большого смысла особенности их использования в различных системах (крутить педали можно учиться на любой марке велосипеда), поэтому в качестве основной математической системы для примеров будем использовать систему Mathcad.

 

7.1           Приближение сигналов рядами Тейлора.

 

Исторически разложение функций в ряд Тейлора явилось одним из первых методов приближения функций в окрестностях точек х0:

f(x) @ f(x0) + (x-x0) + (x-x0)2 + … + (x-x0)n.

f(x) @ f(x0) +(x-x0)i.

         При разложении функции в окрестностях точки х0=0 ряд Тейлора принято называть рядом Маклорена.

Первый член f(x0) ряда представляет собой отсчет функции в точке х0 и грубое приближение к значениям функции в окрестностях этой точки. Все остальные члены ряда детализируют значения функции в окрестностях точки х0 по информации в соседних точках и тем точнее приближают сумму ряда к значениям функции, чем больше членов суммы участвуют в приближении, с одновременным расширением интервала окрестностей точного приближения. Наглядно это можно видеть на примере двух функций, приведенном на рисунке 7.1.1 (копия расчетов в среде Mathcad с усечением отображения членов длинных рядов f2(x) и f4(x)).

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

 

7.2 Интерполяция и экстраполяция сигалов

 

         Линейная и квадратичная интерполяция являются самыми простыми способами обработки таблиц и выполняются по уравнениям:

f(x)лин = а0 + а1х.         f(x)кв = а0 + а1х + а2х2.

При кусочно-линейной интерполяции вычисления дополнительных точек выполняются по линейной зависимости.

 

18

 

Рисунок 7.2.1.- Примеры разложения функций в ряд Маклорена.

        

    Графически это означает простое соединение узловых точек отрезками прямых. В системе Mathcad  для этого используется функция linterp(X,Y,x), где X и Y – вектора узловых точек.  Функция linterp(X,Y,x) возвращает значение функции при её линейной интерполяции по заданным аргументам х. При небольшом числе узловых точек (менее 10) линейная интерполяция оказывается довольно грубой. Первая производная функции аппроксимации испытывает резкие скачки в узловых точках. Для целей экстраполяции функция linterp не предназначена и за пределами области определения сигнала может вести себя непредсказуемо.

Полиномиальная интерполяция. Линейная и квадратичная аппроксимация являются частным случаем полиномиальной интерполяции с помощью аппроксимирующего полинома:

f(x) = а0 + а1х + а2х2 + … + anxn =ai·xi.                      (7.2.1)

Для выполнения полиномиальной интерполяции достаточно по выражению (7.2.1) составить систему линейных уравнений для n узловых точек и определить n значений коэффициентов ai. При глобальной интерполяции, по всем N точкам задания функции, степень полинома равна N-1. Пример выполнения глобальной интерполяции приведен на рисунке 7.2.1. Равномерной дискретизации данных для интерполяции не требуется. Максимальная степень полинома на практике обычно устанавливается не более 8-10, а большие массивы данных интерполируются последовательными локальными частями.

Для практического использования более удобны формулы аппроксимации, не требующие предварительного определения коэффициентов аппроксимирующих полиномов. К числу таких формул относится интерполяционных многочлен по Лагранжу /40/. При аппроксимации функции у(х) многочленом n-й степени Y(x):

 

18

 

Рисунок 7.2.2.- Интерполяция данных.

        

Y(x) =  +  +…

…+ .                              (7.2.2)

         Пример интерполяции по Лагранжу приведен на рисунке 7.2.2.

Кривые Безье. Для задач аппроксимации широко применяются кривые Безье. Это связано с их удобством как для аналитического описания, так и для наглядного геометрического построения. При использовании кривых Безье в компьютерной графике пользователь может задавать форму кривой интерактивно, двигая опорные точки курсором на экране.

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

(t) = (1-t)(t) + t(t),          (7.2.3)

где нижний индекс - номер точки, верхний индекс - уровень разбиения. Уравнение кривой n-ого порядка задается уравнением: P(t) = (t).

Для примера построим кривую для трех опорных точек P0(t=0), P1 и P2(t=1) на интервале t=[0, 1] (см. рисунок 7.2.3).

Для каждого tÎ[0, 1] определим точку :

= (1-t)P0 + tP1,

= (1-t)P1 + tP2,

 

18

 

Рисунок 7.2.3.- Интерполяция по Лагранжу.

        

= (1-t) + t=(1-t)2P0(t)+2t(1-t)P1(t)+t2P2(t),

и тем самым получим кривую второго порядка.

 

 

 

                

 

 

 

 

 

 

 

Рисунок 7.2.3.- Кривые второго порядка 

 

Аналогичным образом построение кривой Безье с четырьмя опорными точками будет определяться следующими выражениями:

= (1-t)P0 + tP1,

= (1-t)P1 + tP2,

= (1-t)P2 + tP3,

= (1-t) + t=(1-t)2P0(t)+2t(1-t)P1(t)+t2P2(t),

= (1-t) + t=(1-t)2P1(t)+2t(1-t)P2(t)+t2P3(t),

= (1-t) + t=(1-t)3P0(t)+3t(1-t)2P1(t)+ 3t2(1-t)P2(t)+t3P3(t).

            Общая аналитическая запись для кривых Безье по N+1 опорной точке:

PN(t) = Pi ,

= ti (1-t)N-i,

= N!/(i! (N-i)!) – биномиальные коэффициенты.

Кривые Безье всегда проходят через начальную P0 и конечную PN точки. Если рассматривать опорные точки в противоположном порядке, то форма кривой не изменяется. Если опорные точки лежат на одной прямой, то кривая Безье вырождается в отрезок, соединяющий эти точки. Степень многочлена, представляющего кривую в аналитическом виде, на 1 меньше числа опорных точек.

 

7.3 Сплайновая интерполяция.

 

         Сплайн - кусочный многочлен степени K с непрерывной производной степени K - 1 в точках соединения сегментов.

При сплайновой интерполяции обычно используются локальные полиномы не выше третьей степени. Так, например, кубические сплайны проходят через три смежные узловые точки (текущие опорные точки вычислений), при этом в граничных точках совпадают как значения полинома и функции, так и значения их первых и вторых производных. Коэффициенты полиномов, проходящих через три смежные узловые точки, рассчитываются так, чтобы непрерывными были первая и вторая его производные. Линия, которую описывает сплайн-функция, напоминает по форме гибкую линейку, закреплённую в узловых точках. Это создает высокую плавность сплайнового полинома по сравнению с другими методами аппроксимации и наглядно видно на рисунке 7.3.1. Полиномы более высоких порядков чрезмерно громоздки для практики.

Сплайновая аппроксимация может применяться для достаточно быстро изменяющихся функций, не имеющих разрывов функции и ее производных. Основной недостаток сплайнов – отсутствие единого аналитического выражения для описания функции. Попутно заметим также, что результаты экстраполяции функций, как это можно видеть на рисунке 7.3.1, существенно зависят от метода аппроксимации, и, соответственно, к их достоверности нужно подходить достаточно осторожно.

Сплайновая интерполяция обычно применяется в составе математических пакетов по определенной технологии. Так, в системе Mathcad при выполнении сплайновой интерполяции по узловым точкам функции (векторам X и Y) сначала вычисляется вектор (обозначим его индексом S) вторых производных входной функции y(x) по одной из программ:

 

18

 

Рисунок 7.3.1. Сплайновая интерполяция и интерполяция по Лагранжу

        

        

- S:= cspline(X,Y) – возвращает вектор S вторых производных при приближении в опорных точках к кубическому полиному;

- S:= pspline(X,Y) – возвращает вектор S при приближении в опорных точках к параболической кривой;

- S:= lspline(X,Y) – возвращает вектор S при приближении в опорных точках к прямой.

По значениям вектора S функцией interp(S,X,Y,x) вычисляются значения аппроксимирующей функции по аргументам х.

        

7.4 Спектральный метод  интерполяции

 

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

Спектр дискретного сигнала. Допустим, что для обработки задается  произвольный аналоговый сигнал s(t), имеющий фурье-образ S(f). Равномерная дискретизация непрерывного сигнала s(t) с частотой F (шаг Dt = 1/F = q) с математических позиций означает умножение функции s(t) на гребневую (решетчатую) функцию Шq(t) = d(t-kDt):

sq(t) = s(t) Шq(t) = s(t)d(t-kDt) =s(kDt)d(t-kDt).             (7.4.1)

С учетом известного преобразования Фурье гребневой функции Шq(t) Û F×ШF(f) фурье-образ дискретной функции sq(t):

SF(f) = S(f) F×ШF(f).                                         (7.4.2)

ШF(f) =d(f-nF).                                           (7.4.3)

Отсюда, для спектра дискретного сигнала имеем:

SF(f) = F×S(f)d(f-nF) = FS(f-nF).                       (7.4.4)

Спектр дискретного сигнала представляет собой непрерывную периодическую функцию с периодом F, совпадающую с функцией F×S(f) непрерывного сигнала s(t) в пределах центрального периода от -fN до fN, где fN = 1/2Dt = F/2 - частота Найквиста. Частота дискретизации сигнала должна быть минимум в два раза выше максимальной частотной составляющей в спектре сигнала (F = 1/Dt ³ 2fmax). Умножая функцию (7.4.2) на прямоугольную весовую функцию ПF(f), равную 1 в пределах главного частотного, получаем непрерывный спектр в бесконечных по частоте границах, равный спектру F×S(f) в пределах главного частотного диапазона:

F×S(f) = F×[S(f) ШF(f)] ПF(f).          (7.4.5)

Обратное преобразование Фурье этого спектра, с учетом коэффициента F, должно восстанавливать непрерывный сигнал, равный исходному аналоговому сигналу s(t).

На рисунке 7.4.1 приведен пример интерполяции и экстраполяции равномерных по аргументу дискретных данных в сравнении с сплайн-методом и методом по Лагранжу. Исходная аналоговая кривая дискретизирована корректно (fmax < 1/2Dt) и восстановленная по дискретным данным кривая fS(z) полностью ее повторяет. Близкие результаты к исходному сигналу дает также и сплайн-интерполяция, но доверять сплайн-экстраполяции, особенно по концевой части интервала задания данного сигнала, не приходится. Что касается интерполяции по Лагранжу, то можно видеть существенную погрешность интерполяции на концевых частях интервала сигнала и полную ее непригодность для задачи экстраполяции.

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

При нарушении корректности дискретизации данных погрешности интерполяции возрастают практически во всех методах интерполяции, а не только в спектральном методе.

Это можно видеть на рисунке 7.4.2, который полностью повторяет рисунок 7.4.1 с изменением значения только одного, пятого отсчета (уменьшение с 7.84 до 2), что вызывает подъем высоких частот в спектре данных. Следует учитывать, что при интерполяции данных, представляющих собой вырезки из сигнальных функций с определенной постоянной составляющей (сигнал не выходит на нулевые значения на концевых участках интервала задания), а равно и любых данных со скачками функций, при спектральном преобразовании на интерполированном сигнале в окрестностях обрезания данных (и скачков) возникает явление Гиббса. Это можно наглядно видеть сравнением рисунков 7.4.1 и 7.4.3. Данные на рисунке 7.4.1 в рисунке 7.4.3 подняты на 20 единиц постоянной составляющей.

Для исключения этого эффекта можно рекомендовать перед интерполяцией производить определение линейного тренда данных по концевым значениям отсчетов и вычитать его из данных, с последующим восстановлением после интерполяции.

Интерполяционный ряд Котельникова-Шеннона. Произведем обратное преобразование обеих частей равенства (7.4.5). Умножение непрерывного и бесконечного спектра на П-импульс в пределах главного диапазона отобразится в динамической области сверткой двух функций:

F×s(t) = F×sq(t) sinc(pFt).

s(t) = sinc(pFt)s(kDt)d(t-kDt),

Отсюда, с учетом равенства d(t-kDt) sinc(pFt) = sinc[pF(t-kDt)], получаем:

s(t) =s(kDt) sinc[pF(t-kDt)].          (7.4.6)

Эта формула носит название интерполяционного ряда Котельникова-Шеннона и, по существу, является разложением сигнала по системе ортогональных функций sinc(pF(t-kDt)) = sinc(p(t/Dt k)). С другой стороны, эта формула представляет собой свертку дискретной функции данных s(kDt) с непрерывной функцией интегрального синуса. Для больших массивов дискретных данных точность восстановления сигнала обычно ограничивается интервалом задания функции интегрального синуса, по которому устанавливается интервал суммирования.

Из совокупности выше приведенных формул следует, что если для частоты дискретизации  сигнала  справедливо  неравенство F ³ 2fmax, где fmax - наибольшая частота в спектре произвольной непрерывной функции s(t), то функция s(t) может представляться в виде числовой последовательности дискретных значений s(kDt),  k = 0,1,2,..., и однозначно по этой последовательности восстанавливаться, в пределе - без потери точности. В этом и состоит сущность теоремы отсчетов Котельникова-Шеннона.

         На рисунке 7.4.4 приведен пример интерполяции входных данных, повторяющих данные рисунок 7.4.1. Результаты интерполяции, как и следовало ожидать, абсолютно идентичны. Аналогичным образом влияют на результаты усечение и скачки функций (явление Гиббса).

        

7.5Децимация и интерполяция цифровых сигналов.

 

Применительно к цифровым сигналам децимация – уменьшение частоты дискретизации данных с сохранением в новом сигнале всей полезной информации. Интерполяция обратна децимации – увеличение частоты дискретизации также без изменения информации. Цифровая децимация и интерполяция широко используется в современных системах обработки данных с различными ширинами полос и частотами дискретизации, для сжатия и восстановления данных, для уменьшения объемов памяти хранения данных, для увеличения скорости передачи данных, для увеличения производительности систем с обработкой данных на разных скоростях, и т.п.

 

 

Рисунок 7.4.1. - Спектральный метод интерполяции и экстраполяции

 

 

Рисунок 7.4.2.- Спектральный метод интерполяции и экстраполяции с изменением значения пятого отсчета

        

 

Рисунок 7.4.3. – Интерполяция дискретных данных

 

         18

 

Рисунок 7.4.4.- Интерполяция по Котельникову-Шеннону

 

Так, например, реализация узкополосных цифровых КИХ-фильтров требует очень большого числа коэффициентов, и децимация данных может существенно снизить размер фильтра и повысить скорость их обработки.

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

Цифровые методы позволяют выполнить эту операцию в более эффективной форме.

Децимация с целым шагом. Кратная компрессия частоты дискретизации снижает частоту дискретизации входного сигнала x(k) с fd до fd/M путем отбрасывания М-1 отсчетов в каждой последовательной серии из М-отсчетов, т.е. из М-отсчетов оставляется только 1:

y(m) =x(k) d(k-mM),  m = 0, 1, 2, … , K/M.                (7.5.1)

Естественно, что частота Найквиста fN входного сигнала x(k) компрессора для выходного сигнала y(k) также уменьшается в М раз и становится равной fN' = fN/M для выходного сигнала. Соответственно, для полного сохранения после компрессии полезной информации, содержащейся в сигнале x(k), максимальная частота полезной информации  во входном сигнале не должна превышать значения fmax £ fN/2M. В противном случае децимация будет некорректной и в новом главном частотном диапазоне выходного сигнала произойдет искажение спектра полезной информации за счет сложения со спектрами боковых диапазонов. Пример корректной децимации сигнала с М=2 и спектры входного и децимированного сигнала приведены на рисунке 7.5.1.

Цифровые методы позволяют выполнить эту операцию в более эффективной форме.

Входные сигналы кроме полезной информации могут содержать статистические шумы и помехи, распределенные по всему частотному диапазону.

 

14.6.1.GIF

 

Рисунок 7.5.1. - Корректная децимация сигнала с М=2 и спектры входного и децимированного сигнала

 

При децимации шумы и помехи в частотном диапазоне от fN/M до fN входного сигнала зеркально отражаются от fN' нового главного частотного диапазона и их суммирование со спектром нового главного диапазона и полезного сигнала может приводить к увеличению уровня шумов и искажению информации. Для исключения этого эффекта перед конверсией сигнала необходимо выполнять его низкочастотную фильтрацию со срезом на частоте fN/2M.

На рисунке 7.5.2 приведен спектр Z2n сигнала zk с рисунка 7.5.1 (спектр Zn), в который для наглядности эффекта зеркального отражения условно введен только высокочастотный шум  на интервале fs - fN, и соответствующий данному спектру сигнал z2k.

 

14.6.2.GIF

 

Рисунок 7.5.2.- Спектр Z2n сигнала

 

При децимации сигнала z2 с М=2 сначала была выполнена его низкочастотная фильтрация с частотой среза fs, что полностью сняло модельный шум, и при восстановлении сигнала (интерполяции) получен сигнал y1k, полностью соответствующий сигналу zk. При выполнении децимации без предварительной фильтрации восстанавливается сигнал y2k, который по своей форме отличается как от сигнала y1k = zk, так и от входного сигнала z2k. Впрочем, следует заметить, что децимация с предварительной фильтрацией никогда не сохраняет формы входного сигнала при обратной интерполяции, если в подавляемой высокочастотной части входного сигнала присутствуют какие-либо значимые гармоники шумов, помех, или, например, от скачков в сигнале.

         При выполнении децимации с большими значениями коэффициента М операцию рекомендуется выполнять последовательными блоками M = M1 M2 …, что снижает требования к частотным характеристикам фильтров нижних частот.

Интерполяция с целым шагом.  Экспандер частоты дискретизации увеличивает частоту дискретизации входного сигнала xk с fd до Lfd путем введения (L-1) нулевых отсчетов после каждого отсчета входного сигнала. При этом форма спектра выходного сигнала yk остается без изменения, но частотная шкала спектра сжимается в L раз и в границы главного диапазона спектра входного сигнала ±fN заходят боковые диапазоны спектра выходного сигнала (зеркальные частоты). Это наглядно можно видеть на рисунке 7.5.3 сравнением спектров Xn  для входного сигнала xk, и Yn для экспандированного сигнала xk с L=2. Следовательно, фактическая частота Найквиста выходного сигнала становится равной fN/L. Для исключения зеркальных частот и распределения энергии отсчетов xk по L выходным интервалам экспандированный сигнал пропускается через фильтр нижних частот со срезом на частоте fN/L и с коэффициентом L для компенсации распределения энергии отсчетов по интервалам L. Результат операции можно видеть на сигнале yk по сравнению с исходным сигналом zk (см.рисунок 7.5.1), децимацией которого с М=2 был получен сигнал xk.

При выполнении интерполяции с большим значением коэффициента L также рекомендуется выполнять операцию последовательными блоками L=L1 L2 … для снижения требований к частотным характеристикам фильтров низких частот.

 

14.6.3.GIF

 

Рисунок 7.5.3. - Интерполяция с целым шагом

        

Преобразование частоты дискретизации с нецелым шагом  на практике обычно выполняют представлением нецелого множителя максимально близким приближением рациональными числами вида L/M, Это позволяет выполнять преобразование частоты дискретизации последовательными операциями сначала интерполяции с шагом L, сохраняющей все частотные составляющие сигнала, и затем децимации с шагом М, при которой часть высокочастотных составляющих и шумов будет подавлена низкочастотной фильтрацией. Поскольку при этом низкочастотные фильтры экспандирования и децимации следуют друг за другом и работают на одной частоте дискретизации, то вместо двух фильтров можно применять один, имеющий минимальную частоту среза с коэффициентом усиления, равным L.

При программной обработке больших пакетов данных децимация и интерполяция может выполняться в спектральной области с использованием БПФ.

 

7.6 Методика аппроксимации эмпирических данных

 

Эмпирические данные, как правило, задаются числовыми рядами значений двух величин: независимой (yk) и зависимой (xk) , каждая из которых в общем случае, кроме определенной регулярной (детерминированной) составляющей, может содержать и случайные составляющие самой различной природы, обусловленные как статистической природой изучаемых процессов, так и внешними факторами процессов измерений и преобразования данных (шумы, помехи, дестабилизирующие факторы и ошибки измерений). Независимая переменная xk обычно полагается детерминированной, а, следовательно, ее случайная составляющая "переносится" на зависимую переменную yk. Полагается также, что значения случайной составляющей зависимой переменной (как собственные, так и "суммарные") распределены по некоторому вероятностному закону (например – нормальному).

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

- выбора меры близости зависимой переменной к искомой функции и метода построения приближения (параметров математической модели);

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

- метода оптимизации порядка модельной функции или числа членов ряда аппроксимирующего выражения.

Отсюда следует, что оптимальная аппроксимация может быть обеспечена только достаточно гибкими интерактивными алгоритмами на основе многоэтапных итерационных процессов с возможностью коррекции на каждом этапе.

Мера приближения. Наиболее распространен критерий наилучшего приближения в виде минимума степенной разности между значениями переменной yk и аппроксимирующей функцией j(xk):

[yk - j(xk)]S ® min,                                      (7.6.1)

где S > 0 - положительное число.

Квадратичная мера реализуется при S = 2 в методе наименьших квадратов (МНК) и обеспечивает максимальное правдоподобие функции приближения при нормальном распределении случайной составляющей зависимой переменной yk. Несмещенной оценкой меры приближения в МНК является дисперсия остатков:

D = {[yk - j(xk)]2}/(k-m),                                (7.6.2)

где

-  m – количество параметров в функции приближения;

-  (k-m) – число степеней свободы.

Однако  эмпирические данные могут содержать выбросы и грубые ошибки, которые вызывают смещения вычисляемых параметров. Их влияние обычно исключается цензурированием данных: вычислением гистограммы разностей yk-j(xk) после определения первого приближения функции аппроксимации и исключением "хвостовых" элементов гистограммы (до 2.5% от количества данных, или резко выделяющихся элементов данных на основании оценок вероятностей с использованием r- или t- распределений).

         Мера наименьших модулей (метод Лагранжа) реализуется при S = 1 и применяется при распределениях случайных составляющих зависимой переменной по законам, близким к закону Лапласа (двустороннее экспоненциальное распределение). Такая мера соответствует площади между графиками эмпирических данных и функции аппроксимации и, по сравнению с квадратической, является более устойчивой, в том числе при наличии случайных составляющих с большими амплитудами (длинные "хвосты" разностных гистограмм). Оценки по модулю получили название "робастных" (robust – устойчивый).

         Свойства квадратичной меры и меры наименьших модулей в определенной степени сочетаются при S = 3/2.

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

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

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

         Наибольшее распространение в практике аппроксимации при отсутствии теоретических аспектов изучаемых явлений получили функциональные ряды, для которых определяющее значение имеет порядок аппроксимирующей функции (модели).

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

         Очевидно, что для функциональных рядов порядок модели (степень ряда для степенных рядов) определяет значение меры приближения. При повышении порядка модели (в пределе до бесконечности) минимум функции (7.6.1) стремится к нулю. Однако это означает, что при повышении порядка модели в функцию аппроксимации входит не только регулярная составляющая данных, но все большая и большая доля случайных составляющих, в пределе до полного соответствия функции jk исходным данным yk. Но повышение степени приближения к исходным данным при наличии в них случайных составляющих с какого-то определенного момента (порядка модели) не только не будет приближать функцию аппроксимации к регулярным составляющим данных, а наоборот – увеличивать расхождение. С этой точки зрения термин "меры приближения" (7.6.1) было бы целесообразнее заменить термином "мера аппроксимации" данных, а под мерой приближения понимать значение меры аппроксимации, при которой обеспечивается максимальная степень приближения функции аппроксимации к регулярной составляющей данных (минимум дисперсии разности функций аппроксимации и регулярной составляющей).

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

[yk-j(xk)]2 =[sk+sk-j(xk)]2 =[sk-j(xk)]2 +2[sk-j(xk)]sk +sk2.

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

[sk-j(xk)]2 ® min,                                      (7.6.3)

[yk-j(xk)]2 ®sk2.                                   (7.6.4)

         В пределе, при идеальной аппроксимации, выражение (7.6.3) стремится к нулю, а выражение (7.6.4) эквивалентно соотношению дисперсий:

{[yk - j(xk)]2}/(k-m) »sk2/k.                    (7.6.5)

         Отсюда следует, что при прочих равных условиях наилучшим является приближение, у которого мера приближения близка к дисперсии шума. Для "белых" шумов оценку их дисперсии в экспериментальных данных можно выполнять в спектральной области, если частота Найквиста данных минимум в 2 раза выше предельных частот регулярной составляющей.

         При отсутствии информации о дисперсии шумов оптимальный порядок модели может определяться методом последовательных уточнений с последовательным нарастанием порядка модели и сравнением по критерию Фишера значимости различия дисперсии остатков каждого нового порядка с предыдущим порядком. При увеличении порядка модели (начиная с 1-го) значимость различия дисперсий сначала является довольно высокой, постепенно уменьшается, и в области оптимальных порядков становится малозначимой. Это объясняется тем, что в этой области при небольших уменьшениях значения числителя выражения (7.6.2) одновременно, за счет увеличения порядка, сокращается число степеней свободы. После прохождения оптимальной зоны значения дисперсий остатков снова начинают увеличиваться с увеличением значимости различий.  

Оценка качества приближения. Для оценки качества математической модели эмпирической зависимости используется коэффициент детерминации (Adjusted R2):

Adjusted R2 = Dj/Dy = 1 – Do/Dy,

где:

- Dj - дисперсия функции приближения;

- Dy – дисперсия данных;

-Do – дисперсия остатков.

Чем выше качество аппроксимации, тем ближе к 1 значение коэффициента детерминации.

8.    РЕГРЕССИЯ

 

Аппроксимация данных с учетом их статистических параметров относится к задачам регрессии. Они обычно возникают при обработке экспериментальных данных, полученных в результате измерений процессов или физических явлений, статистических по своей природе (как, например, измерения в радиометрии и ядерной геофизике), или на высоком уровне помех (шумов). Задачей регрессионного анализа является подбор математических формул, наилучшим образом описывающих экспериментальные данные.

Термин "регрессия" появился при исследовании соотношения роста родителей и их детей, в которых было установлено, что рост "регрессирует" к среднему, т.е. высокие родители имеют более низких детей, а низкие родители – более высоких. 

В качестве основной математической системы для примеров будем использовать систему Mathcad.

 

8.1            Постановка задачи регрессии

 

Математическая постановка задачи регрессии заключается в следующем. Зависимость величины (числового значения) определенного свойства случайного процесса или физического явления Y от другого переменного свойства или параметра Х, которое в общем случае также может относиться к случайной величине, зарегистрирована на множестве точек xk множеством значений yk, при этом в каждой точке зарегистрированные значения yk и xk отображают действительные значения Y(xk) со случайной погрешностью sk, распределенной, как правило, по нормальному закону. По совокупности значений yk требуется подобрать такую функцию f(xk, a0, a1, … , an), которой зависимость Y(x) отображалась бы с минимальной погрешностью. Отсюда следует условие приближения:

yk = f(xk, a0, a1, … , an) + sk.

Функцию f(xk, a0, a1, … , an) называют регрессией величины y на величину х. Регрессионный анализ предусматривает задание вида функции f(xk, a0, a1, … , an) и определение численных значений ее параметров a0, a1, … , an, обеспечивающих наименьшую погрешность приближения к множеству значений yk. Как правило, при регрессионном анализе погрешность приближения вычисляется методом наименьших квадратов (МНК). Для этого выполняется минимизация функции квадратов остаточных ошибок:

s,(a0, a1, … , an) =[f(xk, a0, a1, … , an) - yk]2.

         Для определения параметров a0, a1, … , an функция остаточных ошибок дифференцируется по всем параметрам, полученные уравнения частных производных приравниваются нулю и решаются в совокупности относительно всех значений параметров. Виды регрессии обычно называются по типу аппроксимирующих функций: полиномиальная, экспоненциальная, логарифмическая и т.п

 

 

 

8.2           Линейная  регрессия

 

Общий принцип. Простейший способ аппроксимации по МНК произвольных данных sk - с помощью полинома первой степени, т.е. функции вида y(t) = a+bt, которую обычно называют линией регрессии. С учетом дискретности данных по точкам tk, для функции остаточных ошибок имеем:

s (a, b) =[(a+b tk) - sk]2.

Для вычисления оценок коэффициентов дифференцируем функцию остаточных ошибок по аргументам a и b, приравниваем полученные уравнения нулю и формируем два нормальных уравнения системы:

2((a+b tk)-sk) º a1 + btksk = 0,

2((a+b tk)-sk) tk º atk + btk2sk tk = 0,

Решение данной системы уравнений в явной форме для К-отсчетов:

b = [Ktk sktksk] / [Ktk2 – (tk)2] = (- ) / (- ).

a = [skbtk] /K = - b

Полученные значения коэффициентов используем в уравнении регрессии y(t) = a+bt.  Прямая (s) = b (t - ) называется линией регрессии s по t. Для получения линии регрессии t по s,  (t - )  = b (s), аргумент b в этой формуле заменяется на значение b = (- ) / (- ).

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

Реализация в Mathcad. Линейная регрессия в системе Mathcad выполняется по векторам аргумента Х и отсчетов Y функциями:

- intercept(X,Y) – вычисляет параметр а, смещение линии регрессии по вертикали;

- slope(X,Y) – вычисляет  параметр b, угловой коэффициент линии регрессии.

         Расположение отсчетов по аргументу Х произвольное. Функцией corr(X,Y) дополнительно можно вычислить коэффициент корреляции Пирсона. Чем он ближе к 1, тем точнее обрабатываемые данные соответствуют линейной зависимости.

         Пример выполнения линейной регрессии приведен на рис. 8.2.1.

 

8.3            Полиномиальная  регрессия

 

         Одномерная полиномиальная регрессия с произвольной степенью n полинома и с произвольными координатами отсчетов в Mathcad выполняется функциями:

- regress(X,Y,n) – вычисляет вектор S для функции interp(…), в составе которого находятся коэффициенты ki полинома n-й степени;

- interp(S,X,Y,x) – возвращает значения функции аппроксимации по координатам х.

Функция interp(…) реализует вычисления по формуле:

f(x) = k0 + k1 x1 + k2 x2 + … + kn xnki xi.

18

 

Рисунок 8.2.1.- Линейная регрессия

 

         Значения коэффициентов ki могут быть извлечены из вектора S функцией  submatrix(S, 3, length(S), 0, 0).

На рисунке 8.3.1 приведен пример полиномиальной регрессии с использованием полиномов 2, 3 и 8-й степени. Степень полинома обычно устанавливают не более 4-6 с последовательным повышением степени, контролируя среднеквадратическое отклонение функции аппроксимации от фактических данных. Нетрудно заметить, что по мере повышения степени полинома функция аппроксимации приближается к фактическим данным, а при степени полинома, равной количеству отсчетов минус 1, вообще превращается в функцию интерполяции данных, что не соответствует задачам регрессии. 

Зональная регрессия. Функция regress по всей совокупности точек создает один аппроксимирующий полином. При больших координатных интервалах с большим количеством отсчетов и достаточно сложной динамике изменения данных рекомендуется применять последовательную локальную регрессию отрезками полиномов малых степеней. В Mathcad это выполняется

отрезками полиномов второй степени функцией loess(X, Y, span), которая формирует специальный вектор S для функции interp(S,X,Y,x). Аргумент span > 0 в этой функции (порядка 0.1-2) определяет размер локальной области и подбирается с учетом характера данных и необходимой степени их сглаживания (чем больше span, тем больше степень сглаживания данных).

 

8.4           Нелинейная  регрессия

 

Линейное суммирование произвольных функций. В Mathcad имеется возможность выполнения регрессии с приближением к функции общего вида в виде весовой суммы функций fn(x):

f(x, Kn) = K1 f1(x) + K2 f2(x) + … + KN fN(x),

при этом сами функции fn(x) могут быть любого, в том числе нелинейного типа. С одной стороны, это резко повышает возможности аналитического отображения функций регрессии. Но, с другой стороны, это требует от пользователя определенных навыков аппроксимации экспериментальных данных комбинациями достаточно простых функций.

 

2

 

Рисунок 8.3.1.- Одномерная полиномиальная регрессия.

 

Реализуется обобщенная регрессия по векторам X, Y и f функцией- linfit(X,Y,f), которая вычисляет значения коэффициентов Kn. Вектор f должен содержать символьную запись функций fn(x). Координаты xk в векторе Х могут быть любыми, но расположенными в порядке возрастания значений х (с соответствующими отсчетами значений yk в векторе Y). Пример выполнения регрессии приведен на рисунке 8.4.1. Числовые параметры функций f1-f3 подбирались по минимуму среднеквадратического отклонения.

Регрессия общего типа. Второй вид нелинейной регрессии реализуется путем подбора параметров ki к заданной функции аппроксимации с использованием функции genfit(X,Y,S,F), которая возвращает коэффициенты ki, обеспечивающие минимальную среднюю квадратическую погрешность приближения функции регрессии к входным данным (векторы Х и Y координат и отсчетов).

Символьное выражение функции регрессии и символьные выражения ее производных по параметрам ki записываются в вектор F. Вектор S содержит начальные значения коэффициентов ki для решения системы нелинейных уравнений итерационным методом. Пример использования метода приведен на рисунке 8.4.2.

Типовые функции регрессии Mathcad. Для простых типовых формул аппроксимации предусмотрен ряд функций регрессии, в которых параметры функций подбираются программой Mathcad самостоятельно. Как можно видеть из сопоставления методов по средним квадратическим приближениям к базовой кривой и к исходным данным, известность функции математического ожидания для статистических данных с ее использованием в качестве базовой для функции регрессии дает возможность с более высокой точностью определять параметры регрессии в целом по всей совокупности данных, хотя при этом кривая регрессии не отражает локальных особенностей фактических отсчетов данной реализации. Это имеет место и для всех других методов с заданием функций регрессии.

 

 

2

 

 

 

 

 

 

 

 

 

 

Рисунок  8.4.1. - Обобщенная регрессия.

 

2

 

Рис. 8.4.2.- Итерационный метод решения системы нелинейных уравнений

 

 

 

 

 

 

 

 

 

Список литературы

 

1.                 Баскаков С. И. Радиотехнические цепи и сигналы. – М.: Высшая школа, 2003. – 462 с.

2.                 Гоноровский И. С. Радиотехнические цепи и сигналы. – М.: Радио и связь, 1986. – 512 с.

3.                 Карамов З. С., Колесниченко Г. И. Цифровая обработка сигналов: Учебное пособие/ МИС. – М.: 1990.- 41 с.

4.                 Куприянов М.С. и др. Техническое обеспечение цифровой обработки сигналов: Справочник. – СПб. Форт, 2000.- 752 с .

5.                 Куприянов М.С., Матюшкин Б. Д.  Цифровая обработка сигналов: процессоры, алгоритмы, средства проектирования. – СПб., 1999. – 592 с.

6.                 Солонина А. И., Улахович Д. А. и другие. Основы цифровой обработки сигналов. – С-П. «БХВ-Петербург», 2003. – 608 с.

7.                 Введение в цифровую фильтрацию/ Под ред. Р. Бочнера и А. Константинидиса. – М.: Изд-во” Мир”, 1976. – 215 с.

8.                 Теоретические основы радиотехники/ Под ред. В. Н. Ушакова. – М.: Высшая школа, 2002. – 306 с.

9.                 Радиотехнические цепи и сигналы/ Под ред. К. А. Самойло. – М.: Радио и связь, 1982. – 528 с.

10.            Гольденберг Л.М.  и др. Цифровые устройства и микропроцессорные системы. Задачи и упражнения: Учебное пособие. – М.: Радио и связь, 1992. – 256 с.

11.            Петрищенко С.Н., Мусапирова Г.Д. Основы цифровой обработки сигналов. Методические указания к выполнению лабораторных работ для студентов всех форм обучения специальностей факультета радиотехники и связи. – Алматы, 2004. – 33 с.

12.            Казиева Г. С. Основы цифровой обработки сигналов в телекоммуникационных системах: Конспект лекций. – Алматы: АИЭС, 2006.

13.            Казиева Г.С,. Сарженко Л. И. Методические указания к расчетно- графическим  работам для студентов очной формы обучения. Основы цифровой обработки сигналов в телекоммуникационных сетях. - Алматы: АИЭС, 2007.

14.            Казиева Г.С. Телекоммуникациялық жүйелерде сигналды цифрлық     өңдеу негіздері. Дәрістер жинағы.-А., 2007.

15.            Айфичер Э., Джервис Б. Цифровая обработка сигналов. Практический подход.  М., "Вильямс", 2004-  992 с.

16.             Машеров Е.  Цифровая обработка сигналов – некоторые основные понятия. http://www.nsi.ru/~EMasherow/DSP.htm.

17.            Давыдов А.В. Теория сигналов и систем. http://prodav.narod.ru/signals/index.html.