Дискретные сигналы и их типы. Дискретные системы и их типы. Линейные инвариантные к сдвигу системы

МОСКОВСКИЙ АВИАЦИОННЫЙ ИНСТИТУТ (НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ) Факультет радиоэлектроники летательных аппаратов Кафедра № 402 Материал к лекционным занятиям по дисциплине «Наименование дисциплины» Москва, 2017 г. ЛЕКЦИЯ 1. ДИСКРЕТНЫЕ СИГНАЛЫ И ИХ ТИПЫ. ДИСКРЕТНЫЕ СИСТЕМЫ И ИХ ТИПЫ. ЛИНЕЙНЫЕ ИНВАРИАНТНЫЕ К СДВИГУ СИСТЕМЫ. 1. Модель дискретной линейной системы Теория линейной (аналоговой) цепи базируется на электрических свойствах индуктивностей, емкостей и сопротивлении, которые приводят через законы Кирхгофа к описанию цепей с помощью линейных дифференциальных уравнений с постоянными коэффициентами. В отличие от этого, теория дискретных или цифровых линейных систем базируется на линейных разностных уравнениях с. постоянными коэффициентами, которые могут быть решены с помощью операций над числами на специализированных или универсальных вычислительных машинах. При реализации алгоритма разностного уравнения существенным моментом является то, что входной сигнал представляется в виде совокупности дискретных отсчетов по времени. Для многих насущных задач может быть применена модель, показанная на рис. 1. В этих задачах сигналы, подвергаемые обработке, являются аналоговыми (например, выходные сигналы микрофона или сейсмографа), но сама обработка должна быть осуществлена в цифровом виде. Однако не все задачи относятся к этой категории; например, цифровой синтезатор речи не нуждается в аналоговом входном сигнале. На Рис. 3.1 для того, чтобы получить последовательность х(nТ), аналоговый сигнал дискретизуется через временные интервалы Т. Рис. 1 - Модель цифровой машины для обработки аналоговых сигналов Последствия такого процесса дискретизации можно истолковать различными способами: почти во всей литературе, ссылки на которую приведены в гл. 1, вскрывается их сущность, се условно можно назвать наложением от английских понятий folding или aliasing, встречающихся в литературе. В терминах, относящихся к частотной области, это означает, что спектр аналогового сигнала повторяется периодически в виде верхнего и нижнего боковых спектров около частот 1Т, 2/Т, 3/Т и т. д. Если частота дискретизации взята ниже частоты Найквиста (т.е. если 1/Т меньше удвоенной ширины полосы сигнала), то соседние спектры перекрываются, вызывая частотные искажения, которые делают невозможным восстановление аналогового сигнала с помощью линейной фильтрации. Квантизатор на рис. 1 является обязательной частью любого современного цифрового процессора. Необходимо подчеркнуть, что эффекты квантования наблюдаются также и внутри цифровой машины, показанной на рис. 1. Они встречаются, например, при записи постоянного коэффициента линейного разностного уравнения в виде числа в регистре конечной длины. Такой эффект может считаться статическим в том случае, когда свойства разностного уравнения установлены раз и навсегда. Существуют также динамические эффекты, когда сигналы умножаются на коэффициенты или на другие сигналы и произведение округляется или ограничивается до заданной конечной длины регистра. В этой главе мы сделаем важное предположение, что подобными ошибками можно пренебречь. Если требуется получить выходной сигнал в непрерывной форме то отсчеты выходного сигнала цифровой машины проходят через декодер или затягивающее устройство, которое создает непрерывный сигнал из последовательности импульсов. Вообще говоря, декодер - это преобразователь цифра - аналог, после которого включен линейный аналоговый фильтр, устраняющий лишние частоты, появившиеся в процессе дискретизации. Так, например, если Г=10~4 сек, декодер может быть низкочастотным фильтром с частотой среза 5 кгц. В самом простом случае таким фильтром является схема, запоминающая отсчетные значения в паузах между ними. Исследование непрерывных линейных динамических систем в значительной степени облегчилось благодаря введению операционных методов преобразований Лапласа и Фурье, а также благодаря использованию теории цепей. Точно так же исследованию линейных дискретных систем помогает введение z-преобразования и использование теории цепей. Читателям, знакомым с непрерывными линейными системами, многие теоретические выкладки для дискретных линейных систем покажутся хорошо известными, да они и в самом деле должны казаться такими. Но необходимо помнить, что здесь мы имеем дело с дискретной или цифровой областью, для которой должны быть развиты новые подходы, и попытки распространить на эту область результаты, полученные для аналоговой области, могут привести к серьезным ошибкам. 2. Линейные системы, инвариантные к сдвигу Система определяется математически как однозначное преобразование или оператор, отображающий входную последовательность х(n) (вход) в выходную у(n) (выход), что математически записывается в виде у(n) =Т[х(n)], а графически часто изображается так, как показано на рис. 2. Рис. 2 - Представление преобразования, отображающего входную последовательность х(n) в выходную последовательность у(n) Классы дискретных систем определяются путем наложения ограничений на преобразование Т. В дальнейшем будет широко рассматриваться класс линейных инвариантных относительно сдвига систем, потому что они сравнительно просты в математическом отношении, а также потому, что они дают удобный вид обработки сигналов. Позже мы обсудим более общий класс систем, включающий как частный случай линейные системы. Класс линейных систем определяется принципом суперпозиции. Если y1(n) и y2(n) -отклики на x1(n) и x2(n) соответственно, то система линейна тогда и только тогда, когда (1) для произвольных постоянных a и b. Мы видели, что произвольная последовательность х(n) может быть представлена в виде задержанной и взвешенной суммы единичных импульсов. Это представление вместе с (1) предполагает, что линейная система может быть полностью охарактеризована откликом на единичный импульс - импульсной характеристикой. А именно, пусть hk(n) - отклик системы на единичный импульс в момент n=k. (2) Таким образом, согласно (2) реакцию системы можно выразить через отклики на. Если накладывается только одно условие - линейность, то hk(n) будет зависеть как от n, так и от k, n в этом случае польза от выражения (2) для вычислений невелика. Более полезный результат получится, если мы наложим дополнительное ограничение, состоящее в инвариантности к сдвигу. Класс инвариантных к сдвигу систем характеризуется следующим свойством: если у(n) -отклик на х(n), то у(n-k) будет откликом па х(n-k), где k - положительное или отрицательное целое число. Когда индекс n связывается со временем, свойству инвариантности к сдвигу соответствует свойство инвариантности во времени. Из свойства инвариантности к сдвигу следует, что если h(n) -отклик на , то откликом на будет просто . Поэтому (3.2) принимает вид (3) Значит, любая инвариантная к сдвигу система полностью характеризуется импульсной характеристикой h(n). Выражение (3) обычно называется сверткой. Если у(n) -последовательность, значения которой связаны со значениями, двух последовательностей h(n) и х(n) выражением (3), то мы говорим, что у(n) есть свертка х(n) с h(n), и обозначаем у(n) = х(n)*h(n). Заменяя переменную в (3.3), получим другое выражение (4) Поэтому порядок, в котором две последовательности входят в свертку, не важен. Другими словами, линейная инвариантная к сдвигу система со входом х(n) и импульсной характеристикой h(n) будет иметь тот же выход, что и линейная инвариантная к сдвигу система со входом h(n) и импульсной характеристикой х(n). Две линейные инвариантные к сдвигу системы, включенные каскадно, образуют линейную инвариантную к сдвигу систему с импульсной характеристикой, равной свертке импульсных характеристик исходных систем. Так как порядок в свертке не важен, то импульсная характеристика составной системы не зависит от порядка, в котором включены исходные системы. Это свойство иллюстрируется рис. 3, где изображены три системы, имеющие одинаковые импульсные характеристики. Из (3) и (4) следует, что две инвариантные к сдвигу системы, включенные параллельно, эквивалентны одной системе с импульсной характеристикой, равной сумме импульсных характеристик исходных систем. Это свойство иллюстрируется рис. 4, 5. Рис. 3 - Три линейные инвариантные к сдвигу системы с одинаковыми импульсными характеристиками Рис. 4 - Параллельное включение линейных инвариантных к сдвигу систем и эквивалентная система Хотя выражение свертки в виде суммы аналогично интегралу Рис. 5 - Последовательности, входящие в свертку [h(n-k)], показаны для нескольких значений n Свертки в теории линейных аналоговых систем, следует подчеркнуть, что свертку в виде суммы нельзя понимать как приближение к интегралу свертки. В противоположность интегралу свертки, играющему в основном теоретическую роль в применении к аналоговым линейным системам, свертка в виде суммы, как мы увидим в дальнейшем, вдобавок к своей теоретической значимости может служить для реализации дискретной системы. Поэтому важно глубже понять свойства свертки и получить опыт в применении свертки для вычислений. 3. Устойчивость и физическая реализуемость Мы видели, что требования линейности и инвариантности во времени определяют класс систем, которые представляются сверткой. Добавочные ограничения устойчивости и физической реализуемости определяют практически важный, но более узкий класс линейных инвариантных во времени систем. Устойчивой системой назовем систему, в которой каждый ограниченный входной сигнал создает ограниченный выходной сигнал. Линейная инвариантная к сдвигу система устойчива тогда и только тогда, когда (5) Это можно показать следующим образом. Если (5) справедливо и х ограничено, т. е. |х(n)|<М для всех n, то из (4) следует Поэтому у ограничено. Доказать обратное можно, показав, что если , то существует ограниченный входной сигнал, который создает неограниченный выходной сигнал. Таким входом является последовательность со значениями где h*(n) -комплексно-сопряженная к h(n) величина. Ясно, что х(n) ограничена. Значение на выходе при n=0 Поэтому при выходная последовательность не ограничена. Физически реализуемая система - это система, у которой изменения на выходе не опережают изменения на входе, т. е. в физически реализуемой системе, если x1(n)=x2(n), n1-ряд расходится. Следовательно, система устойчива только при |а|<1. 4. Разностные уравнения Системы, у которых входная и выходная последовательности х(n) и у(n) связаны линейным разностным уравнением с постоянными коэффициентами, образуют подмножество класса линейных систем с постоянными параметрами. Описание ЛПП-систем разностными уравнениями очень важно, так как оно часто позволяет найти эффективные способы построения таких систем. Более того, по разностному уравнению можно определить многие характеристики рассматриваемой системы, включая собственные частоты и их кратность, порядок системы, частоты, соответствующие нулевому коэффициенту передачи, и т.д. В самом общем случае линейное разностное уравнение М-то порядка с постоянными коэффициентами, относящееся к физически реализуемой системе, имеет вид (6) где коэффициенты {bi} и {ai} описывают конкретную систему, причем ам ≠ 0. Каким именно образом порядок системы М характеризует математические свойства разностного уравнения, будет показано ниже. Уравнение (6) записано в виде, удобном для решения методом прямой подстановки. Имея набор начальных условий [например, х (i), у(i) для i = -1, -2, ... , -М] и входную последовательность х(n), по формуле (3.6) можно непосредственно вычислить выходную последовательность у(n) для . Например, разностное уравнение (7) с начальным условием у (-1) = 0 и х(n) = n2 + n можно решить подстановкой, что дает Хотя решение разностных уравнений прямой подстановкой и целесообразно в некоторых случаях, значительно полезнее получить решение уравнения в явном виде. Методы нахождения таких решений подробно освещены в литературе по разностным Уравнениям, и здесь будет дан лишь краткий обзор. Основная идея сводится к получению двух решений разностного уравнения: однородного и частного. Однородное решение получается путем подстановки нулей вместо всех членов, содержащих элементы входной последовательности х(n), и определения отклика при нулевой входной последовательности. Именно этот класс решений описывает основные свойства заданной системы. Частное решение получают, подбирая вид последовательности у(n) на выходе при заданной входной последовательности х(n). Для определения произвольных постоянных однородного решения используются начальные условия. В качестве примера решим этим методом уравнение (7). Однородное уравнение имеет вид (8) Известно, что характеристическими решениями однородных уравнений, соответствующих линейным разностным уравнениям с постоянными коэффициентами, являются решения вида . Поэтому, подставив в уравнение (8) вместо у(n), получим (9) Частное решение, соответствующее входной последовательности х(n) = n2 +n, попробуем найти в виде (10) Из уравнения (7) получаем (11) Поскольку коэффициенты при равных степенях n должны совпадать, В, С и D должны быть равны: (12) Таким образом, общее решение имеет вид (13) Коэффициент А определяется из начального условия у (-1) = 0, откуда А = -9/32 и (14) Выборочная проверка решения (14) при показывает полное его совпадение с приведенным выше прямым решением (рис. 6). Очевидное преимущество решения (14) состоит в том, что оно позволяет весьма просто определить у(n) для любого конкретного n = n0. Рис. 6 - Схема реализации простого разностного уравнения первого порядка Важное значение разностных уравнений состоит в том, что они непосредственно определяют способ построения цифровой системы. Так, разностное уравнение первого порядка самого общего вида (15) можно реализовать с помощью схемы, изображенной на рис. 1.6. Блок «задержка» осуществляет задержку на один отсчет. Рассмотренная форма построения системы, в которой для входной и выходной последовательностей используются раздельные элементы задержки, называется прямой формой 1. Ниже мы обсудим различные методы построения этой и других цифровых систем. Разностное уравнение второго порядка самого общего вида (16) Рис. 7 - Схема реализации разностного уравнения второго порядка может быть реализовано с помощью схемы, приведенной на рис. 7. В этой схеме для входной и выходной последовательностей также используются раздельные элементы задержки. Из последующего изложения материалов этой главы станет ясно, что системы первого и второго порядка могут быть использованы при реализации систем более высокого порядка, так как последние могут быть представлены в виде последовательно или параллельно соединенных систем первого и второго порядка. 5. Частотная характеристика В предыдущих разделах рассматривался отклик ЛПП-систем на произвольные входные последовательности. В данном разделе для описания ЛПП-систем в частотной области будет использован специальный класс входных последовательностей, имеющих вид . Как будет показано, этот класс последовательностей является набором собственных функций ЛПП-систем дискретного времени, т. е. для них выходная последовательность совпадает с входной, умноженной на некоторый комплексный коэффициент, зависящий только от . Рассмотрим класс входных последовательностей вида (17) Если такая последовательность поступает на вход ЛПП-системы с импульсной характеристикой h(n), то на выходе появится последовательность (18) (19) (20) Таким образом, для выбранного класса входных последовательностей отклик совпадает с входной последовательностью с точностью до комплексного множителя, который выражается через импульсную характеристику системы следующим образом: (21) Поскольку последовательность вида функционально эквивалентна дискретизованной синусоиде с частотой , то множитель называют частотной характеристикой системы, так как он представляет коэффициент передачи ЛПП-системы для каждого значения . Рис. 8 - Импульсная и частотная характеристики системы первого порядка Вычислим в качестве примера частотную характеристику ЛПП-системы с импульсной характеристикой. Частотная характеристика имеет вид (22) Так как |а| < 1, то сумма геометрической прогрессии (22) будет равна (23) На рис. 8 графически представлены h(n), а также модуль и фаза как функции частоты в диапазоне . Отметим некоторые свойства частотной характеристики. Нетрудно заметить, что частотная характеристика является периодической функцией , причем ее период равен . Эта периодичность связана со спецификой дискретизованного колебания: входная последовательность с частотой не отличается от входной последовательности с частотой , т. е. (24) Поскольку - периодическая функция, то для полного описания достаточно задать ее на любом интервале длиной . Обычно для этой цели используют интервал. Другим важным свойством частотной характеристики является то, что для действительных h(n) (как обычно и бывает на практике) модуль симметричен, а фаза антисимметрична на интервале . Аналогично действительная часть симметрична, а мнимая - антисимметрична на том же интервале. Поэтому при действительных импульсных характеристиках интервал частот, на котором задают частотную характеристику, обычно сокращают до . 6. Замечания о единицах измерения частоты Часто возникает необходимость выразить спектральный состав последовательности h(nT) в единицах частоты, связанных с интервалом дискретизации Т. В этом случае равенство (21) преобразуются к виду (25) (26) Функция периодична по частоте с периодом, равным . Частота в (25) и (26) выражается в радианах в секунду. Характеристику (25) можно выразить и через частоту f, измеряемую в герцах, если заменить на . КОНТРОЛЬНЫЕ ВОПРОСЫ 1. Опишите модель дискретной линейной системы 2.Кратко охарактеризуйте линейные системы, инвариантные к сдвигу 3.Графически изобразите модель цифровой машины для обработки аналоговых сигналов 4. Какой показатель частоты Найквиста? 5.Дайте определение принципам суперпозиции ЛЕКЦИЯ 2. РАЗЛОЖЕНИЕ В РЯД ФУРЬЕ И ЕГО СВОЙСТВА. ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ ФУРЬЕ (ДПФ). СВОЙСТВА ДИСКРЕТНОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ Общие сведения Как известно, звуковой сигнал в компьютере может представляться в виде некоторого набора отсчётов его амплитуд, производимых через определённые промежутки времени (период дискретизации) и представляемых некоторым количеством двоичных разрядов (разрядность выборки). Такое представление удобно для хранения звукового сигнала и его преобразования обратно в непрерывный сигнал. Однако, некоторые операции обработки звукового сигнала в таком представлении бывает не всегда удобной. Это связано с тем, что реально звуковой сигнал складывается из составляющих его частот с определённой амплитудой и фазой. Таким образом, применение таких операций обработки звукового сигнала как "фильтр нижних частот" или "фильтр верхних частот" требует преобразования представления звукового сигнала в виде отсчётов в его частотный спектр. После этого преобразования звуковой сигнал будет представлен в виде коэффициентов, соответствующим амплитудам и фазам частот, составляющих этот сигнал. Теперь, например, операция "фильтр нижних частот", которая "вырезает" из сигнала все частоты выше некоторой заданной, может просто обнулить коэффициенты, соответствующие частотам, которые необходимо "вырезать". На самом деле, круг областей применения такого преобразования значительно шире: обработка растровых изображений, телекоммуникации, исследование и измерение сигналов, радиолокация и т.д. Примером применения преобразования может служить передача данных в цифровой форме по аналоговым линиям телефонной сети (модем). Для передачи данных в цифровой форме, они сначала преобразуются в некоторый набор частот и передаются по линиям передач, а затем, на приёмной стороне выполняется обратное преобразование и восстанавливаются исходные данные. Далее рассмотрены некоторые фундаментальные понятия из математического анализа, необходимые для понимания работы преобразования Фурье. Ряд Фурье Рядом фурье называется бесконечная математическая последовательность, состоящая из коэффициентов при функциях синуса и косинуса вида: Доказано, что если некоторая периодическая функция с периодом 2l на интервале [-l, l] удовлетворяет условиям Дирихле (непрерывна и имеет конечное число экстремумов и точек разрыва I рода), то она может быть представлена в виде суммы ряда Фурье (разложена в ряд Фурье). Для определения коэффициентов ряда Фурье справедливы следующие формулы: Если раскладываемая функция является чётной ( f(-x) = f(x) ), то ряд Фурье состоит только из косинусов, т. е. все коэффициенты при синусах равны 0. Если раскладываемая функция является нечётной ( f(-x) = -f(x) ), то ряд Фурье состоит только из синусов, т. е. все коэффициенты при косинусах равны 0. В общем случае, коэффициенты при синусах и косинусах не равны 0. Таким образом, любую периодическую функцию, удовлетворяющую условиям Дирихле, можно разложить в ряд Фурье, тем самым представляя её в виде суммы синусов и косинусов. Преобразование Фурье в общем виде Предположим теперь, что исследуемая функция не является периодической, т. е. её период повторения равен бесконечности. В этом случае справедлива интегральная формула Фурье, получаемая путём предельного перехода из ряда Фурье периодической функции с периодом 2l при l : С этой формулой связаны так называемые интегральные преобразования Фурье: Из этих формул можно вывести синус- и косинус-преобразования Фурье соответственно для нечётных и чётных функций. Синус-преобразование Фурье: Косинус-преобразование Фурье: Таким образом, непрерывное преобразование Фурье позволяет представить непериодическую функцию в виде интеграла функции, представляющей в каждой своей точке коэффициент ряда Фурье для непериодической функции. Метод корреляций Метод корреляций позволяет определить тесноту линейной зависимости между исследуемой и базисной функциями. Это легче понять на примере. Пусть имеется импульсная радиолокационная станция. Для обнаружения цели эта станция формирует короткий высокочастотный радиоимпульс, огибающая которого имеет прямоугольную форму. Этот импульс излучается в пространство и отражается от цели. Так как в однородной среде электромагнитные волны распространяются с постоянной скоростью, близкой к скорости света в вакууме, то зная время, через которое отражённый от цели импульс поступил в приёмник радиолокационной станции, можно определить расстояние до цели. Однако здесь возникает следующая проблема: так как только часть импульса отражается от цели и поступает обратно в приёмник, а также в приёмник поступают некоторые помехи и сам приёмник имеет некоторый коэффициент шума, результирующий импульс будет несколько размыт на фоне помех и шумов. Возникает вопрос, как определить какая часть графика, представленного на рис. будет представлять отражённый импульс? Воспользуемся следующим методом: возьмём за основу функцию, которая на некотором интервале будет иметь скачок, представляющий собой идеальный отражённый импульс (при отсутствии помех и ослаблений) . Далее будем в каждой точке перемножать эту базисную функцию на функцию, формируемую приёмником, а затем проинтегрируем полученную функцию. Таким образом, это преобразование позволяет определить насколько велико значение базисной функции в исследуемой. Если исследуемая функция в каждой своей точке равна базисной, то интеграл функции, полученной в результате этого преобразования будет иметь максимальное значение. Если исследуемая функция ни как не отражает базисную, то результат будет равен 0. В промежуточных вариантах значения интеграла результирующей функции будут отражать то, насколько точно исследуемая функция соответствует базисной. Это и есть метод корреляций. Вернёмся к определению расстояния до цели. Теперь будем формировать новую базисную функцию, интервал скачка которой будем смещать вправо по оси абсцисс. Как только этот интервал будет соответствовать интервалу функции, соответствующему отражённому импульсу, значение интеграла результирующей функции будет максимальным. Таким образом может решаться задача об определении расстояния до цели. Дискретное преобразование Фурье В дискретном преобразовании Фурье исследуемая функция является периодической имеет конечный период повторения, и является дискретной. Фактически, дискретное преобразование Фурье позволяет представить дискретную функцию в виде конечного числа частот с определёнными значениями амплитуды и фазы (раскладывает функцию в её спектр). Это основывается на том, что по следствию из теоремы Котельникова в дискретном сигнале период, соответствующий наивысшей представимой частоте соответствует двум периодам дискретизации. Для определения амплитуд и фаз частотных составляющих сигнала, в дискретном преобразовании Фурье используется корреляция с базисными функциями синуса и косинуса. Спектр частот в дискретном преобразовании Фурье определяется из амплитуд синусов и косинусов, с частотами повторения в исследуемой выборке от 0 до N/2 раз, где N - количество элементов выборки. Преобразование Фурье раскладывает дискретизированный сигнал из N выборок на N/2 + 1 синусных и N/2 + 1 косинусных составляющих. Почему именно синусных и косинусных? Не только потому, что это следует из формулы ряда Фурье. Установлено, что результат сложения двух взаимно перпендикулярных колебаний с одинаковыми частотами, но разными амплитудами даёт колебание с той же частотой, но с другой амплитудой и фазой. Поэтому можно сказать, что дискретное преобразование Фурье раскладывает исследуемый сигнал по базисным функциям синуса и косинуса. Они являются аналогами двух взаимно перпендикулярных колебаний, так как по фазе смещены друг относительно друга на 90 градусов. Все вышеприведённые размышления приводят к следующим формулам дискретного преобразования Фурье: Эти формулы описывают прямое преобразование Фурье. ReX[x] - массив, содержащий значения косинусоидальных составляющих. ImX[x] - массив, содержащий значения синусоидальных составляющих. Такие обозначения введены в силу комплексного представления непрерывного преобразования Фурье (см. выше). При этом действительной части соответствуют косинусы, а мнимой - синусы. Это не должно вводить Вас в заблуждение - элементы этих массивов являются действительными числами и коэффициенты при синусах и косинусах являются действительными числами. Также, исследуемая функция является функцией действительного переменного. Комплексная форма преобразования Фурье может вводится для удобства записи двух интегралов - для косинуса и синуса. Массивы Re[x] и Im[x] составляют так называемый частотный домен (frequency domain), в то время как исходная выборка называется временным доменом (time domain). Свойства дискретного преобразования Фурье По приведённым выше формулам производится разложение исследуемого сигнала в его спектр. Выясним теперь свойства преобразования Фурье. Допустим, требуется произвести обратное преобразование - из частотных составляющих сформировать исходный сигнал. Для этого справедливы приведённые ниже формулы: Коэффициенты ReX[k] и ImX[k] определяются по следующим формулам: За исключением двух случаев: Такой процесс преобразования называется синтезом или обратным преобразованием Фурье. Заметим, что формулы обратного преобразования аналогичны формулам прямого преобразования, только теперь подынтегральной функцией являются коэффициенты при синусах и косинусах. Это свойство является очень важным и называется двойственностью преобразования Фурье. Свойство двойственности позволяет объяснить следующий факт: единичный импульс во временном домене (единичное значение одной выборки при нулевых значениях остальных) соответствует синусоиде и косинусоиде в частотном домене и наоборот (рис.). Во втором случае всё понятно - имеется один коэффициент при синусе или косинусе - это значит, что исходный сигнал (выборка) содержит составляющую одной частоты синусоидальной или косинусоидальной формы. Первый же случай может быть объяснён на основе двойственности преобразования Фурье. Описанный факт используется при построении алгоритма быстрого преобразования Фурье. Дело в том, что приведённые выше формулы для прямого и обратного преобразований имеют временную сложности реализующих их алгоритмов порядка O(n^2). Таким образом, при больших объёмах выборки, не удаётся за реальное время произвести преобразование Фурье. Для этой цели в середине 60-х годов был разработан алгоритм быстрого преобразования Фурье, который описывается в конце статьи. Свойства амплитуды и фазы Как уже было сказано, дискретное преобразование Фурье раскладывает исследуемый сигнал (выборку) на синусоидальные и косинусоидальные составляющие. В то же время, хотелось бы получить значения амплитуд и фаз частотных составляющих сигнала. Для перевода пары соответствующих коэффициентов при синусе и косинусе в амплитуду и фазу частотной составляющей справедливы следующие формулы: Обратное преобразование имеет вид: Полученные значения амплитуды и фазы называют полярным представлением (polar notation). Значения же коэффициентов при синусах и косинусах характеризуют прямоугольное представление (rectangular notation). При выполнении первого преобразования возможно несколько вариантов. 1. Значение коэффициента при косинусе равно 0. Тогда значение фазы равно ±/2. Соответствующее значение выбирается в соответствии со знаком коэффициента при синусе - если он положительный, то значение фазы /2, если отрицательный, то -/2. 2. Значение коэффициента при синусе равно 0. Тогда значение фазы равно 0 или . Соответствующее значение выбирается в соответствии со знаком коэффициента при косинусе - если он положительный, то значение фазы 0, если отрицательный, то . 3. "Дрейф" фазы. Подобный эффект наблюдается при очень малом, близком к 0, значении коэффициента при косинусе. При этом возможны резкие скачки значений фазы от -/2 к +/2. Так как значения коэффициента при косинусе очень малы, то, практически, фаза в этом случае не имеет значения и её можно принять равной -/2 или +/2. Как правило, значение фазы сигнала несёт в себе больше информации о форме и изменениях сигнала. Пусть, например, во временном домене имеется импульс некоторой длительности. Преобразуем этот временной домен в частотный и далее в полярное представление. Установим все амплитуды в случайное значение и преобразуем частотный домен обратно во временной. Таким образом, новый временной домен содержит в себе информацию, полученную из значений фаз частотных составляющих. Теперь в том месте, где был расположен импульс на его границах заметны скачки. Это и объясняет то, что фаза сигнала несёт больше информации о форме сигнала. Эффект Гиббса. Предположим, что из спектра прямоугольного импульса, убрали все частоты выше некоторой заданной. Далее применили обратное преобразование Фурье. Теперь на границах импульса будут отчётливо заметны затухающие колебания. Это объясняется тем, что границы идеального импульса представляются бесконечным числом синусоид и косинусоид. При увеличении этого количества от 0 до бесконечности форма фронтов импульса принимает свой истинный вид. При меньшем же количестве синусоид и косинусоид на границах импульса появляются искажения. Этот эффект носит название эффекта Гиббса. Таким образом, применение цифровых фильтров, построенных на основе преобразования Фурье в некоторых случаях может искажать сигнал. Этим ограничивается область применения таких фильтров. Алгоритм дискретного преобразования Фурье Ниже представлен исходный текст программы на C++, реализующей алгоритм преобразования Фурье. Мнимая функция GetValues() записывает во временной домен (массив f) некоторые значения, определяющие исследуемый сигнал. Быстрое преобразование Фурье Как было отмечено выше, обычное дискретное преобразование Фурье имеет большую временну сложность алгоритма (порядка O(n^2)), что ограничивает его применение. Алгоритм быстрого преобразования Фурье основывается на таком математическом преобразовании как свёртка функций. Свёртка функций подробно описывается в учебниках по мат. анализу и здесь не рассматривается. Остановимся на некоторых характерных чертах работы алгоритма быстрого преобразования Фурье. Количество элементов в выборке, поступающих на вход быстрого преобразования Фурье должно быть 2^n. Временная сложность преобразования - O(n · log n). Алгоритм быстрого комплексного преобразования Фурье работает следующим образом. 1. Для всей выборки производится так называемая побитовая обратная сортировка (bit-reversal sorting), в процессе которой меняются элементы с номерами, соответствующими прямой и перевёрнутой двоичной записи номера элемента. 2. Для каждого элемента отсортированной выборки выполняется одноэлементное дискретное преобразование Фурье. Это действие является мнимым, так как спектр одноточечного дискретного сигнала соответствует самому сигналу. Фактически, на этом этапе ничего не выполняется, но подразумевается, что выполняется. 3. Производится операция, обратная побитовой обратной сортировке. Но при этом преобразуются не элементы выборки, а элементы спектра сигнала. Причём на каждом шаге выполняется умножение частотного домена на комплексную синусоиду и сложение частотных доменов. 4. Результирующая выборка из 2^n комплексных элементов представляет спектр комплексного сигнала. Существуют модификации алгоритма для сигнала, выборки которого являются действительными числами. Для полного понимания работы этого алгоритма следует ознакомиться с теорией линейных систем, свёрткой, свойствами и парами преобразований Фурье и модификациями преобразований Фурье для функций комплексного переменного. КОНТРОЛЬНЫЕ ВОПРОСЫ 1.Дайте характеристику ряду Фурье 2.Охарактеризуйте дискретное преобразование Фурье 3.Приведите пример алгоритма дискретного преобразования Фурье 4.Опишите эффект Гиббса 5.Дайте характеристику метод корреляций при дискретном преобразовании ЛЕКЦИЯ 3. КОРРЕЛЯЦИОННАЯ ФУНКЦИЯ. ВЗАИМНАЯ КОРРЕЛЯЦИОННАЯ ФУНКЦИЯ. АВТОКОРРЕЛЯЦИОННАЯ ФУНКЦИЯ, ЕЕ СВОЙСТВА. Задача корреляционного анализа возникла из радиолокации, когда нужно было сравнить одинаковые сигналы, смещённые во времени. Для количественного определения степени отличия сигнала и его смещённой во времени копии принято вводить автокорреляционную функцию (АКФ) сигнала равную скалярному произведению сигнала и его сдвинутой копии. (4.1) Свойства АКФ 1) При автокорреляционная функция становится равной энергии сигнала: (4.2) 2) АКФ – функция чётная (4.3) 3) Важное свойство автокорреляционной функции состоит в следующем: при любом значении временного сдвига модуль АКФ не превосходит энергии сигнала: 4) Обычно, АКФ представляется симметричной линей с центральным максимумом, который всегда положителен. При этом в зависимости от вида сигнала автокорреляционная функция может иметь как монотонно убывающей, так и колеблющийся характер. Существует тесная связь между АКФ и энергетическим спектром сигнала. В соответствии с формулой (4.1) АКФ есть скалярное произведение . Здесь символом обозначена смещённая во времени копия сигнала . Обратившись к теореме Планшереля – можно записать равенство: (4.4) (4.4) Таким образом, приходим к результату (4.5) Квадрат модуля спектральной плотности представляет собой энергетический спектр сигнала. Итак энергетический спектр и автокорреляционная функция связаны парой преобразований Фурье. Ясно, что имеется и обратное соотношение (4.6) Эти результаты принципиально важны по двум причинам: во-первых, оказывается возможным оценивать корреляционные свойства сигналов, исходя из распределения их энергии по спектру. Во-вторых, формулы (4.5), (4.6) указывают путь экспериментального определения энергетического спектра. Часто удобнее вначале получить АКФ, а затем, используя преобразование Фурье, найти энергетический спектр сигнала. Такой приём получил распространение при исследовании свойств сигналов с помощью быстродействующих ЭВМ в реальном масштабе времени. Часто вводят удобный числовой параметр – интервал корреляции , представляющий собой оценку ширины основного лепестка АКФ. 9.. Взаимокорреляционная функция и ее свойства. Связь взаимокорреляционной функции и взаимного энергетического спектра. Взаимокорреляционной функцией (ВКФ) двух вещественных сигналов и называется скалярное произведение вида: (4.8) ВКФ служит мерой «устойчивости» ортогонального состояния при сдвигах сигналов во времени. При прохождении этих сигналов через различные устройства возможно, что сигнал будет сдвинут относительно сигнала на некоторое время . Свойства ВКФ. 1) В отличие от АКФ одиночного сигнала, ВКФ, описывающая свойства системы двух независимых сигналов, не является чётной функцией аргумента : (4.9) 2) Если рассматриваемые сигналы имеют конечные энергии, то их ВКФ ограничена. 3) При значения ВКФ вовсе не обязаны достигать максимума. Примером ВКФ может служить взаимокорреляционная функция прямоугольного и треугольного видеоимпульсов. На основании теоремы Планшереля получаем (4.11) Таким образом, взаимокорреляционная функция и взаимный энергетический спектр связаны между собой парой преобразований Фурье. Анализируя формулу обратного преобразования Фурье, приходим к выводу, что произвольный сигнал с известной спектральной плотностью можно записать как сумму двух составляющих, каждая из которых содержит или только положительные, или только отрицательные частоты: (3.1) Назовём функцию (3.2) аналитическим сигналом, отвечающим колебанию S(t). Итак, аналитический сигнал: (3.6) На комплексной плоскости этот сигнал отображается вектором, модуль и фазовый угол которого изменяются во времени. Проекция аналитического сигнала на вещественную ось в любой момент времени равна исходному сигналу . (3.8) Анализируя (3.7) и (3.8), можно убедиться в том, что спектральная плотность исходного и сопряжённого сигналов связаны между собой следующим образом: . (3.9) Чтобы на практике получить сопряжённый сигнал, необходимо исходное колебание подать на вход некоторой системы, которая осуществляет поворот фаз всех спектральных составляющих на угол -в области положительных частот и на угол в области отрицательных частот, не изменяя по амплитуде. Формула . (3.9) показывает, что спектральная плотность сопряжённого сигнала есть произведение спектра исходного сигнала и функции . В соответствии с обратной теоремой о свёртке сопряжённый сигнал представляет собой свёртку двух функций и которая является обратным преобразованием Фурье по отношении к функции . Таким образом, сопряжённый сигнал связан с исходным сигналом соотношением: (3.11) Можно поступить и по иному, выразив сигнал через , который полагается известным. Для этого достаточно заметить, что из (3.9) вытекает следующая связь между спектральными плотностями: Поэтому соответствующая формула будет отличаться от (3.11) лишь знаком: (3.12) Формулы (3.11) и (3.12) называются прямым и обратным преобразованием Гильберта. Символическая запись его такова: (3.13) Функция называется ядром этих преобразований. КОНТРОЛЬНЫЕ ВОПРОСЫ 1. Перечислите свойства АКФ 2. Опишите взаимокорреляционную функцию двух сигналов 3.Дайте определение понятию – энергетический спектр сигнала ЛЕКЦИЯ 4. СВЯЗЬ МЕЖДУ КОРРЕЛЯЦИОННЫМИ ФУНКЦИЯМИ И СПЕКТРАМИ СИГНАЛОВ. Из свойств преобразований Фурье известно, что спектр произведения двух сигналов  определяется соотношением свертки: В частном случае при  получается: Если положить теперь ,, То ,  В этом случае получаем для автокорреляционной функции следующее выражение Так как , то Величина S2() имеет смысл спектральной плотности энергии сигнала. При этом прямое преобразование Фурье можно записать в виде: Таким образом, спектральная плотность энергии сигнала связана с автокорреляционной функцией парой преобразований Фурье, что означает однозначную связь ширины частотного спектра и длительности АКФ: чем шире полоса частот, занимаемая сигналом, тем меньше интервал корреляции, т.е. сдвиг x, в пределах которого корреляционная функция отлична от нуля, и наоборот. Кроме того, поскольку АКФ не зависит от ФЧХ сигнала, а на форму функции существенно влияет её ФЧХ, то можно сказать, что различным по форме сигналам, обладающим одинаковыми АЧХ и разными ФЧХ, соответствуют одинаковые АКФ. Соотношение между корреляционной функцией и спектральной характеристикой сигнала Воспользуемся выражением (1.63), в котором положим f(t)=s(t),g(t)=s(t+τ) и соответственно F(ω) = S(ω), G(ω) = S(ω) e-iωτ.Тог­да получим Учитывая, что S(ω)S* (ω) =S2(ω), приходим к искомому соотно­шению (1.83) На основании известных свойств преобразований Фурье можно также написать: (1.84) Итак, прямое преобразование Фурье (1.84) корреляционной функции Bs(τ) дает спектральную плотность энергии , а преобразование (2.83) дает корреляционную функциюВs(τ). Из выражений (1.83) и (1.84) вытекают свойства: чем шире спектр S(ω)сигнала,тем меньше интервал кор­реляции, т. е. сдвигτ,в пределах которого корреляционная функция отлич­на от нуля. Соответственно чем больше интервал корреляции заданного сигнала, тем уже его спектр. Из выражений (1.83) и (1.84) также видно, что корреляционная функ­ция Bs(τ) не зависит от ФЧХ спектра сигнала. Так как при заданном ам­плитудном спектреS(ω) форма функцииs(t) существенно зависит от ФЧХ, то можно сделать следующее заключение:различным по форме сигналам s(t), обладающим одинаковыми амплитудными спектрами, соответствуют одинаковые корреляционные функции Bs (τ). КОНТРОЛЬНЫЕ ВОПРОСЫ 1. Что характеризует корреляционная функция? 2. Как определяется корреляционная функция? 3.Как связаны математически корреляционная функция и энергетический спектр сигнала. 4.Назовите основные свойства корреляционной функции. ЛЕКЦИЯ 5. КЛАСС ОРТОГОНАЛЬНЫХ ФУНКЦИЙ. ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ ГИЛЬБЕРТА. СВОЙСТВА ПРЕОБРАЗОВАНИЯ ГИЛЬБЕРТА Преобразование Гильберта для любого произвольного сигнала представляет собой идеальный широкополосный фазовращатель, который осуществляет поворот начальных фаз всех частотных составляющих сигнала на угол, равный 90о (сдвиг на /2). Применение преобразования Гильберта позволяет выполнять квадратурную модуляцию сигналов, в каждой текущей координате модулированных сигналов производить определение огибающей и мгновенной фазы и частоты сигналов, выполнять анализ каузальных систем обработки сигналов. Давид Гильберт (Hilbert, 1862-1943), немецкий математик. Окончил Кенигсбергский университет. В 1895-1930 годах профессор Геттингенского университета. Исследования Гильберта оказали большое влияние на развитие многих разделов математики, а его деятельность в Гёттингенском университете содействовала тому, что Геттинген являлся одним из основных мировых центров математической мысли. Определение преобразования. Прямое преобразование Гильберта произвольной действительной функции x(t), - < t < , результат которого будем отображать со знаком тильды над индексом исходной функции, задается сверткой x(t) с функцией hb(t) = 1/(t): (t) = ТН[x(t)] = x(t) * (1/t), (t) = . (17.1.1) Функция 1/(t-) называется ядром преобразования Гильберта. Обратное преобразование Гильберта определяется выражением: x(t) = - (17.1.1') Рис. 17.1.1. Интегралы преобразования имеет особую точку a = t-  0, в которой при вычислении используется их главное значение по Коши: [... +...]. Оператор Гильберта определен по аргументу от - до  и имеет полюс в точке t=0 с разрывом значений от - до . Основной участок формы оператора Гильберта и пример преобразования сигнала приведены на рис. 17.1.1. Функции x(t) и (t) обычно называют сопряженными по Гильберту. Спектральная характеристика преобразования. Выполним преобразование Фурье функции (17.1.1). В общей форме: (f) = TF[(t)] = X(f)  Hb(f), (17.1.2) (f) =(t) exp(-j2ft) dt. (17.1.2') Заметим, что произведение X(f)Hb(f) не является преобразованием Гильберта спектральной функции X(f). Это не более чем преобразование Фурье свертки функций: x(t)*hb(t)  X(f)Hb(f), которое позволяет вычислить результат преобразования Гильберта во временной области через частотную область: (t) =(f)exp(j2ft) df =X(f)Hb(f)exp(j2ft) df . Рис. 17.1.2. Функция hb(t)=1/t является нечетной, а спектр этой функции, представленный только мнимой частью, является обратной сигнатурной функцией (рис. 17.1.2): Hb(f) = TF[1/t] = -jsgn(f) = (17.1.3) Рис.17.1.3. Соответственно, формулы (17.1.1) задают преобразование сигнала x(t) системой, частотная передаточная характеристика которой отображается функцией -jsgn(f). Фурье-образ функции (t): (f) = -j sgn(f)X(f). (17.1.2") Изменение спектра сигналов при выполнении преобразования Гильберта. На рис. 17.1.3 приведено преобразование радиоимпульсного сигнала x(t) с несущей частотой fo в сигнал (t) во временной области непосредственно через операцию свертки по (17.1.1). Сигнал x(t) является односторонним каузальным. Спектр сигнала содержит реальную и мнимую составляющие, т.е. может быть записан в виде X() = Re(X()) + jIm(X()). Эти составляющие для сигнала x(t) показаны непрерывными кривыми на рис. 17.1.4. Рис. 17.1.4. При выполнении преобразования (17.1.2") реальная и мнимая части спектра X() умножаются на -jsgn(). Функция Re(X()) умножается на 1 при <0, на 0 при =0 и на –1 при >0, и тем самым превращается в нечетную мнимую часть Im(()) спектра () функции (t), показанную пунктиром. Это означает, что все косинусные гармоники сигнала, которым соответствует реальная часть спектра сигнала, превращаются в синусные гармоники. Аналогично на функцию -jsgn() умножается и мнимая функция jIm(X(, при этом сигнатурная функция инвертируется (-jj = 1), что меняет знак левой части функции Im(X( – области отрицательных частот, и превращает ее в реальную четную часть Re(()) спектра (). Синусные гармоники спектра сигнала превращаются в косинусные гармоники. Угол, на который изменяется фаза гармоник, можно определить из следующих соображений. Частотную характеристику Hb(f) (17.1.3) можно записать в следующем виде: Hb(f) = |Hb(f)|exp(jh(f)), где |Hb(f)| = 1. Hb(f) = -jsgn(f) = , (17.1.3') Если спектр функции x(t) также представить в виде X(f) = |X(f)|exp(jx(f)), то выражение (17.1.2) преобразуется к следующей форме: (f) = |X(f)|exp(jx(f))exp(jh(f)) = |X(f)|exp[j(x(f)+h(f))], (17.1.2''') т.е. амплитудный спектр сигнала (t) – как результат преобразования Гильберта сигнала x(t), не изменяется и остается равным амплитудному спектру сигнала x(t). Фазовый спектр сигнала (t) (начальные фазовые углы всех гармонических составляющих сигнала) сдвигается на -90о при f > 0 и на 90о при f < 0 относительно фазового спектра сигнала x(t). Но такой фазовый сдвиг и означает не что иное, как превращение косинусных гармоник в синусные, а синусных в косинусные. Последнее нетрудно проверить на единичной гармонике. Если x(t) = cos(2fot), то имеем следующее преобразование Гильберта через частотную область: (t) = H[x(t)]  TF[H[x(t)]] = -j sgn(f)[(f+fo)+(f-fo)]/2. (17.1.4) (f) = -j[-(f+fo)+(f-fo)]/2 = j·[(f+fo)-(f-fo)]/2. (17.1.5) Но последнее уравнение - спектр синусоиды. При обратном преобразовании Фурье: (t) = TF-1[(f)] = sin(2fot). (17.1.6) При x(t) = sin(2ft) аналогичная операция дает (t) = -cos(2fot). Знак минус демонстрирует отставание (запаздывание) выходного сигнала преобразования, как операции свертки, от входного сигнала. Для гармонических сигналов любой частоты с любой начальной фазой это запаздывание составляет четверть периода колебаний. На рис. 17.1.3 этот сдвиг на четверть периода для единичной гармонической составляющей (несущей частоты радиоимпульса) виден достаточно наглядно. Таким образом, преобразование Гильберта, по существу, представляет собой идеальный фазовращатель, осуществляющий фазовый сдвиг на 900 всех частотных составляющих сигналов одновременно. Сдвиг фазы спектров сигналов x(t) на /2 определяет изменение четности и самих сигналов: четный x(t)  нечетный (t), и наоборот. Спектры каузальных функций. Каузальная (физически осуществимая) линейная система (равно как и произвольная причинно обусловленная функция) задается односторонним импульсным откликом (выражением) h(t), t  0, и имеет частотную характеристику H(f): H(f) = A(f) - jB(f), где A(f) и B(f) - действительная (четная) и мнимая (нечетная) части частотной характеристики. Осуществим обратное преобразование Фурье для всех частей этого выражения: h(t) = a(t) + b(t), a(t) =A(f) cos(2ft) df, b(t) =B(f) sin(2ft) df, где a(t) и b(t) - соответственно четная и нечетная части импульсного отклика h(t). Условие каузальности для импульсного отклика (h(t) = 0 при t<0) будет выполнено, если при t<0 функции a(t) и b(t) компенсируют друг друга. Тогда общее условие каузальности, как можно наглядно видеть на рис. 17.1.5, с учетом нечетности функции b(t), запишется в следующем виде: b(t) = -a(t), b(t) = -h(t)/2, a(t) = h(t)/2, t < 0, (17.1.7) b(t) = 0, a(t) = a(0) = h(0), t = 0, b(t) = a(t) = h(t)/2, t > 0. Рис. 17.1.5. Параметры каузальной функции. Из этих условий следует, что нечетная функция b(t) в каузальной системе однозначно связана с четной функцией: b(t) = sgn(t)a(t), (17.1.8) Осуществляя преобразование Фурье обеих частей данного равенства при известном преобразовании сигнатурной функции (sgn(t)  j/f), получаем: Im(H(f)) = (j/f) * A(f), или, с учетом знака мнимой части: B(f) = -(1/f) * A(f) = -(1/) [A(v)/(f-v)] dv. (17.1.9) Аналогично определяется и действительная компонента спектра по мнимой части: A(f) = (1/f) * B(f) = (1/) [B(v)/(f-v)] dv. (17.1.10) Таким образом, реальная и мнимая части спектра физически осуществимых (односторонних) систем, а равно и произвольных каузальных сигналов, связаны парой преобразований Гильберта. Они позволяют производить определение любой, действительной или мнимой, части частотной характеристики каузальной функции путем свертки другой ее части с функцией 1/f. Для любых произвольных функций x(t) и y(t), имеющих Фурье – образы X(), Y() и преобразования Гильберта (t) = ТН[x(t)] и (t) = ТН[y(t)], действительны следующие свойства: Линейность. ТН[ax(t)+by(t)] = a(t)+b(t) при любых постоянных значениях коэффициентов а и b для любых произвольных функций x(t) и y(t). Сдвиг. ТН[x(t-a)] = (t-a). Преобразование константы, а в силу линейности преобразования, и постоянной составляющей сигнала, равно нулю. Это прямо следует из нечетности ядра преобразования Гильберта. Отсюда следует, что при преобразовании Гильберта из квадратурной составляющей исключается постоянная составляющая. Свойство четности и нечетности определяется сдвигом всех гармоник сигнала на /2, при этом четные сигналы x(t) дают нечетные сигналы (t), и наоборот. Это действительно и для произвольных сигналов относительно их четных и нечетных частей. Последовательное двойное преобразование Гильберта возвращает исходную функцию с обратным знаком ТН[ТН[x(t)]] = ТН[(t)] = -x(t). Это определяется тем, что при двойном преобразовании фазы всех гармоники сигнала сдвигаются на , что изменяет знак их гармоник. Однако в силу исключения из сигнала при первом преобразовании постоянной составляющей, при двойном преобразовании сигнал x(t) восстанавливается с исключенным средним значением по интервалу задания. Обратное преобразование Гильберта, по существу, это второе преобразование в последовательном двойном преобразовании Гильберта с изменением знака результата: x(t) = ТH-1[(t)] = -= (t) * (-1/t). (17.2.1) Альтернативная форма вычисления x(t) из (t): x(t) = TF-1[(j sgn(f)TF[(t)]]. (17.2.1') Подобие при изменении масштаба аргумента: ТН[x(at)] = (at). Энергетическая эквивалентность: x2(t) dt =2(t) dt. (17.2.2) Это следует из теоремы Парсеваля (энергия сигнала равна сумме энергии всех частотных составляющих сигнала) и равенства модулей спектров сигналов x(t) и (t) (энергия сигнала не зависит от его фазовочастотной характеристики). Свойство ортогональности: x(t)(t) dt = 0. (17.2.3) Если все косинусные составляющие сигнала x(t) превращаются в ортогональные им синусные составляющие сигнала , а синусные – в ортогональные им косинусные, то и сигналы x(t) и должны быть ортогональны. Из теоремы Парсеваля следует: x(t)(t) dt = X*(f)(f). Функция X*(f)(f) = -X*j sgn(f)X(f) = -j sgn(f)|X(f)|2 является нечетной, а поэтому определенный интеграл от этой функции по симметричным относительно нуля пределам равен нулю. Ортогональность сигналов наглядно видна на рис. 17.1.1. Свойство свертки: TH[x(t) * y(t)] = (t) * y(t) = x(t) * (t). (17.2.4) Это вытекает из следующих соображений. Примем z(t) = x(t) * y(t), при этом: Z(f) = X(f)Y(f), (f) = -j sgn(f)Z(f) = -j sgn(f) X(f)Y(f). (f) = [-j sgn(f) X(f)]Y(f) = (t)Y(f)  (t) * y(t). (f) = X(f)[-j sgn(f) Y(f)] = X(f)(f)  x(t) *(t). Отсутствие коммутативности с преобразованием Фурье: TF[ТН[x(t)]]  ТН[TF[x(t)]]. (17.2.5) Свойство модуляции: Модулирующие сигналы u(t), как правило, имеют ограниченный спектр, максимальные частоты которого  много меньше значения несущей частоты o, при этом: ТН[u(t)cos(ot)] = u(t)sin(ot). (17.2.6) Для четных функций u(t) это свойство очевидно. При переходе в частотную область: ТН[u(t)cos(o)]  -jsgn()[U() * ((+o)+(-o))]. Множитель -jsgn() является знаковой константой по  и может быть внесен под интеграл свертки и умножен на ((+o)+(-o)), что, как уже рассматривалось ранее (см. 17.1.4 – 17.1.6), при обратном преобразовании Фурье дает u(t)sin(ot). Аналогично можно показать, что ТН[u(t)sin(ot)] = -u(t)cos(ot). (17.2.7) Преобразование Гильберта аналоговых сигналов целесообразно выполнять не по формулам линейной свертки с оператором 1/t, который стремится к  при t  0, а через спектр аналитической функции: z(t) = x(t) + j(t)  X(f) + j(f) = Z(f). (17.3.1) Заменяя в этом выражении функцию (f) = -j sgn(f)X(f), получаем: Z(f) = [1+sgn(f)]X(f), (17.3.2) где функция 1+sgn(f) равна 0 при f < 0, 1 при f = 0 и 2 при f > 0, при этом: Z(f) = , (17.3.2') т.е. спектр функции z(t) является односторонним и устанавливается непосредственно по спектру функции x(t) при f  0 (см. также (17.1.13)). Обратное преобразование Фурье функции Z(f) должно давать комплексную функцию z(t), при этом из (17.3.2') следует: x(t) = Re [2X(f) exp(j2ft) df], (17.3.3) (t) = Im [2X(f) exp(j2ft) df]. (17.3.3') В дискретной форме, при общем числе N отсчетов функции x(t) с шагом t, с шагом по частоте f =1/(Nt): X(nf) = tx(kt)exp(-j2kn/N), n = 0,1,...,N/2. (17.3.4) х(kt) = fRe[Xo+2X(nf)exp(j2kn/N)]. (17.3.5') (kt) = 2fIm[X(nf)exp(j2kn/N)]. (17.3.5) Рис. 17.3.1. На рис. 17.3.1 приведен пример преобразования Гильберта, выполненный через частотную область. Естественно, что при реальном использовании преобразования Гильберта выполнять вычисления по (17.3.5') не требуется. Оператор дискретного преобразования Гильберта hb(kt)  1/t на интервале от -Т до Т с шагом t можно получить обратным преобразованием Фурье частотной характеристики Hb(f) (выражение 17.1.3) в интервале от -fN до fN (fN=1/2t). При t=1: hb(kt) =Hb(f) exp(j2fkt) df =j exp(j2fkt) df -j exp(j2fkt) df = = [1/(2kt)][1-exp(-jkt)-exp(jkt)+1] = [1/(kt)][1-(exp(-jkt)+exp(jkt)/2] = = [1/(kt)](1-cos(kt)) = [2/(kt)] sin2(kt/2). (17.3.6) hb(kt) = 2/(t), k = 1, 3, 5, ... , (17.3.6') hb(kt) = 0, k = 0, 2, 4, ... . Нетрудно убедиться, что коэффициент усиления постоянной составляющей оператора равен нулю, а коэффициент усиления дисперсии помех равен 1. В частотной области при выполнении преобразования Гильберта спектральных функций оператор свертки hb(kf)1/f не отличается от приведенного для временной области. Согласно методу преобразований Гильберта, огибающая и мгновенная частота сигнала жёстко связаны друг с другом и их нельзя выбрать произвольно. КОНТРОЛЬНЫЕ ВОПРОСЫ 1.Перечислите свойства преобразований Гильберта. 2.Охарактеризуйте преобразования Гильберта для узкополосного сигнала 3.Напишите формулу оператора дискретного преобразования Гильберта 4.Перечислите свойство свертки 5.Дайте характеристику свойству ортогональности ЛЕКЦИЯ 6. ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЕ. ВИДЫ И СРЕДСТВА ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ, ОРТОНОРМИРОВАННЫЕ БАЗИСЫ ВЕЙВЛЕТОВ. Вейвлетное преобразование сигналов является обобщением спектрального анализа, типичный представитель которого – классическое преобразование Фурье. Термин "вейвлет" (wavelet) в переводе с английского означает "маленькая (короткая) волна". Вейвлеты - это обобщенное название семейств математических функций определенной формы, которые локальны во времени и по частоте, и в которых все функции получаются из одной базовой (порождающей) посредством ее сдвигов и растяжений по оси времени. Вейвлет-преобразования рассматривают анализируемые временные функции в терминах колебаний, локализованных по времени и частоте. Как правило, вейвлет-преобразования (WT) подразделяют на дискретное (DWT) и непрерывное (CWT). DWT используется для преобразований и кодирования сигналов, CWT - для анализа сигналов. Вейвлет-преобразования в настоящее время принимаются на вооружение для огромного числа разнообразных применений, нередко заменяя обычное преобразование Фурье. Это наблюдается во многих областях, включая молекулярную динамику, квантовую механику, астрофизику, геофизику, оптику, компьютерную графику и обработку изображений, анализ ДНК, исследования белков, исследования климата, общую обработку сигналов и распознавание речи. Вейвлетный анализ представляет собой особый тип линейного преобразования сигналов и физических данных. Базис собственных функций, по которому проводится вейвлетное разложение сигналов, обладает многими специфическими свойствами и возможностями. Вейвлетные функции базиса позволяют сконцентрировать внимание на тех или иных локальных особенностях анализируемых процессов, которые не могут быть выявлены с помощью традиционных преобразований Фурье и Лапласа. К таким процессам в геофизике относятся поля различных физических параметров природных сред. В первую очередь это касается полей температуры, давления, профилей сейсмических трасс и других физических величин. Вейвлеты имеют вид коротких волновых пакетов с нулевым средним значением, локализованных по оси аргументов (независимых переменных), инвариантных к сдвигу и линейных к операции масштабирования (сжатия/растяжения). По локализации во временном и частотном представлении вейвлеты занимают промежуточное положение между гармоническими функциями, локализованными по частоте, и функцией Дирака, локализованной во времени. Теория вейвлетов не является фундаментальной физической теорией, но она дает удобный и эффективный инструмент для решения многих практических задач. Основная область применения вейвлетных преобразований – анализ и обработка сигналов и функций, нестационарных во времени или неоднородных в пространстве, когда результаты анализа должны содержать не только частотную характеристику сигнала (распределение энергии сигнала по частотным составляющим), но и сведения о локальных координатах, на которых проявляют себя те или иные группы частотных составляющих или на которых происходят быстрые изменения частотных составляющих сигнала. По сравнению с разложением сигналов на ряды Фурье вейвлеты способны с гораздо более высокой точностью представлять локальные особенности сигналов, вплоть до разрывов 1-го рода (скачков). В отличие от преобразований Фурье, вейвлет-преобразование одномерных сигналов обеспечивает двумерную развертку, при этом частота и координата рассматриваются как независимые переменные, что дает возможность анализа сигналов сразу в двух пространствах. Одна из главных и особенно плодотворных идей вейвлетного представления сигналов на различных уровнях декомпозиции (разложения) заключается в разделении функций приближения к сигналу на две группы: аппроксимирующую - грубую, с достаточно медленной временной динамикой изменений, и детализирующую - с локальной и быстрой динамикой изменений на фоне плавной динамики, с последующим их дроблением и детализацией на других уровнях декомпозиции сигналов. Это возможно как во временной, так и в частотной областях представления сигналов вейвлетами. Историческая справка. История спектрального анализа восходит к И. Бернулли, Эйлеру и Фурье, который впервые построил теорию разложения функций в тригонометрические ряды. Однако это разложение долгое время применялось как математический прием и не связывалось с какими-либо физическими понятиями. Спектральные представления применялись и развивались сравнительно узким кругом физиков–теоретиков. Однако, начиная с 20-х годов прошлого века, в связи с бурным развитием радиотехники и акустики, спектральные разложения приобрели физический смысл и практическое применение. Основным средством анализа реальных физических процессов стал гармонический анализ, а математической основой анализа - преобразование Фурье. Преобразование Фурье разлагает произвольный процесс на элементарные гармонические колебания с различными частотами, а все необходимые свойства и формулы выражаются с помощью одной базисной функции exp(jt) или двух действительных функций sin(t) и cos(t). Гармонические колебания имеют широкое распространение в природе, и поэтому смысл преобразования Фурье интуитивно понятен независимо от математической аналитики. Преобразование Фурье обладает рядом замечательных свойств. Областью определения преобразования является пространство L2 интегрируемых с квадратом функций, и многие физические процессы в природе можно считать функциями, принадлежащими этому пространству. Для применения преобразования разработаны эффективные вычислительные процедуры типа быстрого преобразования Фурье (БПФ). Эти процедуры входят в состав всех пакетов прикладных математических программ и реализованы аппаратно в процессорах обработки сигналов. Было также установлено, что функции можно разложить не только по синусам и косинусам, но и по другим ортогональным базисным системам, например, полиномам Лежандра и Чебышева, функциям Лагерра и Эрмита. Однако практическое применение они получили только в последние десятилетия ХХ века благодаря развитию вычислительной техники и методов синтеза цифровых линейных систем обработки данных. Непосредственно для целей спектрального анализа подобные ортогональные функции не нашли широкого применения из-за трудностей интерпретации получаемых результатов. По тем же причинам не получили развития в спектральном анализе функции типа "прямоугольной волны" Уолша, Радемахера, и пр. Теоретические исследования базисных систем привели к созданию теории обобщенного спектрального анализа, которая позволила оценить пределы практического применения спектрального анализа Фурье, создала методы и критерии синтеза ортогональных базисных систем. Иллюстрацией этому является активно развивающаяся с начала 80-х годов прошлого столетия теория базисных функций типа вейвлет. Благодаря прозрачности физической интерпретации результатов анализа, сходной с "частотным" подходом в преобразовании Фурье, ортогональный базис вейвлетов стал популярным и эффективным средством анализа сигналов и изображений в акустике, сейсмике, медицине и других областях науки и техники. Вейвлет-анализ является разновидностью спектрального анализа, в котором роль простых колебаний играют функции особого рода, называемые вейвлетами. Базисная функция вейвлет – это некоторое "короткое" колебание, но не только. Понятие частоты спектрального анализа здесь заменено масштабом, и, чтобы перекрыть "короткими волнами" всю временную ось, введен сдвиг функций во времени. Базис вейвлетов – это функции типа ((t-b)/a), где b - сдвиг, а – масштаб. Функция (t) должна иметь нулевую площадь и, еще лучше, равными нулю первый, второй и прочие моменты. Фурье-преобразование таких функций равно нулю при =0 и имеет вид полосового фильтра. При различных значениях масштабного параметра 'a' это будет набор полосовых фильтров. Семейства вейвлетов во временной или частотной области используются для представления сигналов и функций в виде суперпозиций вейвлетов на разных масштабных уровнях декомпозиции (разложения) сигналов. Первое упоминание о подобных функциях (которые вейвлетами не назывались) появилось в работах Хаара (Haar) еще в начале прошлого века. Вейвлет Хаара - это короткое прямоугольное колебание на интервале [0,1], показанное на рис. 1.1.1. Однако он интересен больше теоретически, так как не является непрерывно дифференцируемой функцией и имеет длинные "хвосты" в частотной области. В 30-е годы физик Paul Levy, исследуя броуновское движение, обнаружил, что базис Хаара лучше, чем базис Фурье, подходит для изучения деталей броуновского движения. Сам термин "вейвлет", как понятие, ввели в своей статье J. Morlet и A. Grossman, опубликованной в 1984 г. Они занимались исследованиями сейсмических сигналов с помощью базиса, который и назвали вейвлетом. Весомый вклад в теорию вейвлетов внесли Гуппилауд, Гроссман и Морлет, сформулировавшие основы CWT, Ингрид Добеши, разработавшая ортогональные вейвлеты (1988), Натали Делпрат, создавшая время-частотную интерпретацию CWT (1991), и многие другие. Математическая формализация в работах Mallat и Meyer привела к созданию теоретических основ вейвлет-анализа, названного мультиразрешающим (кратномасштабным) анализом. В настоящее время специальные пакеты расширений по вейвлетам присутствуют в основных системах компьютерной математики (Matlab, Mathematica, Mathcad, и др.), а вейвлет-преобразования и вейвлетный анализ используются во многих областях науки и техники для самых различных задач. Многие исследователи называют вейвлет-анализ "математическим микроскопом" для точного изучения внутреннего состава и структур неоднородных сигналов и функций. Не следует рассматривать вейвлет-методы обработки и анализа сигналов в качестве новой универсальной технологии решения любых задач. Возможности вейвлетов еще не раскрыты полностью, однако это не означает, что их развитие приведет к полной замене традиционных средств обработки и анализа информации, хорошо отработанных и проверенных временем. Вейвлеты позволяют расширить инструментальную базу информационных технологий обработки данных. Преобразование Фурье (ПФ). В основе спектрального анализа сигналов лежит интегральное преобразование и ряды Фурье. Напомним некоторые математические определения. В пространстве функций, заданных на конечном интервале (0,T), норма, как числовая характеристика произвольной функции s(t), вычисляется как корень квадратный из скалярного произведения функции. Для комплексных функций, квадрат нормы (энергия сигнала) соответствует выражению: ||s(t)||2 = s(t), s(t) =s(t) s*(t) dt, (1.1.1) где s*(t) – функция, комплексно сопряженная с s(t). Если норма функции имеет конечное значение (интеграл сходится), то говорят, что функция принадлежит пространству функций L2[R], R=[0,T], интегрируемых с квадратом (пространство Гильберта), и имеет конечную энергию. В пространстве Гильберта на основе совокупности ортогональных функций с нулевым скалярным произведением v(t), w(t) =v(t) w*(t) dt = 0 может быть создана система ортонормированных "осей" (базис пространства), при этом любой сигнал, принадлежащий этому пространству, может быть представлен в виде весовой суммы проекций сигнала на эти "оси" – базисных векторов. Значения проекций определяются скалярными произведениями сигнала с соответствующими функциями базисных "осей". Базис пространства может быть образован любой ортогональной системой функций. Наибольшее применение в спектральном анализе получила система комплексных экспоненциальных функций. Проекции сигнала на данный базис определяются выражением: Sn = (1/T)s(t) exp(-jnt) dt, n  (-∞, ∞), (1.1.2) где =2/T – частотный аргумент векторов. При известных выражениях базисных функций сигнал s(t) однозначно определяется совокупностью коэффициентов Sn и может быть абсолютно точно восстановлен (реконструирован) по этим коэффициентам: s(t) =Sn exp(jnt). (1.1.3) Уравнения (1.1.2) и (1.1.3) называют прямым и обратным преобразованием Фурье сигнала s(t). Любая функция гильбертова пространства может быть представлена в виде комплексного ряда Фурье (1.1.3), который называют спектром сигнала или его Фурье-образом. Ряд Фурье ограничивается определенным количеством членов N, что означает аппроксимацию с определенной погрешностью бесконечномерного сигнала N – мерной системой базисных функций спектра сигнала. Ряд Фурье равномерно сходится к s(t): ||s(t) -Sn exp(jnt)|| = 0. (1.1.4) Таким образом, ряд Фурье - это разложение сигнала s(t) по базису пространства L2(0,T) ортонормированных гармонических функций exp(jnt) с изменением частоты, кратным частоте первой гармоники 1=. Отсюда следует, что ортонормированный базис пространства L2(0,T) построен из одной функции exp(jt) = cos(t)+j·sin(t) с помощью масштабного преобразования независимой переменной. Для коэффициентов ряда Фурье справедливо равенство Парсеваля сохранения энергии сигнала в различных представлениях: (1/T)|s(t)|2 dt =|Sn|2. (1.1.5) Разложение в ряд Фурье произвольной функции y(t) корректно, если функция y(t) принадлежит этому же пространству L2(0,T), т.е. квадратично интегрируема с конечной энергией: |y(t)|2 dt <  , t  (0,T), (1.1.6) при этом она может быть периодически расширена и определена на всей временной оси пространства R(-, ) так, что y(t) = y(t-T), t  R, при условии сохранения конечности энергии в пространстве R(-, ). С позиций анализа произвольных сигналов и функций в частотной области и точного восстановления после преобразований можно отметить ряд недостатков разложения сигналов в ряды Фурье, которые привели к появлению оконного преобразования Фурье и стимулировали развитие вейвлетного преобразования. Основные из них: Ограниченная информативность анализа нестационарных сигналов и практически полное отсутствие возможностей анализа их особенностей (сингулярностей), т.к. в частотной области происходит «размазывание» особенностей сигналов (разрывов, ступенек, пиков и т.п.) по всему частотному диапазону спектра. Гармонические базисные функции разложения не способны отображать перепады сигналов с бесконечной крутизной типа прямоугольных импульсов, т.к. для этого требуется бесконечно большое число членов ряда. При ограничении числа членов ряда Фурье в окрестностях скачков и разрывов при восстановлении сигнала возникают осцилляции (явление Гиббса). Преобразование Фурье отображает глобальные сведения о частотах исследуемого сигнала и не дает представления о локальных свойствах сигнала при быстрых временных изменениях его спектрального состава. Так, например, преобразование Фурье не различает стационарный сигнал с суммой двух синусоид от нестационарного сигнала с двумя последовательно следующими синусоидами с теми же частотами, т.к. спектральные коэффициенты (1.1.2) вычисляются интегрированием по всему интервалу задания сигнала. Преобразование Фурье не имеет возможности анализировать частотные характеристики сигнала в произвольные моменты времени. Оконное преобразование Фурье. Частичным выходом из этой ситуации является оконное преобразование Фурье с движущейся по сигналу оконной функцией, имеющей компактный носитель. Временной интервал сигнала разделяется на подинтервалы и преобразование выполняется последовательно для каждого подинтервала в отдельности. Тем самым осуществляется переход к частотно-временному (частотно-координатному) представлению сигналов, при этом в пределах каждого подинтервала сигнал "считается" стационарным. Результатом оконного преобразования является семейство спектров, которым отображается изменение спектра сигнала по интервалам сдвига окна преобразования. Это позволяет выделять на координатной оси и анализировать особенности нестационарных сигналов. Размер носителя оконной функции w(t) обычно устанавливается соизмеримым с интервалом стационарности сигнала. По существу, таким преобразованием один нелокализованный базис разбивается на определенное количество базисов, локализованных в пределах функции w(t), что позволяет представлять результат преобразования в виде функции двух переменных - частоты и временного положения окна. Оконное преобразование выполняется в соответствии с выражением: S(,bk) = s(t) w*(t-bk) exp(-jt) dt. (1.1.7) Функция w*(t-b) представляет собой функцию окна сдвига преобразования по координате t, где параметром b задаются фиксированные значения сдвига. При сдвиге окон с равномерным шагом значения bk принимаются равными kb. В качестве окна преобразования может использоваться как простейшее прямоугольное окно, так и специальные весовые окна (Бартлетта, Гаусса, и пр.), обеспечивающие малые искажения спектра при вырезке оконных отрезков сигналов (нейтрализация явления Гиббса). Рис. 1.1.2. Пример оконного преобразования для нестационарного сигнала на большом уровне шума приведен на рис. 1.1.2. По спектру сигнала можно судить о наличии в его составе гармонических колебаний на трех частотах, определять соотношение между амплитудами этих колебаний и конкретизировать локальность колебаний по интервалу сигнала. Координатная разрешающая способность оконного преобразования определяется шириной оконной функции и обратно пропорциональна частотной разрешающей способности. При ширине оконной функции, равной b, частотная разрешающая способность определяется значением  = 2/b. При требуемой величине частотного разрешения  соответственно ширина оконной функции должна быть равна b = 2/. Для оконного преобразования Фурье эти ограничения являются принципиальными. Так, для рис. 1.1.2 при размере массива данных N = 300 и ширине оконной функции b = 100 частотная разрешающая способность результатов преобразования уменьшается в N/b = 3 раза по сравнению с исходными данными, и графики Sw(nSw) по координате n для наглядного сопоставления с графиком S(nS построены с шагом по частоте Sw = 3S, т.е. по точкам n = 0, 3, 6, … , N. Частотно-временное оконное преобразование применяется для анализа нестационарных сигналов, если их частотный состав изменяется во времени. Функция оконного преобразования (23.1.1) может быть переведена в вариант с независимыми переменными и по времени, и по частоте: S(t,) = s(t-) w() exp(-j) d. (23.1.2) Координатная разрешающая способность оконного преобразования определяется шириной оконной функции и, в силу принципа неопределенности Гейзенберга, обратно пропорциональна частотной разрешающей способности. Хорошая разрешающая способность по времени подразумевает небольшое окно времени, которому соответствует плохая частотная разрешающая способность и наоборот. Оптимальным считается ОПФ с гауссовым окном, которое получило название преобразование Габора (Gabor). Пример преобразования приведен на рис. 23.1.2 (в дискретном варианте вычислений. Рис. 23.1.2. На рис. 23.1.3 приведен пример вычисления и представления (модуль правой части главного диапазона спектра) результатов спектрограммы при дискретном задании зашумленного входного сигнала sq(n). Сигнал представляет собой сумму трех последовательных радиоимпульсов с разными частотами без пауз, с отношением сигнал/шум, близким к 1. Оконная функция wi задана в одностороннем варианте с эффективной шириной окна b  34 и полным размером М = 50. Установленный для результатов шаг по частоте  = 0.1 несколько выше фактической разрешающей способности 2/M = 0.126. Для обеспечения работы оконной функции по всему интервалу сигнала задавались начальные и конечные условия вычислений (продление на M точек обоих концов сигнала нулевыми значениями). Рис. 23.1.3. Как видно из приведенных примеров, оконное преобразование позволяет выделить информативные особенности сигнала и по времени, и по частоте. Разрешающая способность локализации по координатам и по частоте определяется принципом неопределенности Гейзенберга. В силу этого принципа невозможно получить произвольно точное частотно-временное представление сигнала. На рис. 23.1.4 приведен пример частотно-временного оконного преобразования сигнала, состоящего из 4-х непересекающихся интервалов, в каждом из которых сумма двух гармоник разной частоты. В качестве окна применена гауссова функция разной ширины. Узкое окно обеспечивает лучшее временное разрешение и четкую фиксацию границ интервалов, но широкие пики частот в пределах интервалов. Широкое окно напротив – четко отмечает частоты интервалов, но с перекрытием границ временных интервалов. При решении практических задач приходится выбирать окно для анализа всего сигнала, тогда как разные его участки могут требовать применения разных окон. Рис. 23.1.4. На рис. 23.1.3 приведен пример частотно-временного оконного преобразования сигнала, состоящего из 4-х непересекающихся интервалов, в каждом из которых сумма двух гармоник разной частоты. В качестве окна применена гауссова функция разной ширины. Узкое окно обеспечивает лучшее временное разрешение и четкую фиксацию границ интервалов, но широкие пики частот в пределах интервалов. Широкое окно напротив – четко отмечает частоты интервалов, но с перекрытием границ временных интервалов. При решении практических задач приходится выбирать окно для анализа всего сигнала, тогда как разные его участки могут требовать применения разных окон. Если сигнал состоит из далеко отстоящих друг от друга частотных компонент, то можно пожертвовать спектральным разрешением в пользу временного, и наоборот. Функции оконного спектрального анализа в Mathcad находятся в пакете Signal Processing. Они позволяют разбивать сигнал на поддиапазоны (с перекрытием или без перекрытия) и выполнять следующие операции: cspectrum(x,n,r[,w]) – расчет кросс-спектра сигнала х; pspectrum(x,n,r[,w]) – расчет распределения спектральной мощности сигнала; coherence(x,y,n,r[,w]) – расчет когерентности сигналов х и у; snr(x,y,n,r[,w]) – расчет отношения сигнал/шум для векторов х и у. Здесь: х и у – вещественные или комплексные массивы данных (векторы), n – число поддиапазонов разбиения входного сигнала х (от 1 до N – размера массива), к – фактор перекрытия поддиапазонов (от 0 до 1), w - код окна (1- прямоугольное, 2- трапеция, 3- треугольное, 4- окно Хеннинга, 5- окно Хемминга, 6- окно Блекмана). Принцип вейвлет-преобразования. Гармонические базисные функции преобразования Фурье предельно локализованы в частотной области (до импульсных функций Дирака при Т  ) и не локализованы во временной (определены во всем временном интервале от - до ). Их противоположностью являются импульсные базисные функции типа импульсов Кронекера, которые предельно локализованы во временной области и "размыты" по всему частотному диапазону. Вейвлеты по локализации в этих двух представлениях можно рассматривать как функции, занимающие промежуточное положение между гармоническими и импульсными функциями. Они должны быть локализованными как во временной, так и в частотной области представления. Однако при проектировании таких функций мы неминуемо столкнемся с принципом неопределенности, связывающим эффективные значения длительности функций и ширины их спектра. Чем точнее мы будем осуществлять локализацию временного положения функции, тем шире будет становиться ее спектр, и наоборот, что наглядно видно на рис. 1.1.5. Отличительной особенностью вейвлет-анализа является то, что в нем можно использовать семейства функций, реализующих различные варианты соотношения неопределенности. Соответственно, исследователь имеет возможность гибкого выбора между ними и применения тех вейвлетных функций, которые наиболее эффективно решают поставленные задачи. Вейвлетный базис пространства L2(R), R(-, ), целесообразно конструировать из финитных функций, принадлежащих этому же пространству, которые должны стремиться к нулю на бесконечности. Чем быстрее эти функции стремятся к нулю, тем удобнее использовать их в качестве базиса преобразования при анализе реальных сигналов. Допустим, что такой функцией является psi - функция t, равная нулю за пределами некоторого конечного интервала и имеющая нулевое среднее значение по интервалу задания. Последнее необходимо для задания локализации спектра вейвлета в частотной области. На основе этой функции сконструируем базис в пространстве L2(R) с помощью масштабных преобразований независимой переменной. Функция изменения частотной независимой переменной в спектральном представлении сигналов отображается во временном представлении растяжением/сжатием сигнала. Для вейвлетного базиса это можно выполнить функцией типа (t) => (amt), a = const, m = 0, 1, … , M, т.е. путем линейной операции растяжения/сжатия, обеспечивающей самоподобие функции на разных масштабах представления. Однако локальность функции (t) на временной оси требует дополнительной независимой переменной последовательных сдвигов функции (t) вдоль оси, типа (t) => (t+k), для перекрытия всей числовой оси пространства R(-, ). C учетом обеих условий одновременно структура базисной функции может быть принята следующей: (t) => (amt+k). (1.1.10) Для упрощения дальнейших выкладок значения переменных m и k примем целочисленными. При приведении функции (1.1.10) к единичной норме, получаем: mk(t) = am/2 (amt+k). (1.1.11) Если для семейства функций mk(t) выполняется условие ортогональности: nk(t), lm(t) =nk(t)·*lm(t) dt =nl·km, (1.1.12) то семейство mk(t) можно использовать в качестве ортонормированного базиса пространства L2(R). Произвольную функцию этого пространства можно разложить в ряд по базису mk(t): s(t) =Smk mk(t), (1.1.13) где коэффициенты Smk – проекции сигнала на новый ортогональный базис функций, как и в преобразовании Фурье, определяются скалярным произведением Smk = s(t), mk(t) =s(t)mk(t) dt, (1.1.14) при этом ряд равномерно сходиться: ||s(t) –Smk mk(t),|| = 0. При выполнении этих условий базисная функция преобразования (t) называется ортогональным вейвлетом. Простейшим примером ортогональной системы функций такого типа являются функции Хаара. Базисная функция Хаара определяется соотношением (t) = ( 1.1.15) Легко проверить, что при а = 2, m = 0, 1, 2, ..., k = 0, 1,2, … две любые функции, полученные с помощью этого базисного вейвлета путем масштабных преобразований и переносов, имеют единичную норму и ортогональны. На рис. 1.1.6 приведены примеры функций для первых трех значений m и b при различных их комбинациях, где ортогональность функций видна наглядно. Рис. 1.1.6. Функции Хаара. Вейвлетный спектр, в отличие от преобразования Фурье, является двумерным и определяет двумерную поверхность в пространстве переменных m и k. При графическом представлении параметр растяжения/сжатия спектра m откладывается по оси абсцисс, параметр локализации k по оси ординат – оси независимой переменной сигнала. Математику процесса вейвлетного разложения сигнала в упрощенной форме рассмотрим на примере разложения сигнала s(t) вейвлетом Хаара с тремя последовательными по масштабу m вейвлетными функциями с параметром а=2, при этом сам сигнал s(t) образуем суммированием этих же вейвлетных функций с одинаковой амплитудой с разным сдвигом от нуля, как это показано на рис. 1.1.7. Рис. 1.1.7. Скалярные произведения сигнала с вейвлетами. Для начального значения масштабного коэффициента сжатия m определяется функция вейвлета (1(t) на рис. 1.1.7), и вычисляется скалярное произведение сигнала с вейвлетом 1(t), s(t+k) с аргументом по сдвигу k. Для наглядности результаты вычисления скалярных произведений на рис. 1.1.7 построены по центрам вейвлетных функций (т.е. по аргументу k от нуля со сдвигом на половину длины вейвлетной функции). Как и следовало ожидать, максимальные значения скалярного произведения отмечаются там, где локализована эта же вейвлетная функция. После построения первой масштабной строки разложения, меняется масштаб вейвлетной функции (2 на рис. 1.1.7) и выполняется вычисление второй масштабной строки спектра, и т.д. Как видно на рис. 1.1.7, чем точнее локальная особенность сигнала совпадает с соответствующей функцией вейвлета, тем эффективнее выделение этой особенности на соответствующей масштабной строке вейвлетного спектра. Можно видеть, что для сильно сжатого вейвлета Хаара характерной хорошо выделяемой локальной особенностью является скачок сигнала, причем выделяется не только скачок функции, но и направление скачка. На рис. 1.1.8 приведен пример графического отображения вейвлетной поверхности реального физического процесса /4/. Вид поверхности определяет изменения во времени спектральных компонент различного масштаба и называется частотно-временным спектром. Поверхность изображается на рисунках, как правило, в виде изолиний или условными цветами. Для расширения диапазона масштабов может применяться логарифмическая шкала. Рис. 1.1.8. Пример вейвлетного преобразования. 1.2. основы Вейвлет - преобразования /1, 3, 7, 9, 11/. В основе вейвлет-преобразований, в общем случае, лежит использование двух непрерывных, взаимозависимых и интегрируемых по независимой переменной функций: Вейвлет-функции (t), как psi-функции времени с нулевым значением интеграла и частотным фурье-образом (ω). Этой функцией, которую обычно и называют вейвлетом, выделяются локальные особенности сигнала. В качестве вейвлетов обычно выбираются функции, хорошо локализованные и во временной, и в частотной области. Пример временного и частотного образа функции приведен на рис. 1.2.1. Масштабирующей функции (t), как временной скейлинг-функции phi с единичным значением интеграла, которой выполняется грубое приближение (аппроксимация) сигнала. Рис. 1.2.1. Вейвлетные функции в двух масштабах. Phi-функции присущи не всем, а, как правило, только ортогональным вейвлетам. Они необходимы для преобразования нецентрированных и достаточно протяженных сигналов при раздельном анализе низкочастотных и высокочастотных составляющих. Роль и использование phi-функции рассмотрим несколько позже. Непрерывное вейвлет-преобразование (НВП, CWT- Continious Wavelet Transform). Допустим, что мы имеем функции s(t) с конечной энергией в пространстве L2(R), определенные по всей действительной оси R(-, ). Для финитных сигналов с конечной энергией средние значения сигналов должны стремиться к нулю на ±. Непрерывным вейвлет-преобразованием (или вейвлетным образом) функции s(t)  L2(R) называют функцию двух переменных: С(a,b) = s(t), (a,b,t) = s(t)(а,b,t) dt, a, b  R, a ≠ 0. (1.2.1) где вейвлеты (a,b,t)  ab(t) – масштабированные и сдвинутые копии порождающего вейвлета (t)  L2(R), совокупность которых создает базис пространства L2(R). Порождающими функциями могут быть самые различные функции с компактным носителем - ограниченные по времени и местоположению на временной оси, и имеющие спектральный образ, локализованный на частотной оси. Базис пространства L2(R) целесообразно конструировать из одной порождающей функции, норма которой должна быть равна 1. Для перекрытия функцией вейвлета всей временной оси пространства используется операция сдвига (смещения по временной оси): (b,t) = (t-b), где значение b для НВП является величиной непрерывной. Для перекрытия всего частотного диапазона пространства L2(R) используется операция временного масштабирования вейвлета с непрерывным изменением независимой переменной: (a,t) = |а|-1/2(t/а). На рис. 1.2.1. видно, что если временной образ вейвлета будет расширяться (изменением значения параметра 'а'), то его "средняя частота" будет понижаться, а частотный образ (частотная локализация) перемещаться на более низкие частоты. Таким образом, путем сдвига по независимой переменной (t-b) вейвлет имеет возможность перемещаться по всей числовой оси произвольного сигнала, а путем изменения масштабной переменной 'а' (в фиксированной точке (t-b) оси) "просматривать" частотный спектр сигнала по определенному интервалу окрестностей этой точки. С использованием этих операций вейвлетный базис функционального пространства образуется путем масштабных преобразований и сдвигов порождающего вейвлета(t): (a,b,t) = |а|-1/2[(t-b)/а], a, b  R, a ≠ 0, (t) Î L2(R). (1.2.2) Нетрудно убедиться, что нормы вейвлетов (a,b,t) равны норме (t), что обеспечивает нормировочный множитель |а|-1/2. При нормировке к 1 порождающего вейвлета (t) все семейство вейвлетов также будет нормированным. Если при этом выполняется требование ортогональности функций, то функции (a,b,t) образуют ортонормированный базис пространства L2(R). Понятие масштаба ВП имеет аналогию с масштабом географических карт. Большие значения масштаба соответствуют глобальному представлению сигнала, а низкие значения масштаба позволяют различить детали. В терминах частоты низкие частоты соответствуют глобальной информации о сигнале, а высокие частоты - детальной информации и особенностям, которые имеют малую протяженность, т.е. масштаб вейвлета, как единица шкалы частотно-временного представления сигналов, обратен частоте. Масштабирование, как математическая операция, расширяет или сжимает сигнал. Большие значения масштабов соответствуют расширениям сигнала, а малые значения - сжатым версиям. В определении вейвлета коэффициент масштаба а стоит в знаменателе. Соответственно, а > 1 расширяет сигнал, а < 1 сжимает его. Процедура преобразования стартует с масштаба а=1 и продолжается при увеличивающихся значениях а, т.e. анализ начинается с высоких частот и проводится в сторону низких частот. Первое значение 'а' соответствует наиболее сжатому вейвлету. При увеличении значения 'а' вейвлет расширяется. Вейвлет помещается в начало сигнала (t=0), перемножается с сигналом, интегрируется на интервале своего задания и нормализуется на 1/. Результат вычисления С(a,b) помещается в точку (a=1, b=0) масштабно-временного спектра преобразования. Сдвиг b может рассматриваться как время с момента t=0, при этом координатная ось b повторяет временную ось сигнала. Для полного включения в обработку всех точек входного сигнала требуется задание начальных и конечных условий преобразования (определенных значений входного сигнала при t<0 и t>tmax на полуширину окна вейвлета). При одностороннем задании вейвлетов результат относится, как правило, к временному положению средней точки окна вейвлета. Затем вейвлет масштаба а=1 сдвигается вправо на значение b и процедура повторяется. Получаем значение, соответствующее t=b в строке а=1 на частотно-временном плане. Процедура повторяется до тех пор, пока вейвлет не достигнет конца сигнала. Таким образом получаем строку точек на масштабно-временном плане для масштаба а=1. Для вычисления следующей масштабной строки значение а увеличивается на некоторое значение. При НВП в аналитической форме b0 и a0. При выполнении преобразования в компьютере выполняется увеличение обоих параметров с определенным шагом. Тем самым осуществляется дискретизация масштабно-временной плоскости. Начальное значение масштабного коэффициента может быть и меньше 1. Для детализации самых высоких частот сигнала минимальных размер окна вейвлета не должен превышать периода самой высокочастотной гармоники. Если в сигнале присутствуют спектральные компоненты, соответствующие текущему значению а, то интеграл произведения вейвлета с сигналом в интервале, где эта спектральная компонента присутствует, дает относительно большое значение. В противном случае - произведение мало или равно нулю, т.к. среднее значение вейвлетной функции равно нулю. С увеличением масштаба (ширины окна) вейвлета преобразование выделяет все более низкие частоты. Рис. 1.2.2. На рис. 1.2.2 приведен пример модельного сигнала и спектра его непрерывного вейвлет-преобразования. Значения параметров 'а' и 'b' в (1.2.2) являются непрерывными, и множество базисных функций является избыточным. Сигналу, определенному на R, соответствует вейвлетный спектр R × R. Отсюда следует, что вейвлетный спектр НПВ имеет огромную избыточность. Обратное преобразование. Так как форма базисных функций (a,b,t) зафиксирована, то вся информация о сигнале в (1.2.1) переносится на значения функции С(a,b). Точность обратного интегрального вейвлет-преобразования зависит от выбора базисного вейвлета и способа построения базиса, т.е. от значений базисных параметров a, b. Строго теоретически вейвлет может считаться базисной функцией L2(R) только в случае его ортонормированности. Для практических целей непрерывного преобразования часто бывает вполне достаточна устойчивость и "приблизительность" ортогональности системы разложения функций. Под устойчивостью понимается достаточно точная реконструкция произвольных сигналов. Для ортонормированных вейвлетов обратное вейвлет-преобразование записывается с помощью того же базиса, что и прямое: s(t) = (1/C) (1/a2) С(a,b) (a,b,t) da db. (1.2.3) где C - нормализующий коэффициент: C = (|()|2 /) d < . (1.2.4) Условие конечности C ограничивает класс функций, которые можно использовать в качестве вейвлетов. В частности, при ω=0, для обеспечения сходимости интеграла (1.2.4) в нуле, значение () должно быть равно нулю. Это обеспечивает условие компактности фурье-образа вейвлета с локализацией вокруг некоторой частоты o – средней частоты вейвлетной функции. Следовательно, функция (t) должна иметь нулевое среднее значение по области его определения (интеграл функции по аргументу должен быть нулевым): (t) dt =0. Однако это означает, что не для всех сигналов возможна их точная реконструкция вейвлетом (t), т.к. при нулевом первом моменте вейвлета коэффициент передачи постоянной составляющей сигнала в преобразовании (1.2.3) равен нулю. Условия точной реконструкции сигналов будут рассмотрены при описании кратномасштабного анализа. Кроме того, даже при выполнении условия (1.2.4) далеко не все типы вейвлетов могут гарантировать реконструкцию сигналов, как таковую. Однако и такие вейвлеты могут быть полезны для анализа особенностей сигналов, как дополнительного метода к другим методам анализа и обработки данных. В общем случае, при отсутствии строгой ортогональности вейвлетной функции (1.2.1), для обратного преобразования применяется выражение: s(t) = (1/C)(1/a2) С(a,b) #(a,b,t) da db, (1.2.3') где индексом #(a,b,t) обозначен ортогональный "двойник" базиса (a,b,t), о котором будет изложено ниже. Рис. 1.2.3. Таким образом, непрерывное вейвлет-преобразование представляет собой разложение сигнала по всем возможным сдвигам и сжатиям/растяжениям некоторой локализованной финитной функции - вейвлета. При этом переменная 'a' определяет масштаб вейвлета и эквивалентна частоте в преобразованиях Фурье, а переменная 'b' – сдвиг вейвлета по сигналу от начальной точки в области его определения, шкала которого повторяет временную шкалу анализируемого сигнала. Вейвлетный анализ является частотно-пространственным анализом сигналов. В качестве примера рассмотрим вейвлет-преобразование чистого гармонического сигнала s(t), приведенного на рис. 1.2.3. На этом же рисунке ниже приведены вейвлеты a(t) симметричного типа разных масштабов. Скалярное произведение (1.2.1) "просмотра" сигнала вейвлетом определенного масштаба 'a' может быть записано в следующей форме: Ca(b)= s(t), a(t+b) =s(t)a(t+b) dt. (1.2.5) Но выражение (1.2.5) эквивалентно взаимной корреляционной функции Ra(b) сигналов s(t) и а(t). Если сигнал s(t) представляет собой гармонику, а второй сигнал симметричен, задан на компактном носителе и имеет нулевое среднее значение, то, как известно, форма взаимной корреляционной функции таких сигналов также является центрированным гармоническим сигналом. В частотной области скалярное произведение двух функций отображается произведением Фурье-образов этих функций, которые приведены на рисунке в правом столбце спектров. Масштабы спектров a() и Ra() для наглядности сопоставления нормированы к спектру s(t). Максимальная амплитуда гармоники Rа(b) будет наблюдаться при совпадении средней частоты локализации вейвлета а(t) определенного масштаба 'а' в частотной области с частотой сигнала s(t), что и можно видеть на рис. 1.2.3 для функции Ra(b) при масштабе вейвлета a=20. Результирующий вейвлетный спектр непрерывного вейвлет-преобразования гармоники приведен на левом нижнем графике и показывает точное положение на временной оси 'b' максимумов и минимумов гармонического сигнала. Дискретное вейвлет-преобразование. В принципе, при обработке данных на ПК может выполняться дискретизированная версия непрерывного вейвлет-преобразования с заданием дискретных значений параметров (a, b) вейвлетов с произвольным шагом a и b. В результате получается избыточное количество коэффициентов, намного превосходящее число отсчетов исходного сигнала, которое не требуется для реконструкции сигналов. Дискретное вейвлет-преобразование (ДВП) обеспечивает достаточно информации, как для анализа сигнала, так и для его синтеза, являясь вместе с тем экономным по числу операций и по требуемой памяти. ДВП оперирует с дискретными значениями параметров а и b, которые задаются, как правило, в виде степенных функций: a = ао-m, b = k·ао-m, ao > 1, m, k  I, где I – пространство целых чисел {-, }, m – параметр масштаба, k – параметр сдвига. Базис пространства L2(R) в дискретном представлении: mk(t) = |ао|m/2(аоmt-k), m,k  I, (t) Î L2(R). (1.2.6) Вейвлет-коэффициенты прямого преобразования: Cmk =s(t)mk(t) dt. (1.2.7) Значение 'a' может быть произвольным, но обычно принимается равным 2, при этом преобразование называется диадным вейвлет-преобразованием. Для диадного преобразования разработан быстрый алгоритм вычислений, аналогичный быстрому преобразованию Фурье, что предопределило его широкое использование при анализе массивов цифровых данных. Обратное дискретное преобразование для непрерывных сигналов при нормированном ортогональном вейвлетном базисе пространства: s(t) = Cmkmk(t). (1.2.8) Число использованных вейвлетов по масштабному коэффициенту m задает уровень декомпозиции сигнала, при этом за нулевой уровень (m = 0) обычно принимается уровень максимального временного разрешения сигнала, т.е. сам сигнал, а последующие уровни (m < 0) образуют ниспадающее вейвлет-дерево. В программном обеспечении вычислений для исключения использования отрицательной нумерации по m знак 'минус' обычно переносится непосредственно в (1.2.6), т.е. используется следующее представление базисных функций: mk(t) = |ао|-m/2(ао-mt-k), m,k  I, (t) Î L2(R). (1.2.6') Устойчивость дискретного базиса определяется следующим образом. Функция (t)Î L2(R) называется R-функцией, если базис на ее основе по (1.2.6) является базисом Рисса (Riesz). Для базиса Рисса существуют значения А и В, 0 < A ≤ B < , для которых выполняется соотношение A||Cmk||2 ≤ || Cmkmk(t)||2 ≤ B||Cmk||2, если энергия ряда Cmk конечна. При этом для любой R-функции существует базис #mk(t), который ортогонален базису mk(t). Его называют ортогональным "двойником" базиса mk(t), таким, что mk(t), #nl(t) = mn·kl. Если A = B = 1 и ао = 2, то семейство базисных функций {mk(t)} является ортонормированным базисом и возможно полное восстановление исходного сигнала, при этом mk(t) ≡ #mk(t) и для реконструкции сигналов используется формула (1.2.8). Если (t) не ортогональный вейвлет, но имеет "двойника", то на базе "двойника" вычисляется семейство #mk(t), которое и используется при обратном преобразовании вместо mk(t), при этом точное восстановление исходного сигнала не гарантировано, но оно будет близко к нему в среднеквадратическом смысле. Как и для непрерывного вейвлет-преобразования, обратное дискретное преобразование (1.2.8) не может выполнить восстановление нецентрированных сигналов в силу нулевого первого момента вейвлетных функций и, соответственно, центрирования значения вейвлет-коэффициентов Cmk при прямом вейвлет-преобразовании. Поэтому при обработке числовых массивов данных дискретные вейвлеты используются, как правило, в паре со связанными с ними дискретными скейлинг-функциями. Скейлин-функции имеют с вейвлетами общую область задания и определенное соотношение между значениями, но первый момент скейлин-функций по области определения равен 1. Если вейвлеты рассматривать, как аналоги полосовых фильтров сигнала, в основном, высокочастотных при выделении локальных особенностей в сигнале, то скейлин-функции вейвлетов представляет собой аналоги низкочастотных фильтров, которыми из сигнала выделяются в отдельный массив составляющие, не прошедшие вейвлетную фильтрацию. Так, например, порождающая скейлинг-функция вейвлета Хаара (1.1.15) задается следующим выражением: При обозначении скейлинг-функций индексом mk(t) аналитика скейлин-функций повторяет выражения (1.2.6-1.2.7) и образует дополнительный базис пространства L2(R). Сумма вейвлет-коэффициентов и скейлинг-коэффициентов разложения сигналов соответственно дает возможность выполнять точную реконструкцию сигналов, при этом вместо (1.2.8) используется следующее выражение обратного вейвлет-преобразования: s(t) =Сak k(t) + Сdmkmk(t), (1.2.9) где Cak – скейлин-коэффициенты, которые обычно называют коэффициентами аппроксимации сигнала, Cdmk – вейвлет-коэффициенты или коэффициенты детализации. Более подробно использование скейлинг-функций будет рассмотрено в теме вейвлетного кратномасштабного анализа. Частотно-временная локализация вейвлет-анализа. Реальные сигналы, как правило, конечны и принадлежат пространству L2(R). Частотный спектр сигналов обратно пропорционален их длительности. Соответственно, достаточно точный низкочастотный анализ сигнала должен производиться на больших интервалах его задания, а высокочастотный – на малых. Если частотный состав сигнала претерпевает существенные изменения на интервале его задания, то преобразование Фурье дает только усредненные данные частотного состава сигнала с постоянным частотным разрешением. Определенная частотно-временная локализация анализа создается применением оконного преобразования Фурье, что дает семейства частотных спектров, локализованных во времени, но в пределах постоянной ширины окна оконной функции, а, следовательно, также с постоянным значением и частотного, и временного разрешения. В отличие от оконного преобразования Фурье, вейвлет-преобразование, при аналогичных дискретных значениях сдвигов b, дает семейства спектров масштабных коэффициентов а сжатия-растяжения С(a,b) =s(t)|а|-1/2о[(t-b)/а]dt. Если считать, что каждый вейвлет имеет определенную "ширину" своего временного окна, которому соответствует определенная "средняя" частота спектрального образа вейвлета, обратная его масштабному коэффициенту а, то семейства масштабных коэффициентов вейвлет-преобразования можно считать аналогичными семействам частотных спектров оконного преобразования Фурье, но с одним принципиальным отличием. Масштабные коэффициенты изменяют "ширину" вейвлетов и, соответственно, "среднюю" частоту их фурье-образов, а, следовательно, каждой частоте соответствует своя длительность временного окна анализа, и наоборот. Так малые значения параметра а, характеризующие быстрые составляющие в сигналах, соответствуют высоким частотам, а большие значения – низким частотам. За счёт изменения масштаба вейвлеты способны выявлять различия на разных частотах, а за счёт сдвига (параметр b) проанализировать свойства сигнала в разных точках на всём исследуемом временном интервале. Многоразмерное временное окно вейвлет-преобразования адаптировано для оптимального выявления и низкочастотных, и высокочастотных характеристики сигналов. Для произвольной оконной функции z(t)Î L2(R) ее центр и радиус определяются формулами: to =t |z(t)|2 dt, z =  Если по этим функциям определить центры и радиусы вейвлетов и их фурье-образов, то временная локализация происходит с центрами окон b+ato шириной wint = 2atа частотная – с центрами ωо/а, и с шириной окна winω = 2/а. При этом значение отношения центральной частоты к ширине окна не зависит от местоположения центральной частоты. Частотно-временное окно wint·winω = 4tсужается при высокой центральной частоте, и расширяется при низкой. Схематическое изображение частотно-временных окон преобразования приведено на рис. 1.2.4. Таким образом, на высоких частотах лучше разрешение по времени, а на низких - по частоте. Для высокочастотной компоненты сигнала мы можем точнее указать ее временную позицию, а для низкочастотной - ее значение частоты. Изменение частотно-временного окна вейвлета определяет угол влияния значений функции в произвольных точках ti на значения коэффициентов С(а,b). И наоборот, угол влияния из точки С(ai,bi) на ось t определяет интервал значений функции, которые принимают участие в вычислении данного коэффициента С(ai,bi) – область достоверности. Схематически это показано на рис. 1.2.5. По углу влияния наглядно видно, что высокочастотная (мелкомасштабная) информация вычисляется на основе малых интервалов сигналов, а низкочастотная – на основе больших. Поскольку анализируемые сигналы всегда конечны, то при вычислении коэффициентов на границах задания сигнала область достоверности выходит за пределы сигнала, и для уменьшения погрешности вычислений сигнал дополняется заданием начальных и конечных условий. Образное представление преобразования. Представим себе длинный и узкий стеклянный ларь, произвольно заполненный шарами трех разных диаметров: 5, 10 и 15 см. Взглянем на ларь сбоку, и линию высоты насыпки будем считать значением сигнала в зависимости от расстояния от одного из торцов ларя (условно – нулевого). Возьмем первый "вейвлет" – идеальное дифференциальное сито с диаметром отверстий d=5 см, через которое проходят только пятисантиметровые шары (аналог значения ao). Передвигаясь вдоль ларя, "просеем" через это сито шары в ларе, не перемешивая их по расстоянию от нулевого торца ларя и размещая отсеиваемые шары в таком же ларе, сохраняя расстояние от начала ларя. Сменим масштаб "вейвлета" и повторим эту операцию ситом с диаметром отверстий 10, а затем 15 см. Если все три ларя расположить радом, мы получим двумерную "поверхность" насыпки отсеянных шаров, которая наглядно покажет распределение шаров в ларе и по размерам, и по их концентрации в различных участках ларя. Данная модель разложения является довольно грубой, но интуитивно понятно, что обратная сборка шаров в ларь с сохранением их местоположения с определенной точностью восстановит высоту насыпки. Замените шары короткими фрагментами электронных сигналов произвольной, но одной и той формы в пределах диаметра шаров, например такими, как (t) на рис.1.2.1, сложите все значения сигналов по текущим значениям t, и Вы получите сложный суммарный сигнал. Используя прямое вейвлет-преобразование с вейвлетами этих же составляющих, Вы можете разложить суммарный сигнал (и любой другой произвольный сигнал) на составляющие в масштабно-временной плоскости. Замените масштабную ось ширины вейвлетов на обратную ей частотную ось, и Вы представите результаты в частотно-временной плоскости. Заметим только, что точность, представительность и информативность результатов анализа во многом будут зависеть как от формы и особенностей анализируемого сигнала, так и от формы выбранных вами вейвлетов и параметров масштабирования и сдвига. Это определяется тем, что дифференциальное сито в примере с шарами – идеальная операция разделения, в то время как при вейвлет-преобразовании "идентификация" составляющих выполняется по скалярному произведению сигнала и функции вейвлета. Скалярное произведение в принципе не может давать однозначного ответа типа "да-нет", а только "наносит" на масштабно-временную плоскость определенные значения величины скалярного произведения. С одной стороны, выбор типа вейвлета вносит определенную субъективность исследователя в методику исследования сигналов, но, с другой стороны, дает исследователю новые возможности и свободу в поиске наиболее эффективных и оптимальных методов обработки сигналов и извлечения из них необходимой информации. Достоинства и недостатки вейвлетных преобразований. Вейвлетные преобразования обладают всеми достоинствами преобразований Фурье. Вейвлетные базисы могут быть хорошо локализованными как по частоте, так и по времени. При выделении в сигналах хорошо локализованных разномасштабных процессов можно рассматривать только те масштабные уровни разложения, которые представляют интерес. Вейвлетные базисы, в отличие от преобразования Фурье, имеют много разнообразных базовых функций, свойства которых ориентированы на решение различных задач. Базисные вейвлеты могут реализоваться функциями различной гладкости. Недостатком вейвлетных преобразований является их относительная сложность. Практическое использование вейвлет-преобразований связано, в основном, с дискретными вейвлетами как в силу повсеместного использования цифровых методов обработки данных, так и в силу ряда различий дискретного и непрерывного вейвлет-преобразований. Непрерывные вейвлеты дают несколько более наглядное представление результатов анализа в виде поверхностей вейвлет-коэффициентов по непрерывным переменным. На рис. 1.2.6 анализируемый сигнал состоит из двух модулированных гауссианов. Преобразование вейвлетом Морлета четко показывает их пространственную и частотную локализацию, в то время как спектр Фурье дает только частотную локализацию. Однако базисы на основе непрерывных вейвлетов, как правило, не являются строго ортонормированными, поскольку элементы базиса бесконечно дифференцируемы и экспоненциально спадают на бесконечности. У дискретных вейвлетов эти проблемы легко снимаются, что обеспечивает более точную реконструкцию сигналов. Выбор конкретного вида и типа вейвлетов во многом зависит от анализируемых сигналов и задач анализа, при этом немалую роль играет интуиция и опыт исследователя. Для получения оптимальных алгоритмов преобразования разработаны определенные критерии, но их еще нельзя считать окончательными, т.к. они являются внутренними по отношению к самим алгоритмам преобразования и, как правило, не учитывают внешних критериев, связанных с сигналами и целями их преобразований. Отсюда следует, что при практическом использовании вейвлетов необходимо уделять достаточное внимание проверке их работоспособности и эффективности для поставленных целей по сравнению с известными методами обработки и анализа. КОНТРОЛЬНЫЕ ВОПРОСЫ 1.Перечислите достоинства и недостатки вейвлетных преобразований 2.Опишите образное представление преобразования. 3.Что такое масштаб ВП? 4.Дайте характеристику вейвлетный спектр 5.Дайте характеристику частотно-временное оконное преобразование ЛЕКЦИЯ 7. СИСТЕМЫ, ОПИСЫВАЕМЫЕ РАЗНОСТНЫМИ УРАВНЕНИЯМИ. ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА. ПЕРЕДАТОЧНАЯ ФУНКЦИЯ. ЧАСТОТНАЯ ХАРАКТЕРИСТИКА. Этот класс моделей описывает детерминированные системы с непрерывным состоянием и дискретным временем перехода. Если интервалы между моментами перехода равны, то удобно пронумеровать рассматриваемые моменты времени, и тогда такую систему можно записать в следующем общем виде  – уравнение состояния, (6)    – уравнение выхода. (7) Соответственно для линейных систем: A(k) – переходная матрица состояний на  – ом шаге вычисления решения. Уравнение с дискретным временем может, например, быть применено в случаях, когда измерения и воздействия на непрерывную систему осуществляются в дискретные моменты времени. Если принять интервал дискретности равным  , то можно в соответствии с (5) записать для линейной непрерывной системы:   ;  . (8) Если принять  на интервале постоянным и равным  , то можно (8) представить в виде     (9)    (10) Если на вход линейного звена или системы звеньев подать гармонический сигнал, то на выходе после окончания переходного процесса также будет гармонический сигнал. Причем частота выходного сигнала, всегда будет равна частоте входного. (Вспомните результаты лабораторной работы, когда гармонический (синусоидальный) сигнал подавали на вход безинерционного звена (масштабирующий блок Gain)). Рассмотрим пример. Подадим на вход блока, моделирующего статическое колебательное инерционное звено синусоидальный сигнал  и посмотрим какой сигнал y(t) получится на выходе данного звена. Рис. 10.1 – Схема моделирования     Рис. 10.3 – Сигналы на входе x(t) и выходе y(t) колебательного звена   Как видно по рис. 10.3, после прохождения через колебательное звено и окончания переходного процесса мы наблюдаем гармонический процесс с частотой равной частоте входного сигнала. Амплитуда и сдвиг фаз выходного процесса изменились по сравнению с входным. Частотными характеристиками звена (системы) называются характеристики изменения параметров динамически установившихся гармонических колебаний на выходе системы в функции значения частоты (w) гармонического сигнала входного воздействия, когда w Î (-¥; ¥). Для графического отображения характеристик, как правило, используют положительную полуось частоты. Существую следующие виды частотных характеристик: - амплитудочастотная характеристика (АЧХ); - фазочастотная характеристика (ФЧХ); - амплитудофазочастотная характеристика (АФЧХ). АЧХ (А(w)) звена (системы) - отношение амплитуды выходного гармонического сигнала к амплитуде входного в функции его частоты w в установившемся режиме. Т.е. если x(t) = xmaxsinwt – входной сигнал; y(t) = ymaxsin(wt + φ) – выходной сигнал, то  , w Î (-¥, ¥). (10.1) ФЧХ (j(w)) звена (системы) называется разность фаз между выходным и входным сигналами в установившемся режиме. j(w) = jy(w) – jx; (10.2)  jx = const ≡ 0;  ;  w Î (-¥, ¥). Рис. 10.4 – Графики входного и выходного сигналов с построениями необходимыми для определения параметров, для расчета АЧХ и ФЧХ АФЧХ – это годограф радиус-вектора, длина которого – АЧХ, а угол поворота – ФЧХ. Частотная передаточная функция – отношение изображений Фурье выходного и входного сигналов.  ,  где  - оператор Фурье. Для линейных устойчивых систем при нулевых начальных условиях возможен взаимный однозначный переход между передаточными функциями систем в непрерывном времени. Для таких процессов преобразование Лапласа переходит в преобразование Фурье, если произвести подстановку s = jw. Частотная ПФ связана с частотными характеристиками следующим образом: частотная передаточная функция  представляет собой комплексное число, модуль которого равен отношению амплитуды выходной величины к амплитуде входной (АЧХ), а аргумент – сдвигу фаз выходной величины по отношению к входной(ФЧХ).  ; (10.4) Частотную передаточную функцию можно представить в алгебраической форме: W(jω) = U(ω) + jV(ω) = ReW(jω) + jImW(jω). Взаимосвязь между различными формами представления W(jω):  ;  . U(ω) и V(ω) можно получить непосредственно из W(jω). Для этого необходимо представить частотную ПФ в следующем виде:  . Затем умножить числитель и знаменатель на выражение комплексно сопряженное знаменателю. Затем выделить действительную и мнимую части. Пример: получение частотной передаточной функции из передаточной функции в форме изображений по Лапласу для статического апериодического инерционного звена 1-го порядка. W(s) =  Общая форма частотной ПФ получается формальной заменой оператора s на  : W(jw) =  . Чтоб перейти к алгебраической форме помножим числитель и знаменатель на выражение комплексно-сопряженное знаменателю Вывести частотные ПФ самостоятельно для всех звеньев. Итак, частотные характеристики можно построить различными способами: 1. Можно аналитически получить частотную передаточную функцию звена (системы), а затем промоделировать выражения (10.4). 2. Можно подавать на вход звена (системы) гармонические входные воздействия с различной частотой и вычислять АЧХ и ФЧХ по формулам (10.1), (10.2). Эти два метода являются достаточно трудоемкими, особенно когда речь идет о построении частотных характеристик сложных систем. 3. Можно воспользоваться средствами компьютерного моделирования. Специальный блок, разработанный в Simulink – «Анализатор частотных характеристик» производит построение АЧХ, ФЧХ и АФЧХ звена или системы. КОНТРОЛЬНЫЕ ВОПРОСЫ: 1. Дайте определение понятию - частотная передаточная функция  2. Напишите формулу - частотной передаточной функции 3. Перечислите и охарактеризуйте часточные характеристики звена. 1. ЛЕКЦИЯ 8. МЕТОДЫ СИНТЕЗА КИХ-ФИЛЬТРОВ. СИНТЕЗ ФИЛЬТРОВ ПО МЕТОДУ ОКНА. ВЕСОВЫЕ ФУНКЦИИ В МЕТОДЕ ОКНА Различают два общих класса сигналов: аналоговые и дискретные. Аналоговым сигналом называется сигнал, определенный для каждого момента времени. Дискретным сигналом - сигнал, определенный только в дискретные моменты времени. Как дискретный, так и аналоговый сигналы могут быть однозначно представлены некоторыми функциями частоты, которые называются их частотными спектрами. Фильтрацией называется процесс изменения частотного спектра сигнала в некотором желаемом направлении. Этот процесс может привести к усилению или ослаблению частотных составляющих в некотором диапазоне частот, к подавлению или выделению какой-нибудь конкретной составляющей и т.п. Рассмотрим аналоговый НЧ-фильтр. Импульсная характеристика: (1) Отфильтрованный сигнал: Переходя к дискретному времени получим: Ограничивая сумму до N и изменяя свертку получается КИХ-фильтр: Рассмотрим z-преобразование данного уравнения: , где . Подставим (1) Следовательно , проводя обратное z-преобразование получаем - это рекурсивный фильтр. Цифровые фильтры и их свойства. Цифровым фильтром называется цифровая система, которую можно использовать для фильтрации дискретных сигналов. Он может быть реализован программным методом или с помощью специальной аппаратуры, и в каждом из этих случаев цифровой фильтр можно применить для фильтрации сигналов в реальном времени или для фильтрации предварительно записанных сигналов. Цифровой фильтр можно представить структурной схемой, изображенной на рис.1. На этой схеме x(n) и y(n) - соответственно входное воздействие и реакция фильтра на это воздействие. Функционально они связаны соотношением : где вид оператора зависит от свойств конкретной системы. рис.1. Реакцию цифрового фильтра на произвольное воздействие можно представить с помощью импульсной характеристики фильтра. Допустим, что x(n) - входная , а y(n) - выходная последовательности фильтра, и пусть h(n) - отклик на единичный импульс, называемый импульсной характеристикой. Тогда Таким образом, x(n) и y(n) связаны соотношением типа свертки. Частотная характеристика фильтра определяется следующим выражением: (1.0) Поскольку частотная характеристика является периодической функцией частоты w, равенство (1.0) можно рассматривать как разложение в ряд Фурье, причем коэффициенты являются одновременно отсчетами импульсной характеристики. Согласно теории рядов Фурье, коэффициенты h(n) могут быть выражены через : Из этого соотношения видно, что h(n) по существу является суперпозицией синусоид с амплитудами . можно представить следующим образом: называют амплитудной характеристикой фильтра, а - фазовой характеристикой фильтра. Свойства цифровых фильтров. 1. Цифровой фильтр называется стационарным, если его параметры не изменяются во времени, т.е. предварительно невозбужденный фильтр, в котором x(n)=y(n)=0 при всех n<0 называют стационарным тогда и только тогда, когда для всех возможных воздействий. 2. Цифровой фильтр называют линейным тогда и только тогда, когда для всех a и b - произвольных постоянных и всех допустимых воздействий x1(n) и x2(n). 3. Цифровой фильтр называют физически реализуемым, если величина отклика при n=n0 зависит только от значений входной последовательности с номерами n£n0. Это означает, что импульсная характеристика h(n) равна нулю при n<0. 4. Цифровой фильтр называется устойчивым тогда и только тогда, когда реакция на ограниченное воздействие ограничена, т.е. если из при всех n следует при всех n. Необходимым и достаточным условием устойчивости фильтра является следующее требование к его импульсной характеристике: . Один наиболее важных классов линейных инвариантных относительно сдвига систем описывается разностным уравнением , (2.11) где x(n) – выборки входной последовательности; y(n) – выборки выходной последовательности, а a(k) и b(k) – коэффициенты, определяющие систему или, точнее, фильтр. В z-области система (2.11) может быть представлена своей передаточной функцией H(z): . (2.12) В случае, когда b(0)=1; b(k)=0, , выборки на выходе фильтра зависят лишь от входных выборок без какого-либо обратного влияния предыдущих выборок на текущие. Такие фильтры называются КИХ-фильтрами. Если по крайней мере один из коэффициентов b(k) при отличен от нуля, то такие фильтры называются фильтрами с бесконечной импульсной характеристикой (БИХ). КИХ-фильтры. Методы синтеза. Класс последовательностей конечной длины (КИХ-последовательности) обладает некоторыми свойствами, желательными с точки зрения построения фильтров. Например, никогда не возникает вопрос об устойчивости и физической реализуемости фильтров, поскольку КИХ-последовательности гарантируют устойчивость. Более того, КИХ-последовательности можно выбрать так, чтобы фильтры имели строго линейные фазовые характеристики. Поэтому, используя КИХ-последовательности, можно проектировать фильтры с произвольной амплитудной характеристикой. Существуют три основных метода синтеза КИХ-фильтров: метод взвешивания (метод «окна»); метод частотной выборки; метод оптимальных фильтров. Свойства КИХ-фильтров Имеется много причин, побуждающих к изучению способов проектирования КИХ-фильтров. Основными достоинствами этих фильтров являются: 1. Легко создавать КИХ-фильтры со строго линейной фазовой характеристикой (постоянной групповой задержкой). Во многих случаях, когда проектируется фильтр с произвольной амплитудной характеристикой, это упрощает задачу аппроксимации. 2. КИХ-фильтры, реализуемые нерекурсивно, т.е. с помощью свертки, всегда устойчивы. 3. При нерекурсивной реализации КИХ-фильтров шумы округления, возникающие за счет выполнения арифметических операций с конечной точностью, легко минимизировать. 4. КИХ-фильтры можно эффективно реализовывать при помощи методов быстрой свертки, основанных на применении алгоритма БПФ. Недостатки КИХ-фильтров следующие: 1. Для аппроксимации фильтров, частотные характеристики которых имеют острые срезы, требуется импульсная характеристика с большим числом отсчетов N. Следовательно, возрастает объем вычислительных операций. 2. Задержка в КИХ-фильтрах с линейной фазовой характеристикой не всегда равна целому числу интервалов дискретизации. Характеристики КИХ-фильтров с линейной фазовой характеристикой Пусть {h(n)} – физически реализуемая последовательность конечной длины, заданная на интервале . Ее z-преобразование равно . (3.1) Преобразование Фурье от {h(n)} (3.2) является периодическим по частоте с периодом , т.е. . (3.3) Рассматривая только действительные последовательности {h(n)}, получим дополнительные ограничения на функцию , представив ее через амплитуду и фазу: . (3.4) Потребуем при расчете КИХ-фильтров строго линейной фазовой характеристики, и рассмотрим, при каких условиях импульсная характеристика фильтра h(n) будет это обеспечивать. Требование линейности фазовой характеристики имеет вид , (3.5) где - постоянная фазовая задержка, выраженная через число интервалов дискретизации. Используя (3.4) и (3.5) соотношение (3.2) переписывается в виде: . (3.6) Приравнивая действительные и мнимые части, и деля друг на друга правые и левые части полученных равенств, можно получить уравнение, решением которого будут следующие значения: , (3.7) . (3.8) Смысл их заключается в следующем. Условие (3.7) означает, что для каждого N существует только одна фазовая задержка , при которой может достигаться строгая линейность фазовой характеристики фильтра. Из условия (3.8) следует, что при заданном , удовлетворяющем условию (3.7), импульсная характеристика должна обладать симметрией. Рассмотрим использование условий (3.7) и (3.8) для случаев четного и нечетного N. Если N - нечетно, то задержка в фильтре равна целому числу интервалов дискретизации. Типичная импульсная характеристика фильтра с линейной фазой для случая N=11 приведена на рис. 3.1. Типичная импульсная характеристика фильтра с линейной фазой при четном N показана на рис. 3.2. Здесь N=10. Рис.3.1. N – нечетно. Рис.3.2. N – четно. Расчет КИХ-фильтров методом взвешивания Как было сказано ранее, частотную характеристику любого цифрового фильтра можно представить рядом Фурье: , (3.9) где . (3.10) Видно, что коэффициенты Фурье совпадают с коэффициентами импульсной характеристики цифрового фильтра. Использование этих соотношений для проектирования КИХ-фильтра связано с двумя трудностями. Во-первых, импульсная характеристика фильтра имеет бесконечную длину, поскольку суммирование в (3.9) производится в бесконечных пределах. Во-вторых, фильтр физически нереализуем, так как импульсная характеристика начинается в , т.е. никакая конечная задержка не сделает фильтр физически реализуемым. Один из возможных методов получения КИХ-фильтра, аппроксимирующего заданную функцию , заключается в усечении бесконечного ряда Фурье (3.9) за . Однако простое усечение ряда приводит к явлению Гиббса, которое проявляется в виде выбросов и пульсаций до и после точек разрыва в аппроксимируемой частотной характеристики. Причем, максимальная амплитуда пульсаций частотной характеристики не уменьшается с увеличением длины импульсной характеристики, т.е. учет все большего числа членов ряда Фурье не приводит к уменьшению максимальной амплитуды пульсаций. Вместо этого уменьшается ширина выброса. Поэтому простое усечение ряда Фурье (3.9) не приводит к приемлемой аппроксимации идеального фильтра нижних частот (к чему необходимо стремиться). Этот метод непригоден для проектирования КИХ-фильтров. Лучшие результаты дает метод, основанный на использовании весовой последовательности конечной длины w(n), называемой окном, для модификации коэффициентов Фурье в формуле (3.9) с тем, чтобы управлять сходимостью ряда Фурье. Для большинства приемлемых окон преобразование Фурье последовательности w(n) имеет главный лепесток, содержащий почти всю энергию окна, и боковые лепестки, которые обычно быстро затухают. Чтобы получить КИХ-аппроксимацию функции , формируется последовательность , равная нулю за пределами интервала . Поскольку результирующая характеристика фильтра равна свертке идеальной частотной характеристики и частотной характеристика окна, то ширина переходных полос зависит от ширины главного лепестка функции . Кроме того, на всех частотах возникают ошибки аппроксимации, имеющие вид пульсаций частотной характеристики, которые обусловлены боковыми лепестками функции . Из приведенного выше следует, что оптимальное окно должно иметь во-первых, минимальную ширину главного лепестка частотной характеристики, во-вторых, минимальную площадь под боковыми лепестками. К сожалению, эти два требования несовместимы и необходим компромиссный вариант. Далее рассмотрим некоторые виды окон. Прямоугольное окно N-точечное прямоугольное окно, соответствующее простому усечению ряда Фурье, описывается весовой функцией (3.11) Частотная характеристика прямоугольного окна описывается соотношением . (3.12) Ее график представлен на рис. 3.3. Рис. 3.3. Частотная характеристика прямоугольного окна (N=25). «Обобщенное» окно Хэмминга Обобщенное окно Хэмминга имеет вид (3.13) причем . Случай соответствует окну Ханна (hanning), случай - окну Хэмминга. Частотную характеристику этого окна можно выразить через частотную характеристику прямоугольного окна (рис. 3.4): . (3.14) Рис.3.4. Частотная характеристика окна Хэмминга . Окно Кайзера Задача расчета хороших окон фактически сводится к математической задаче отыскания ограниченных во времени функций, преобразования Фурье которых наилучшим образом аппроксимируют функции, ограниченные по частоте, т.е. имеют минимальную энергию за пределами заданного интервала частот. Одним из решений такой задачи является окно Кайзера: (3.15) где - константа, определяющая компромисс между максимальным уровнем боковых лепестков и шириной главного лепестка частотной характеристики окна, а - функция Бесселя нулевого порядка. Окно Кайзера является по существу оптимальным окном в том смысле, что оно представляет последовательность конечной длины, которая имеет минимум энергии спектра за пределами некоторой заданной частоты. Прочие окна Существует еще много различных окон. Вот некоторые из них. Окно Блэкмана: . (3.16) Окно Фейера (треугольное окно): (3.17) Окно Ланцоша: , (3.18) где L – положительное целое число. Последние два окна основываются на методах суммирования рядов Фурье для ускорения их сходимости. Так же существуют окна, полученные согласно некоторым критериям оптимальности (окна Долфи-Чебышева, окно Каппелини, и т.д.). Особенности использования метода взвешивания Метод взвешивания весьма удобен для проектирования КИХ-фильтров, однако он обладает некоторыми особенностями, которые могут препятствовать применению окон. Прежде всего необходимо иметь выражения для коэффициентов ряда Фурье: . (3.19) Когда характеристика имеет сложный вид или не может быть просто преобразована в математическое выражение, формула (3.19) может оказаться громоздкой и неудобной для интегрирования. Еще одна особенность метода взвешивания заключается в отсутствии достаточной гибкости при проектировании фильтров. Например, при расчете ФНЧ обычно трудно определить граничную частоту полосы пропускания, поскольку окно «размывает» разрыв идеальной характеристики. Вышеупомянутые ограничения метода взвешивания не препятствуют его широкому применению. Чтобы определить невзвешенные коэффициенты Фурье в том случае, когда аналитическое выражение для h(n) громоздко или неудобно для интегрирования, интеграл можно аппроксимировать суммой по следующей формуле . (3.20) Ясно, что значения (3.20) можно эффективно вычислять с помощью M-точечного ОДПФ последовательности . Поскольку формула (3.20) является дискретным аналогом формулы (3.19), легко показать, что ростом M различие между h(n) и уменьшается. Расчет КИХ-фильтров методом частотной выборки КИХ-фильтр может быть однозначно задан коэффициентами импульсной характеристики {h(n)}, так и коэффициентами ДПФ частотной характеристики {H(k)}. Обе последовательности связаны соотношениями: ДПФ, (3.21) ОДПФ. (3.22) Из формулы (3.22) сразу вытекает прямой способ получения импульсной характеристики фильтра. Такой фильтр имеет частотную характеристику, значения которой равны H(k) в N равноотстоящих на оси частот точках. К сожалению, эта прямая процедура не представляет практического интереса, так как невозможно предсказать поведение частотной характеристики между частотными выборками. Для того, чтобы получить фильтры с линейной фазой, частотные выборки должны быть симметричными по амплитуде и иметь линейную антисимметричную фазу в интервале . Учитывая то, что будет использоваться ОДПФ, удобнее выразить условия симметрии на интервале . Ели частотные выборки записаны в виде , то условия симметрии при нечетном N можно записать в виде (3.23) Для улучшения качества аппроксимации часть частотных отсчетов имеет смысл сделать независимыми переменными. Значения этих независимых переменных обычно рассчитывают методами оптимизации, таким образом, чтобы минимизировать некоторую простую функцию ошибки аппроксимации. Можно сформулировать основную идею метода частотной выборки. Искомую частотную характеристику можно аппроксимировать ее отсчетами, взятыми в N равноотстоящих точках, а затем путем интерполяции получить результирующую частотную характеристику, которая будет проходить через исходные отсчеты. Ошибка интерполяции для фильтров с достаточно гладкими частотными характеристиками обычно имеет небольшую величину. В случае селективных фильтров частотные отсчеты в переходных полосах остаются незаданными переменными, значения которых подбираются с помощью алгоритма оптимизации. Для выполнения необходимой минимизации можно использовать простые методы линейного программирования. Проектирование оптимальных фильтров Смысл слова оптимальный заключается в следующем: оптимальными являются те фильтры, для которых максимальная ошибка в полосе пропускания и (или) в полосе задерживания минимальна по сравнению с любыми другими фильтрами, которые можно получить, изменяя значения n выборок в переходной полосе. Рассмотрим подход к расчету КИХ-фильтров при котором минимизируется максимальная ошибка аппроксимации. Пусть - заданная (желаемая) частотная характеристика, - аппроксимирующая функция, - положительная весовая функция, позволяющая определять ошибки для различных интервалов. Тогда взвешенная ошибка аппроксимации по определению равна . (3.24) Задачу чебышевской аппроксимации можно сформулировать как задачу поиска (а точнее коэффициентов, через которые можно выразить ), которая минимизирует максимум модуля ошибки в тех частотных полосах, где выполняется аппроксимация: , (3.25) где A – совокупность всех интересующих частотных полос. БИХ-фильтры. Методы синтеза. Фильтром с бесконечной импульсной характеристикой (БИХ - фильтром) называют фильтр, длина импульсной характеристики которого не ограничена справа или слева. Будем рассматривать БИХ-фильтры при условии, что они являются физически реализуемыми и устойчивыми. Для импульсных характеристик таких фильтров h(n) справедливы следующие ограничения: Наиболее общая форма записи z-преобразования импульсной характеристики БИХ-фильтров имеет вид: (3) Будем предполагать, что M£N. Системы, удовлетворяющие этому условию, называют системами N-го порядка. Решение задачи расчета фильтров сводится к нахождению значений его коэффициентов ai и bi, обеспечивающих аппроксимацию заданных характеристик фильтра. Таким образом, задача расчета фильтра в значительной мере сводится к задаче аппроксимации и может быть решена чисто математическими методами. Наиболее распространенным методом расчета цифровых БИХ-фильтров является метод дискретизации аналогового фильтра, удовлетворяющего заданным требованиям. При расчете цифровых фильтров верхних частот, полосовых и режекторных используются два подхода, представленные на рис. 2. рис.2. В первом случае нормализованный аналоговый фильтр предварительно преобразуется в другой аналоговый фильтр, из которого путем дискретизации рассчитывается фильтр с заданными характеристиками. Во втором случае нормализованный фильтр нижних частот дискретизуется сразу же, а затем путем преобразования его полосы частот формируется цифровой фильтр с заданными характеристиками. Аналоговые фильтры-прототипы. Приведем расчетные формулы для нескольких стандартных типов аналоговых фильтров. Пусть нужно рассчитать аналоговый фильтр нижних частот с частотой среза W=1 рад/с. В качестве аппроксимируемой функции будет использоваться квадрат амплитудной характеристики (исключением является фильтр Бесселя). Будем считать, что передаточная функция аналогового фильтра является рациональной функцией переменной S следующего вида: Фильтры Баттерворта. Фильтры Баттерворта нижних частот характеризуются тем, что имеют максимально гладкую амплитудную характеристику в начале координат в S-плоскости. Квадрат амплитудной характеристики нормированного (т.е. имеющего частоту среза 1 рад/с) фильтра Баттерворта равен: , где n - порядок фильтра. Аналитически продолжая функцию на всю S-плоскость, получим Все полюсы этой функции находятся на единичной окружности на одинаковом расстоянии друг от друга в S-плоскости. Выразим передаточную функцию H(S) через полюсы, располагающиеся в левой полуплоскости S: , где , k=1,2,...,n k0 - константа нормирования. Можно сформулировать несколько свойств фильтров Баттерворта нижних частот. · Фильтры Баттерворта имеют только полюсы (все нули передаточных функций этих фильтров расположены на бесконечности). · На частоте W=1 рад/с коэффициент передачи фильтра равен (т.е. на частоте среза их амплитудная характеристика спадает на 3 дБ). · Порядок фильтра n полностью определяет весь фильтр. Фильтры Чебышева. Отличительной чертой фильтров Чебышева является наименьшая величина максимальной ошибки аппроксимации в заданной полосе частот. В действительности ошибка аппроксимации представляется в заданной полосе частот равновеликими пульсациями, т.е. она флуктуирует между максимумами и минимумами равной величины. В зависимости от того, где минимизируется ошибка аппроксимации - в полосе пропускания или в полосе непропускания - различают фильтры Чебышева типа 1 и 2 . Фильтры Чебышева типа 1 имеют только полюсы и обеспечивают равновеликие пульсации амплитудной характеристики в полосе пропускания и монотонное изменение ослабления в полосе непропускания. Квадрат амплитудной характеристики фильтра Чебышева типа 1 n-го порядка описывается выражением: где - полином Чебышева n-го порядка, по определению равный e - параметр, характеризующий пульсации в полосе пропускания. Свойство оптимальности фильтров Чебышева типа 1 порядка n заключается в том, что не существует какого-либо другого фильтра n-го порядка, содержащего только полюсы, который имел бы такие же или лучшие характеристики и в полосе пропускания, и в полосе непропускания. Фильтры Чебышева типа 2 (иногда их называют обратными фильтрами Чебышева) обеспечивают монотонное изменение ослабления в полосе пропускания и равновеликие пульсации в полосе непропускания. Нули фильтров этого типа располагаются на мнимой оси в S-плоскости, а полюсы - в левой полуплоскости. Квадрат амплитудной характеристики фильтров Чебышева типа 2 порядка n можно представить следующим образом: где Wr - наинизшая частота, на которой в полосе непропускания достигается заданный уровень ослабления. Эллиптические фильтры. Эллиптические фильтры характеризуются тем, что их амплитудная характеристика имеет равновеликие пульсации и в полосе пропускания, и в полосе непропускания. Эллиптические фильтры являются оптимальными с точки зрения минимальной ширины переходной полосы. Квадрат амплитудной характеристики эллиптического фильтра нижних частот записывается в виде где - рациональная функция Чебышева, L - параметр, характеризующий пульсации функции . Фильтры Бесселя. Фильтры Бесселя характеризуются максимально гладкой характеристикой групповой задержки в начале координат в S-плоскости. Переходная характеристика фильтров Бесселя имеет весьма малый выброс. Однако при дискретизации непрерывных фильтров Бесселя методами, рассмотренными ниже, характерное для этих фильтров свойство максимальной гладкости характеристики групповой задержки не сохраняется. Передаточная функция фильтров Бесселя записывается в виде где B0(S) - функция Бесселя n-го порядка, d0 - константа нормирования: Функции Бесселя удовлетворяют следующему рекуррентному соотношению: с начальными условиями B0(S)=1 и B0(S)=S+1. Эти функции можно представить в виде: где k=1,2,...,n Фильтры Бесселя имеют только полюсы, которые расположены на окружности с центром на действительной положительной полуоси S-плоскости. Методы дискретизации аналогового фильтра. Предположим, что передаточная функция аналогового фильтра (представляющая собой преобразование Лапласа от импульсной характеристики) равна: (4) причем коэффициенты ai и bi (или ci и di) известны. Наиболее распространенными методами дискретизации аналогового фильтра с передаточной функцией (4) являются следующие: · метод отображения дифференциалов; · метод инвариантного преобразования импульсной характеристики; · метод билинейного преобразования; · метод согласованного z-преобразования. Рассмотрим некоторые из них. Метод инвариантного преобразования импульсной характеристики. Отличительной особенностью этого метода является то, что в качестве импульсной характеристики рассчитываемого цифрового фильтра используется дискретная импульсная характеристика соответствующего аналогового фильтра. В результате частотная характеристика цифрового фильтра образуется путем наложения частотной характеристики дискретизуемого аналогового фильтра. Разложим выражение (4) на простые дроби: где причем каждый коэффициент di определяет положение i-го полюса. При записи разложения (4) предполагалось, что порядок числителя M меньше порядка знаменателя N и все полюсы H(S) простые. Импульсная характеристика h(t) аналогового фильтра с передаточной функцией вида (4) описывается соотношением: Дискретизуя ее, получим импульсную характеристику цифрового фильтра: где T - период дискретизации. Найдем ее z-преобразование: Изменив порядок суммирования и просуммировав по n, получим: , (5) Сравним формулы (4) и (5). Видно, что для простых полюсов переход от H(S) к H(z) осуществляется с помощью отображения, при котором используется замена . (6) Если полюсы dj комплексные, то остатки cj в (4) также будут комплексными. Функция h(t) действительная, поэтому должны существовать также комплексно сопряженные полюс и остаток . Просуммируем эти комплексно сопряженные члены в (4): (7) Положим di=si+Wi и ci=gi+jhi, тогда (8) Использование отображающей замены применительно к каждому слагаемому в формуле (6) дает (9) Из формул (8) и (9) получаем (индекс i здесь опущен, а числители поделены на 2g) Частотная характеристика цифрового фильтра, рассчитываемого методом инвариантного преобразования импульсной характеристики, образуется путем наложения частотной характеристики дискретизуемого аналогового фильтра. Таким образом, можно записать где - угловая частота дискретизации цифрового фильтра. На рис.3 показано соответствующее инвариантному преобразованию импульсной характеристики отображение из S-плоскости в z-плоскость. Каждая горизонтальная полоса шириной из S-плоскости отображается на z-плоскость. Поэтому все смежные полосы из S-плоскости будут при отображении накладываться друг на друга в z-плоскости. рис.3. Значит, для того, чтобы частотные характеристики исходного аналогового фильтра и рассчитываемого методом инвариантного преобразования импульсной характеристики цифрового фильтра соответствовали друг другу, необходимо, чтобы полоса пропускания аналогового фильтра находилась в пределах диапазона . Для выполнения этого условия необходимо до начала преобразования вводить дополнительный фильтр нижних частот, гарантирующий соответствующее ограничение полосы пропускания аналогового фильтра. Метод инвариантного преобразования импульсной характеристики дает хорошие результаты для баттервортовских, бесселевых или чебышевских фильтров нижних частот и полосовых фильтров. Достоинство этого метода заключается в сохранении в цифровом фильтре таких же фазовой характеристики и характеристики затухания, как и у исходного аналогового фильтра. Метод билинейного преобразования. Билинейное z-преобразование использует следующую замену: (10) На рис. 4 показано, каким образом S-плоскость отображается в z-плоскость. рис.4. Вся ось jW из S-плоскости отображается в единичную окружность на z-плоскости, левая полуплоскость S отображается в единичный круг, а правая полуплоскость S - в область, расположенную вне единичного круга на z-плоскости. Из формулы (10) найдем выражение для z: (11) При S=jW получим: Отсюда видно, что . При =0 имеем z=1 и при W=¥ z=-1, в промежутке z монотонно меняется от 0 до . Подставив в формулу (11) S=s+jW, получим: При s<0 (для левой полуплоскости) , т.е. точки расположены внутри единичной окружности. При билинейном преобразовании передаточная функция цифрового фильтра H(z) рассчитывается с помощью алгебраической подстановки (10), т.е. Из этого соотношения видно, что порядки знаменателей функций H(z) и H(S) совпадают, но порядки числителей могут отличаться. Например, передаточная функция имеет числитель нулевого порядка, а знаменатель - первого порядка. В то же время получаемая методом билинейного преобразования функция равна где и числитель и знаменатель 1-го порядка. Причиной этого является то, что H(S) имеет нуль на бесконечности (S=¥), который при билинейном преобразовании отображается в точку z=-1. Так как в единичную окружность на z-плоскости отображается вся ось jW из S-плоскости, то эффекты, связанные с наложениями в частотной характеристике цифрового фильтра, характерные для метода инвариантного преобразования импульсной характеристики, в данном случае будут отсутствовать. Однако соотношение между частотами аналогового фильтра и цифрового фильтра w оказывается существенно нелинейным. Рассмотрим характер этой нелинейности. Пусть в (10) и S=jW, тогда или откуда При небольших w отображение почти линейно, однако для основной части частотной шкалы оно существенно нелинейно и сильно ограничивает область применения билинейного преобразования. Существует, правда, довольно большой класс фильтров, для которых частотная деформация может быть скомпенсирована. К ним относятся фильтры нижних частот, высоких частот, полосовые и режекторные. Эффекты нелинейности соотношения между частотными шкалами аналогового и цифрового фильтров удается учесть лишь в том случае, когда частотная характеристика аналогового фильтра имеет вид ступенчатой функции. Кроме того, при билинейном преобразовании ни импульсная, ни фазовая характеристики аналогового фильтров не будут совпадать. Преобразования полосы частот для аналоговых фильтров. Существует много различных методов преобразования фильтров нижних частот с частотой среза, равной 1 рад/с, в другой фильтр нижних частот (имеющий другую частоту среза), а также в фильтр верхних частот, полосовой или режекторный. Перечислим наиболее простые преобразования. 1. Фильтр нижних частот Фильтр нижних частот 2. Фильтр нижних частот Фильтр верхних частот 3. Фильтр нижних частот Полосовой фильтр 4. Фильтр нижних частот Режекторный фильтр Здесь Wr - нижняя частота среза, Wu - верхняя частота среза. Преобразование полосы для цифровых фильтров. По аналогии с фильтрами непрерывного времени существует несколько простых преобразований цифрового фильтра нижних частот (с частотой среза wc) в другой фильтр нижних частот (с другой частотой среза wu), а также в цифровой фильтр верхних частот, полосовой или режекторный. Ниже приведены формулы для этих преобразований. 1. Фильтр нижних частот Фильтр нижних частот где wu - заданная частота среза фильтра нижних частот. 2. Фильтр нижних частот Фильтр верхних частот где wu - заданная частота среза фильтра верхних частот 3. Фильтр нижних частот Полосовой фильтр Здесь w0 - центральная частота полосового фильтра 4. Фильтр нижних частот Режекторный фильтр Здесь wo - центральная частота режекторного фильтра Методы реализации цифровых фильтров. Цифровые фильтры с заданной передаточной функцией можно построить различными способами. В любом реальном фильтре шумы и погрешности, появляющиеся при квантовании, существенно зависят от структуры фильтра. Прежде всего все фильтры можно разделить на два больших класса: · рекурсивные; · нерекурсивные. Для рекурсивных фильтров соотношение между входной последовательностью {x(n)} и откликом фильтра {y(n)} может быть записано следующим образом: т.е. текущий отсчет отклика y(n) определяется не только текущим и предшествующим значениями входной последовательности, но и предшествующими отсчетами отклика. В нерекурсивных фильтрах связь между входной последовательностью и откликом имеет вид: т.е. текущий отсчет отклика зависит от текущего и предшествующих значений входной последовательности. Реализация может осуществляться на основе следующих форм построения схем фильтра: · прямой, · канонической прямой, · каскадной, · параллельной. Прямая форма. Рассмотрим передаточную функцию N-го порядка вида: (12) причем b0=1. Приведя это равенство к общему знаменателю, получим или Если рассматривать члены вида z-kY(z) как обратные z-преобразования последовательностей y(n-k), то, взяв обратные z-преобразования обеих частей последнего равенства, можно получить искомое разностное уравнение: Поскольку b0=1, разностное уравнение можно решить относительно y(n) : Простая структура реализации данного разностного уравнения показана на рис.5. Она носит название прямой формы. В ней для образования цепей, соответствующих числителю и знаменателю формулы (12), используются раздельные элементы задержки. Характерными чертами этой структуры являются ее простота и непосредственная связь с z-преобразованием. Однако эта структура очень чувствительна к квантованию коэффициентов. рис.5. Прямая каноническая форма. Запишем формулу (12) в следующем виде: Цифровой фильтр, соответствующий этой формуле, состоит из двух последовательно соединенных фильтров с коэффициентами передачи соответственно H1(z) и H2(z). Первый фильтр имеет только полюсы, а второй - только нули. Если записать то получится пара разностных уравнений (в предположении, что b0=1): которые можно реализовать, как показано на рис.6. Поскольку в цепях, соответствующих H1(z) и H2(z), сигнал w(n) задерживается одинаково, то для построения фильтра достаточно использовать один набор элементов задержки (рис.7). Такую структуру называют канонической формой, т.к. в ней используется минимальное количество сумматоров, умножителей и элементов задержки. рис.6. рис.7. Каскадная форма. Записав формулу (12) в виде: (13) получим третью структуру построения цифрового фильтра. Множители Hi(z) соответствуют либо блокам второго порядка, т.е. либо блокам первого порядка, т.е. а K равно целой части числа . Схему, реализующую формулу (13), называют каскадной (или последовательной) формой (рис.8). Каждый из блоков, образующих последовательную форму, можно реализовать в прямой или канонической форме. рис.8. Параллельная форма. Другим способом описания передаточной функции может быть ее представление разложением на простые дроби: (14) где слагаемые Hi(z) соответствуют или блокам второго порядка: или блокам первого порядка: причем К равно целой части от и . На рис.9 приведена структурная схема, реализующая соотношение (14). Ее называют параллельной формой. Блоки 1-го и 2-го порядка строятся по схеме одной из прямых форм. рис.9. КОНТРОЛЬНЫЕ ВОПРОСЫ: 1.Перечислите свойства КИХ-фильтров 2.Что такое фильтрация? Какая разница фильтрации разных сигналов 3.Перечислите и охарактеризуйте Методы синтеза 4.Что такое Фильтры Баттерворта? 5.Охарактеризуйте эллиптические фильтры 1. ЛЕКЦИЯ 9. ФИЛЬТРЫ НА ОСНОВЕ ЧАСТОТНОЙ ВЫБОРКИ. СИНТЕЗ ОПТИМАЛЬНОГО ФИЛЬТРА. ОСНОВНЫЕ МЕТОДЫ СИНТЕЗА КИХ-ФИЛЬТРОВ. СРЕДСТВА ДЛЯ СИНТЕЗА ЦИФРОВЫХ ФИЛЬТРОВ. Класс последовательностей конечной длины (КИХ-последовательности) обладает некоторыми свойствами, желательными с точки зрения построения фильтров. Например, никогда не возникает вопрос об устойчивости и физической реализуемости фильтров, поскольку КИХ-последовательности гарантируют устойчивость. Более того, КИХ-последовательности можно выбрать так, чтобы фильтры имели строго линейные фазовые характеристики. Поэтому, используя КИХ-последовательности, можно проектировать фильтры с произвольной амплитудной характеристикой. Существуют три основных метода синтеза КИХ-фильтров: • метод взвешивания (метод «окна»); • метод частотной выборки; • метод оптимальных фильтров. Свойства КИХ-фильтров Имеется много причин, побуждающих к изучению способов проектирования КИХ-фильтров. Основными достоинствами этих фильтров являются: 1. Легко создавать КИХ-фильтры со строго линейной фазовой характеристикой (постоянной групповой задержкой). Во многих случаях, когда проектируется фильтр с произвольной амплитудной характеристикой, это упрощает задачу аппроксимации. 2. КИХ-фильтры, реализуемые нерекурсивно, т.е. с помощью свертки, всегда устойчивы. 3. При нерекурсивной реализации КИХ-фильтров шумы округления, возникающие за счет выполнения арифметических операций с конечной точностью, легко минимизировать. 4. КИХ-фильтры можно эффективно реализовывать при помощи методов быстрой свертки, основанных на применении алгоритма БПФ. Недостатки КИХ-фильтров следующие: 1. Для аппроксимации фильтров, частотные характеристики которых имеют острые срезы, требуется импульсная характеристика с большим числом отсчетов N. Следовательно, возрастает объем вычислительных операций. 2. Задержка в КИХ-фильтрах с линейной фазовой характеристикой не всегда равна целому числу интервалов дискретизации. Характеристики КИХ-фильтров с линейной фазовой характеристикой Пусть {h(n)} – физически реализуемая последовательность конечной длины, заданная на интервале. Ееz-преобразование равно . (3.1) Преобразование Фурье от {h(n)} (3.2) является периодическим по частоте с периодом , т.е. . (3.3) Рассматривая только действительные последовательности {h(n)}, получим дополнительные ограничения на функцию, представив ее через амплитуду и фазу: . (3.4) Потребуем при расчете КИХ-фильтров строго линейной фазовой характеристики, и рассмотрим, при каких условиях импульсная характеристика фильтра h(n)будет это обеспечивать. Требование линейности фазовой характеристикиимеет вид , (3.5) где - постоянная фазовая задержка, выраженная через число интервалов дискретизации. Используя (3.4) и (3.5) соотношение (3.2) переписывается в виде: . (3.6) Приравнивая действительные и мнимые части, и деля друг на друга правые и левые части полученных равенств, можно получить уравнение, решением которого будут следующие значения: , (3.7) . (3.8) Смысл их заключается в следующем. Условие (3.7) означает, что для каждого Nсуществует только одна фазовая задержка, при которой может достигаться строгая линейность фазовой характеристики фильтра. Из условия (3.8) следует, что при заданном, удовлетворяющем условию (3.7), импульсная характеристика должна обладать симметрией. Рассмотрим использование условий (3.7) и (3.8) для случаев четного и нечетногоN. ЕслиN- нечетно, то задержка в фильтре равна целому числу интервалов дискретизации. Типичная импульсная характеристика фильтра с линейной фазой для случаяN=11 приведена на рис. 3.1. Типичная импульсная характеристика фильтра с линейной фазой при четномNпоказана на рис. 3.2. ЗдесьN=10. Рис.3.1. N – нечетно. Рис.3.2. N – четно. Расчет КИХ-фильтров методом взвешивания Как было сказано ранее, частотную характеристику любого цифрового фильтра можно представить рядом Фурье: , (3.9) где . (3.10) Видно, что коэффициенты Фурье совпадают с коэффициентами импульсной характеристики цифрового фильтра. Использование этих соотношений для проектирования КИХ-фильтра связано с двумя трудностями. Во-первых, импульсная характеристика фильтра имеет бесконечную длину, поскольку суммирование в (3.9) производится в бесконечных пределах. Во-вторых, фильтр физически нереализуем, так как импульсная характеристика начинается в , т.е. никакая конечная задержка не сделает фильтр физически реализуемым. Один из возможных методов получения КИХ-фильтра, аппроксимирующего заданную функцию , заключается в усечении бесконечного ряда Фурье (3.9) за. Однако простое усечение ряда приводит к явлению Гиббса, которое проявляется в виде выбросов и пульсаций до и после точек разрыва в аппроксимируемой частотной характеристики. Причем, максимальная амплитуда пульсаций частотной характеристики не уменьшается с увеличением длины импульсной характеристики, т.е. учет все большего числа членов ряда Фурье не приводит к уменьшению максимальной амплитуды пульсаций. Вместо этого уменьшается ширина выброса. Поэтому простое усечение ряда Фурье (3.9) не приводит к приемлемой аппроксимации идеального фильтра нижних частот (к чему необходимо стремиться). Этот метод непригоден для проектирования КИХ-фильтров. Лучшие результаты дает метод, основанный на использовании весовой последовательности конечной длины w(n), называемой окном, для модификации коэффициентов Фурье в формуле (3.9) с тем, чтобы управлять сходимостью ряда Фурье. Для большинства приемлемых окон преобразование Фурье последовательностиw(n)имеет главный лепесток, содержащий почти всю энергию окна, и боковые лепестки, которые обычно быстро затухают. Чтобы получить КИХ-аппроксимацию функции, формируется последовательность, равная нулю за пределами интервала. Поскольку результирующая характеристика фильтра равна свертке идеальной частотной характеристики и частотной характеристика окна, то ширина переходных полос зависит от ширины главного лепестка функции. Кроме того, на всех частотахвозникают ошибки аппроксимации, имеющие вид пульсаций частотной характеристики, которые обусловлены боковыми лепестками функции. Из приведенного выше следует, что оптимальное окно должно иметь во-первых, минимальную ширину главного лепестка частотной характеристики, во-вторых, минимальную площадь под боковыми лепестками. К сожалению, эти два требования несовместимы и необходим компромиссный вариант. Ких-фильтры. Методы синтеза. Класс последовательностей конечной длины (КИХ-последовательности) обладает некоторыми свойствами, желательными с точки зрения построения фильтров. Например, никогда не возникает вопрос об устойчивости и физической реализуемости фильтров, поскольку КИХ-последовательности гарантируют устойчивость. Более того, КИХ-последовательности можно выбрать так, чтобы фильтры имели строго линейные фазовые характеристики. Поэтому, используя КИХ-последовательности, можно проектировать фильтры с произвольной амплитудной характеристикой. Существуют три основных метода синтеза КИХ-фильтров: • метод взвешивания (метод «окна»); • метод частотной выборки; • метод оптимальных фильтров. Свойства КИХ-фильтров Имеется много причин, побуждающих к изучению способов проектирования КИХ-фильтров. Основными достоинствами этих фильтров являются: 1. Легко создавать КИХ-фильтры со строго линейной фазовой характеристикой (постоянной групповой задержкой). Во многих случаях, когда проектируется фильтр с произвольной амплитудной характеристикой, это упрощает задачу аппроксимации. 2. КИХ-фильтры, реализуемые нерекурсивно, т.е. с помощью свертки, всегда устойчивы. 3. При нерекурсивной реализации КИХ-фильтров шумы округления, возникающие за счет выполнения арифметических операций с конечной точностью, легко минимизировать. 4. КИХ-фильтры можно эффективно реализовывать при помощи методов быстрой свертки, основанных на применении алгоритма БПФ. Недостатки КИХ-фильтров следующие: 1. Для аппроксимации фильтров, частотные характеристики которых имеют острые срезы, требуется импульсная характеристика с большим числом отсчетов N. Следовательно, возрастает объем вычислительных операций. 2. Задержка в КИХ-фильтрах с линейной фазовой характеристикой не всегда равна целому числу интервалов дискретизации. Характеристики КИХ-фильтров с линейной фазовой характеристикой Пусть {h(n)} – физически реализуемая последовательность конечной длины, заданная на интервале. Ееz-преобразование равно . (3.1) Преобразование Фурье от {h(n)} (3.2) является периодическим по частоте с периодом , т.е. . (3.3) Рассматривая только действительные последовательности {h(n)}, получим дополнительные ограничения на функцию, представив ее через амплитуду и фазу: . (3.4) Потребуем при расчете КИХ-фильтров строго линейной фазовой характеристики, и рассмотрим, при каких условиях импульсная характеристика фильтра h(n)будет это обеспечивать. Требование линейности фазовой характеристикиимеет вид , (3.5) где - постоянная фазовая задержка, выраженная через число интервалов дискретизации. Используя (3.4) и (3.5) соотношение (3.2) переписывается в виде: . (3.6) Приравнивая действительные и мнимые части, и деля друг на друга правые и левые части полученных равенств, можно получить уравнение, решением которого будут следующие значения: , (3.7) . (3.8) Смысл их заключается в следующем. Условие (3.7) означает, что для каждого Nсуществует только одна фазовая задержка, при которой может достигаться строгая линейность фазовой характеристики фильтра. Из условия (3.8) следует, что при заданном, удовлетворяющем условию (3.7), импульсная характеристика должна обладать симметрией. Рассмотрим использование условий (3.7) и (3.8) для случаев четного и нечетногоN. ЕслиN- нечетно, то задержка в фильтре равна целому числу интервалов дискретизации. Типичная импульсная характеристика фильтра с линейной фазой для случаяN=11 приведена на рис. 3.1. Типичная импульсная характеристика фильтра с линейной фазой при четномNпоказана на рис. 3.2. ЗдесьN=10. Рис.3.1. N – нечетно. Рис.3.2. N – четно. Расчет КИХ-фильтров методом взвешивания Как было сказано ранее, частотную характеристику любого цифрового фильтра можно представить рядом Фурье: , (3.9) где . (3.10) Видно, что коэффициенты Фурье совпадают с коэффициентами импульсной характеристики цифрового фильтра. Использование этих соотношений для проектирования КИХ-фильтра связано с двумя трудностями. Во-первых, импульсная характеристика фильтра имеет бесконечную длину, поскольку суммирование в (3.9) производится в бесконечных пределах. Во-вторых, фильтр физически нереализуем, так как импульсная характеристика начинается в , т.е. никакая конечная задержка не сделает фильтр физически реализуемым. Один из возможных методов получения КИХ-фильтра, аппроксимирующего заданную функцию , заключается в усечении бесконечного ряда Фурье (3.9) за. Однако простое усечение ряда приводит к явлению Гиббса, которое проявляется в виде выбросов и пульсаций до и после точек разрыва в аппроксимируемой частотной характеристики. Причем, максимальная амплитуда пульсаций частотной характеристики не уменьшается с увеличением длины импульсной характеристики, т.е. учет все большего числа членов ряда Фурье не приводит к уменьшению максимальной амплитуды пульсаций. Вместо этого уменьшается ширина выброса. Поэтому простое усечение ряда Фурье (3.9) не приводит к приемлемой аппроксимации идеального фильтра нижних частот (к чему необходимо стремиться). Этот метод непригоден для проектирования КИХ-фильтров. Лучшие результаты дает метод, основанный на использовании весовой последовательности конечной длины w(n), называемой окном, для модификации коэффициентов Фурье в формуле (3.9) с тем, чтобы управлять сходимостью ряда Фурье. Для большинства приемлемых окон преобразование Фурье последовательностиw(n)имеет главный лепесток, содержащий почти всю энергию окна, и боковые лепестки, которые обычно быстро затухают. Чтобы получить КИХ-аппроксимацию функции, формируется последовательность, равная нулю за пределами интервала. Поскольку результирующая характеристика фильтра равна свертке идеальной частотной характеристики и частотной характеристика окна, то ширина переходных полос зависит от ширины главного лепестка функции. Кроме того, на всех частотахвозникают ошибки аппроксимации, имеющие вид пульсаций частотной характеристики, которые обусловлены боковыми лепестками функции. Из приведенного выше следует, что оптимальное окно должно иметь во-первых, минимальную ширину главного лепестка частотной характеристики, во-вторых, минимальную площадь под боковыми лепестками. К сожалению, эти два требования несовместимы и необходим компромиссный вариант. Под синтезом ЦФ понимают задачу определения алгоритма его работы, удовлетворяющего заданным условиям. Эти условия могут формулироваться как во временной области, так и в частотной. Обычно разделяют задачи синтеза рекурсивных и нерекурсивных ЦФ. Рассмотрим вначале задачи синтеза рекурсивных ЦФ по заданный временным характеристикам. Синтез ЦФ, обладающих заданными характеристиками, может осуществляться двумя способами: 1) посредством построения ДЛС, эквивалентной в некотором смысле некоторой НЛС с заданными характеристиками; причем синтез НЛС осуществляется с помощью теории синтеза аналоговых фильтров; 2) путем построения цифровой системы, обладающей характеристиками, близкими в некотором смысле к заданным. Вся процедура синтеза в этом случае проводится в дискретном ряде точек временной оси и не требует знания непрерывного фильтра прототипа. Часто синтез ЦФ осуществляют по методу замены НДС её импульсным эквивалентом,в частности, по заданной импульсной реакции НЛС, по заданной переходной функции НЛС и др. Здесь входную последовательность {xi}, рассматриваемую как результат дискретизации некоторого аналогового сигнала, пропускают через фильтр-интерполятор, осуществляющий дискретно-аналоговое преобразование. Затем полученный аналоговый сигнал пропускают через заданную НЛС. Систему, состоящую из интерполирующего фильтра (ИФ) и НЛС, называют приведенной непрерывной системой. Её отклик в дискретные моменты времени определяется как дискретная свертка исходной последовательности {xi} и импульсной реакции приведенной НЛС. Эта процедура, следовательно, позволяет получить ДЛС по исходной НЛС. При этом различному типу ИФ соответствует свой метод синтеза. Если в качестве ИФ использован безынерционный усилитель, то метод называется z-преобразованием или методом инвариантной импульсной реакции. Здесь импульсную реакцию искомой ДЛС (ЦФ) определяют в результате дискретизации импульсной реакции НЛС. Затем, вычисляя Z-преобразование от последней, находят передаточную функцию ЦФ в виде дробно-рациональной функции, от которой переходят к разностному уравнению, задающему алгоритм работы ЦФ. В случае, когда ИФ представляет собой ступенчатый интерполятор с импульсной реакцией, постоянной на интервале дискретизации Т, приходят к синтезу ЦФ методом инвариантной переходной характеристики (метод Цыпкина-Гольденберга). Если же ИФ - линейный интерполятор, то метод синтеза ЦФ называют методом Рагаззини-Бергена. Перейдем к вопросам синтеза рекурсивных ЦФ в частотной области. Синтез ЦФ с близкими к заданным АЧХ и ФЧХ аналогового прототипа затруднителен и может быть выполнен лишь с помощью специальных методов оптимизации, требующих громоздких вычислений. Поэтому прибегают к более простым, но конструктивным методам синтеза, не требующим знания всех характеристик НЛС. В частности, для синтеза ЦФ по заданной АЧХ НЛС широко используется метод билинейного преобразования. Суть метода состоит в переходе от переменной р в передаточной функции К(р) НЛС к переменной z, определяющей передаточную функцию H(z) синтезируемого ЦФ. Для того, чтобы получить возможность непосредственного перехода от К(р) к Н(z), необходимо аппроксимировать функцию р = lnz/T рациональной функцией. Это можно сделать, если воспользоваться разложением этой функции в ряд Тейлора. Ограничиваясь только первым членом этого разложения, получают зависимость вида р = k(z - 1)/(z + 1), k = 2/T. Данное преобразование называют билинейным и находит широкое применение при проектировании стандартных ФНЧ, ФВЧ и полосовых фильтров. Основные особенности этого метода состоят в неучете ФЧХ аналогового прототипа и в наблюдаемой деформации масштаба частот ЦФ по отношению НЛС. Поэтому данный метод в основном применяется тогда, когда ЦФ может иметь произвольную ФЧХ. Второй эффект влияет незначительно при уменьшении величины T. Синтез нерекурсивных ЦФ также осуществляется как во временной, так и в частотной областях. Наиболее простым способом синтеза ЦФ является синтез по известной импульсной реакции НДС. С этой целью последнюю дискретизируют и затем усекают до конечного числа отсчетов N так, что получается импульсная реакция ЦФ конечной длительности. Алгоритм работы такого ЦФ описывается дискретной сверткой исходной сигнальной последовательности с полученной импульсной реакцией. Заметим, что такой же результат получается при методе инвариантной импульсной реакции, применяемом для синтеза рекурсивных ЦФ. При проектировании нерекурсивных ЦФ методом усечения дискретизированной импульсной реакции НЛС возникает вопрос о выборе весовой функции усечения (временного окна). Дело в том, что при простом прямоугольном окне в спектре отклика ЦФ наблюдаются значительные пульсации, искажающие спектр и сам отклик. В теории рядов Фурье это явление известно под названием эффекта Гиббса. Важно подчеркнуть, что уровень пульсаций практически не зависит от длительности N окна. Возникает вопрос о том, какими свойствами должно обладать временное окно, чтобы уровень пульсаций по возможности был минимальным. В настоящее время предложено много различных окон, решающих в той или иной мере эту задачу. Известны окна Бартлета, Тыэки, Хэмминга. Нерекурсивный ЦФ или КИХ-фильтр может обладать строго линейной ФЧХ. Для этого достаточно, чтобы его импульсная реакция была симметрична относительно точки (N - 1)/2 при нечетном N. Этого нельзя достичь в БИХ-фильтре, так как из-за бесконечной протяженности импульсной реакции она не может быть строго симметричной. Для синтеза нерекурсивного ЦФ в частотной области можно использовать метод частотной выборки. Для этого в H(z) подставляют z = еJωT, где частота ω принадлежит непрерывному множеству, ограниченному частотой дискретизации ω = 2πfd. Для перехода к дискретным значениям частоты в спектральной области осуществляется дискретизация с шагом по частоте Δω = 2π/NT = ωd/N. В результате находят комплексные амплитуды Gk = H(ejk2π/N) = H(WNk) спектра импульсной реакции синтезируемого ЦФ. Взяв от последовательности {Gk} обратное ДПФ или применив алгоритм БПФ, находят значения импульсной реакции (gj). Синтез нерекурсивных ЦФ методом частотной выборки особенно эффективен в случае узкополосных частотно-избирательных систем. Выбирая N достаточно большим, можно получить ЦФ с любыми требуемыми характеристиками. Однако такой путь не всегда ведет к наиболее простым в реализации ЦФ. Поэтому на практике применяют и другие более совершенные методы синтеза нерекурсивных ЦФ в частотной области. Часто нерекурсивная фильтрация в частотной области осуществляется на основе использования алгоритма БПФ. Здесь последовательность исходного сигнала {xi} регистрируется в блоке памяти ЦФ. Затем с помощью прямого БПФ определяется её дискретный спектр {Xk}. Умножением этого спектра на дискретную функцию {Gk} комплексного коэффициента передачи синтезируемого ЦФ находится дискретный спектр отклика Yk = NXkGk, k = 0, 1, ..., N-1. Применяя теперь обратное БПФ к {Yk}, находится сам отклик как дискретная последовательность {yi} . Для повышения эффективности этого метода синтеза ЦФ применяют процедуру перекрытия спектров с накоплением. КОНТРОЛЬНЫЕ ВОПРОСЫ: 1. Что понимается под задачей синтеза цифрового фильтра? 2. Поясните задачи синтеза рекурсивного ЦФ по заданным временным характеристикам. 3. Как осуществляется синтез ЦФ по методу замены НЛС её импульсным эквивалентом? 4. Поясните задачу синтеза ЦФ методом билинейного преобразования. 6. В чем суть задачи синтеза ЦФ методом частотной выборки? ЛЕКЦИЯ 10. АНАЛОГОВЫЕ ФИЛЬТРЫ, ИХ ТИНЫ И ХАРАКТЕРИСТИКИ. МЕТОДЫ АППРОКСИМАЦИИ АНАЛОГОВЫХ ФИЛЬТРОВ. При измерении параметров ТП в распределенных АСОИУ возникают различные ошибки. Причинами ошибок измерения, влияющих на качество управления ТП, являются: – нелинейность функций передаточных звеньев (из-за механических ударов, трения, гистерезисных явлений, старения, загрязнения); – случайные помехи (наводки внешних полей, импульсы помех, обрыв линий связи); – инерционность сигналов, изменяющихся во времени (задержки, зоны нечувствительности). Ошибки измерения параметров ТП можно разделить на две группы: а) систематические (они приводят к изменению полезного сигнала по амплитуде и фазе); б) случайные (приводят к наложению сигналов помех, шумов). Первый вид ошибок устраняется сравнением с эталоном или моделью, калибровкой, линеаризацией, корректировкой измеренных значений; Второй – проверкой данных на достоверность, фильтрацией, выбором соответствующей частоты опроса датчиков. Рассмотрим подробнее аналоговую фильтрацию сигналов, типы фильтров, их параметры и расчет. Виды аналоговых фильтров и их характеристики Фильтрация сигналов играет важную роль в цифровых системах управления. В них фильтры используются для устранения случайных ошибок измерения (наложения сигналов помех, шумов) (рис. 16.1). 16.1. Фильтрация зашумленного синусоидального сигнала Различают аппаратную (схемную) и цифровую (программную) фильтрацию. В первом случае используют электронные фильтры из пассивных и активных элементов, во втором случае применяют различные программные методы выделения и устранения помех. Аппаратная фильтрация применяется в модулях УСО (устройств связи с объектом) контроллеров и распределенных систем сбора данных и управления. Цифровая фильтрация используется в УВМ верхнего уровня АСУ ТП. В первую очередь, нас интересуют частотно-избирательные фильтры, которые пропускают сигналы определенных частот и задерживают, ослабляют сигналы других частот. Существуют следующие типы аналоговых фильтров: – фильтры нижних частот (ФНЧ), которые пропускают низкие частоты и задерживают высокие частоты; – фильтры верхних частот (ФВЧ), которые пропускают высокие частоты и задерживают низкие частоты; – полосно-пропускающие фильтры (ФПП), которые пропускают полосу частот и задерживают частоты, расположенные выше и ниже этой полосы; – полосно-заграждающие фильтры (ФПЗ), которые задерживают полосу частот и пропускают частоты, расположенные выше и ниже этой полосы. Принцип действия ФНЧ лучше всего проиллюстрировать на примере фильтрации синусоидального сигнала (рис. 16.1), на который наложен высокочастотный шум. Как видно из этого рисунка, сигнал после фильтрации шумов имеет существенный фазовый сдвиг, а амплитуда сигнала практически не изменилась. Передаточная функция (ПФ) фильтра описывается выражением H(s) = V2(s) / V1(s), (16.1) где V1 и V2 – входное и выходное напряжения фильтра. Для s = j можно записать Н( j) = |Н( j)| exp [ j ()], (16.2) где |Н( j)| – модуль ПФ или АЧХ;  () – ФЧХ;  – угловая частота (рад/с), связанная с частотой f (Гц) отношением  = 2 f. Диапазоны или полосы частот, в которых сигналы проходят, – это полосы пропускания и в них значение АЧХ |Н(j)| велико, а в идеальном случае постоянно. Диапазоны частот, в которых сигналы подавляются, – это полосы задерживания и в них значение АЧХ |Н( j)| мало, а в идеальном случае равно нулю. АЧХ фильтров различного вида приведены на рис. 16.2-16.5. Рис. 16.2 Рис. 16.3 Рис. 16.4 Рис. 16.5 На этих АЧХ параметр с – это частота среза, при которой заканчивается полоса пропускания фильтра, а з – частота задерживания, с которой начинается полоса задерживания. Полоса пропускания – это диапазон частот, где значение АЧХ больше заданной величины А1. Полоса задерживания – это диапазон частот, в котором АЧХ меньше значения – А1. В реальных фильтрах с неидеальной АЧХ интервал частот перехода от полосы пропускания к полосе задержания, называют переходной областью. Затухание в децибелах (дБ) определяется из АЧХ по формуле a = –20 log10 |Н( j)|. (2.3) Значению амплитуды А = 1 соответствует затухание a = 0. Если A1 = A/= 1/, то затухание на частоте с: а1 = –20 log10 (1/) = 10 log10 2 = 3 дБ. Обычно принимают а1 = 0,1; 0,5; 1; 2 или 3 дБ. АЧХ и затухание идеального и реального фильтров показаны на рис. 16.6 и 16.7. Рис. 16.6. Идеальная и реальная АЧХ ФНЧ Рис. 16.7. Идеальная и реальная кривые затухания ФНЧ В отличие от ФНЧ всепропускающий фильтр имеет почти постоянное значение АЧХ для всех частот (т.е. сигналы всех частот проходят одинаково хорошо), а ФЧХ является функцией частоты. Этот фильтр называется фазосдвигающим, так как его ФЧХ может изменяться, или сдвигаться, в зависимости от частоты. ФЧХ фильтра линейна и определяется из выражения: , (16.4) где  – постоянное число. Чем более нелинейна ФЧХ, тем сильнее искажается выходной сигнал. Чем лучше АЧХ, тем хуже ФЧХ, и наоборот. При расчете фильтра нужно найти компромисс между хорошими АЧХ и ФЧХ. Время замедляющий фильтр В нем основная характеристика – время замедления Т() – отрицательное значение наклона ФЧХ: Т() = – d  ()/d . (16.5) Для линейной ФЧХ из (2.4) Т() = , где  – постоянное число. Фильтр рассчитывается так, что Т() почти постоянно для заданного диапазона частот. Передаточная функция реализуемого фильтра описывается выражением (16.6) где а и b – постоянные величины, а m, n = 1, 2, 3 ... (m  n). (16.7) Степень полинома знаменателя n определяет порядок фильтра. Чем он выше, тем лучше АЧХ, но сложнее схема, а стоимость выше. При разработке фильтров нужно получить АФХ, ближе к идеальной с заданной степенью точности при наименьших затратах. Типы аналоговых фильтров. Различают пассивные и активные аналоговые, не интегральные фильтры. Пассивные фильтры создаются на основе пассивных R, L, C элементов. Однако катушки индуктивности не приспособлены для интегрального исполнения. На низких частотах (ниже 0,5 МГц), параметры катушек индуктивности неудовлетворительны из-за больших размеров и отклонения характеристик от идеальных. Схемы пассивных фильтров приведены на рис. 16.8, 16.9. Активные фильтры создаются на основе R, C элементов и активных элементов – операционных усилителей, которые должны иметь: • высокий коэффициент усиления (в 50 раз больше, чем у фильтра); • высокую скорость нарастания выходного напряжения (до 100-1000 В/мкс). Схемы активных ФНЧ первого и второго порядков показаны на рис. 16.10 и 16.11. Построение фильтров n-го порядка осуществляется каскадным соединением звеньев N1, N2, ... , Nm с ПФ Н1(s), H2(s), ..., Нm(s). Фильтр четного порядка с п > 2 содержит n/2 звеньев второго порядка, соединенных каскадно. Рис. 16.8. Пассивный ФНЧ первого порядка (а) и его АЧХ (б) Пассивные ФНЧ второго порядка и выше могут иметь Т- и П-образные схемы. Рис. 16.9. Пассивные ФНЧ: а) Т-образный; б) П-образный; Рис. 16.10. Активный ФНЧ первого порядка Рис. 16.11 Активный ФНЧ второго порядка Фильтр нечетного порядка с п > 2 содержит (п – 1)/2 звеньев второго порядка и одно звено первого порядка с ПФ (16.8)-(16.9). Передаточные функции для фильтров первого порядка имеют вид (16.8) где С – постоянное число; P(s) – полином первой или нулевой степени. Передаточные функции для фильтров второго порядка таковы (16.9) где В и С – постоянные числа; P(s) – полином второй или меньшей степени. Типы аналоговых фильтров низкой частоты У ФНЧ максимальное затухание в полосе пропускания a1 не превышает 3 дБ (A1 = 0,707), а типовое значение затухания в полосе задерживания a2 больше и находится в пределах от 20 до 100 дБ (0,1 ≥ A2 ≥ 0,00001). Коэффициент усиления фильтра НЧ – это значение его передаточной функции при s = 0 или значение его АЧХ при  = 0, т.е. он равен А. В зависимости от вида АЧХ и ПФ различают следующие типы ФНЧ: • Баттерворта – с монотонной АЧХ (рис. 16.12); • Чебышева (типа I) – АЧХ содержит пульсации в полосе пропускания и монотонна в полосе задерживания (рис. 16.13; • инверсные Чебышева (типа II) – АЧХ монотонна в полосе пропускания и обладает пульсациями в полосе задерживания (рис. 16.14; • эллиптические – АЧХ имеет пульсации и в полосе пропускания и в полосе задерживания (рис. 16.15). Рис. 16.12. АЧХ фильтра Баттерворта Рис. 16.13. АЧХ фильтра Чебышева I типа Рис. 16.14. АЧХ фильтра Чебышева II типа Рис. 16.15. АЧХ эллиптического фильтра Фильтры Баттерворта АЧХ фильтра низких частот n-го порядка (16.10) где n = 1,2,3… . АЧХ монотонно спадает при росте частоты. Увеличение порядка n ведет к улучшению характеристики, как это показано на рис. 16.16. Рис. 16.16. АЧХ фильтров Баттерворта n–го порядка Передаточная функция фильтра Баттерворта как полиномиального фильтра равна , (16.11) где К – постоянное число. Для нормированного фильтра Баттерворта, т. е. при с = 1 рад/с ПФ для п = 2, 4, 6 … в виде произведения сомножителей равна . (16.12) Для п = 3, 5, 7 … передаточная функция нормированного фильтра Баттерворта равна . (16.13) Коэффициенты ak и bk задаются при b0 = 1 и k = 1, 2 ... выражениями: (16.14) Коэффициент усиления К фильтра в (16.11) равен произведению коэффициентов усиления звеньев Аk и/или А0 в (16.12) или (16.13). Фильтры Чебышева I типа Эти фильтры имеют АЧХ следующего вида: (16.15) где параметры  и К – постоянные числа, а Сп – полином Чебышева первого рода степени п, равный . (16.16) Аргумент x в (16.16) для фильтра Чебышева равен отношению /с . Фильтр Чебышева называют равноволновым, так как все пульсации равны по значению. Для К = 1 размах пульсаций в полосе пропускания (ripple passband) . (16.17) Размах Rр можно уменьшить, выбрав значение параметра  достаточно малым. Минимально допустимое затухание в полосе пропускания – постоянный размах пульсаций – выражается в децибелах как . (16.18) Откуда . (16.19) Рис. 16.17. АЧХ фильтров Чебышева n–го порядка Рис. 16.18. ФЧХ фильтров Чебышева и Баттерворта n–го порядка Сравнение фильтров Баттерворта и Чебышева. ПФ фильтров НЧ Чебышева и Баттерворта идентичны по форме, описываются выражениями (16.12)-(16.13). АЧХ фильтра Чебышева лучше АЧХ фильтра Баттерворта такого же порядка, так как у первого уже ширина переходной области. Однако у фильтра Чебышева ФЧХ хуже (более нелинейна), чем ФЧХ у фильтра Баттерворта. Выбор минимального порядка n фильтра осуществляют по заданным значениям: – максимально допустимого затухания в полосе пропускания a1 (дБ); – минимально допустимого затухания в полосе задерживания a2 (дБ); – частоты среза с (рад/с) или fc (Гц); – максимальной допустимой ширины переходной области TW , равной TW = 1– с. (16.20) Для фильтра Баттерворта с a1 = 3 дБ минимальный порядок п из (1.18) равен . (16.21) Для фильтра Чебышева с К = 1 минимальный порядок п из (1.18) равен , (16.22) где arch (читается ареакосинус) – обратная функция от гиперболического косинуса ch x = (ex + e–х) / 2 . Из (16.20) получаем 1/с = (TW /с) + 1. Пример. Пусть задано: a1 = 3 дБ, a2 = 20 дБ, fс = 1000 Гц, TW = =300 Гц. Тогда 1/с = (300/1000) + 1 =1,3. Из (16.21) для фильтра Баттерворта . (16.23) Из (16.22) для фильтра Чебышева . Как видим, при заданных параметрах фильтр Чебышева типа имеет в два раза меньший порядок, чем фильтр Баттерворта. Фильтры Чебышева II типа. Эти фильтры называются также инверсные фильтрами Чебышева и имеют АЧХ следующего вида: (16.24) где параметр  – постоянное число, а Сп – полином Чебышева первого рода степени п (16.16). АЧХ инверсного фильтра Чебышева показана на рис. 16.19. Рис. 16.19. АЧХ фильтра Чебышева II типа Частота среза с по уровню 3 дБ равна (16.25) АЧХ фильтра Чебышева монотонна в полосе пропускания и имеет пульсации в полосе задерживания, значения которые равны (16.26) Величина  при а2 = –20 log10 A2 определяется выражением (16.27) Расчет параметров фильтров. Для заданных значений – порядка п , затухания а2 и частоты 1 можно из (16.27) найти значение  , а из (16.24) – требуемую АЧХ. Тогда частоту wс и ширину TW можно определить из (16.25) и (16.26). Можно установить wс , что легче, и из (16.25) найти частоту 1. При одинаковых параметрах фильтры Чебышева I и II типов имеют равный порядок, но меньший порядка фильтра Баттерворта. ПФ инверсного фильтра. Для четного порядка п она равна (16.28) Для нечетного порядка п она равна (16.29) Так как ПФ не полиномиальна, инверсный фильтр сложнее в реализации. При выборе типа фильтра следует руководствоваться следующими соображениями. 1. Если требуется монотонная АЧХ в полосе пропускания, то инверсный фильтр Чебышева по параметрам лучше фильтра Баттерворта того же порядка. 2. Если допустимы пульсации в полосе пропускания, то лучше фильтр Чебышева, так как его ПФ проще, чем у инверсного фильтра Чебышева. 3. Если нужна монотонная АЧХ, тогда хорош и фильтр Баттерворта, так как его ПФ проще, чем у инверсного фильтра. Эллиптические фильтры. Эти фильтры имеют АЧХ с пульсациями в полосе и пропускания и задерживания. Для заданных порядка и допустимых отклонений АЧХ в этих полосах они обладают самой узкой шириной переходной области. АЧХ эллиптического фильтра n = 5 показана на рис. 16.20. Рис. 16.20. АЧХ эллиптического фильтра Пульсации в полосе пропускания и в полосе задерживания у эллиптического фильтра равны по значению между собой, но не обязательно равны размаху друг другу. Пульсации в полосах характеризуют: • максимальным допустимым затуханием – неравномерностью передачи в полосе пропускания (Rp – от англ. ripple passband), дБ, Rp = a1 = – 20 log A1. • минимальным затуханием – неравномерностью передачи в полосе задерживания (Rs – от англ. ripple stopband), дБ, Rs = a2 = – 20 log A2. Для заданных Rp и Rs повышение порядка n ведет к увеличению числа пульсаций в полосах пропускания и задерживания и уменьшению ширины переходной области TW. Зависимость порядка эллиптических фильтров и фильтров Чебышева от ширины переходных областей (Rp = 0,1 и Rs = 60 дБ) приведена на рис. 16.21. Рис. 16.21. Зависимость порядка фильтров эллиптических (1) и Чебышева (2) от ширины переходных областей TW при Rp = 0,1 и Rs = 60 дБ Как видно из рис. 16.21 при ширине TW /с нормированной переходной области не более 0,1, а 1/с – меньше или равным 1,1 нужен эллиптический фильтр с порядком п = 10 или фильтр Чебышева с порядком п = 22. Преимущество эллиптических фильтров над фильтрами Чебышева более заметно при узкой ширине TW. Так, если ширина TW /с = 0,03, то нужен эллиптический фильтр 12-го порядка или фильтр Чебышева уже 39-го порядка. Передаточная функция эллиптического фильтра идентична функции инверсного фильтра Чебышева, описываемой (4.5) и (4.6). Однако постоянные аi , bi и сi находят с помощью эллиптических функций Якоби (эллиптических интегралов, синуса и косинуса, дельта-амплитуды). Расчет эллиптических фильтров, как и расчет инверсных фильтров Чебышева сложен и его обычно ведут с помощью табулированных функций или на ЭВМ с помощью специальных программ – Classic Control, MATLAB. В последнем программном продукте имеется специальные пакеты программ для расчета фильтров различного типа высокого порядка. Все пропускающие фильтры. Если в частотно-избирательных фильтрах наиболее важной характеристикой является АЧХ, то во все пропускающих фильтрах – ФЧХ и/или соответствующее ей время замедления. Все пропускающий или фазосдвигающий фильтр обладает постоянной АЧХ |Н(j)|= K и ФЧХ (), являющейся функцией частоты (рис. 16.22). Линейная ФЧХ все пропускающих фильтров характеризуется постоянным временем замедления Т(w) = – d j(w)/d w = t. Такой фильтр является фильтром с линейной фазой или с постоянным временем замедления. Наилучшим из полиномиальных фильтров с постоянным временем замедления является фильтр Бесселя. Рис. 16.22. ФЧХ все пропускающего фильтра Фильтр Бесселя имеет ПФ, описываемую выражением (16.30) где К – коэффициент усиления фильтра, а Bn(s) – полином n-й степени (16.31) где для k = 0, 1, 2, ..., n (16.32) Полином Bn(s) при с = 1 относится к полиномам Бесселя, отсюда название этого фильтра. ФЧХ фильтров Бесселя намного лучше ФЧХ фильтров Баттерворта и Чебышева. Однако АЧХ фильтра Бесселя хуже АЧХ фильтров Баттерворта или Чебышева. При увеличении порядка фильтра время замедления Т() приближается к постоянному значению (рис. 16.23). Рис. 16.23. ФЧХ фильтра Бесселя в зависимости от порядка фильтра АЧХ фильтра Бесселя монотонно спадает от максимального значения на нулевой частоте. Она имеет сходство с АЧХ фильтра Баттерворта, но крутизна нарастания затухания гораздо меньше. Частота с для фильтра Бесселя в (16.32) – это не частота среза, а частота, определяющая диапазон постоянного времени замедления. Для упрощения расчета перечисленных выше фильтров в справочных таблицах приведены значения различных параметров, входящих в передаточные функции фильтров для заданных значений – порядка п, затуханий a1 и a2, частоты с , ширины переходной полосы TW. По этим значениям находят номиналы резисторов и конденсаторов, образующих различные фильтры в соответствии с выбранной электрической схемой. КОНТРОЛЬНЫЕ ВОПРОСЫ 1. Поясните понятие - время замедляющий фильтр 2. Изобразите график ФЧХ фильтра Бесселя в зависимости от порядка фильтра 3.Дайте характеристику - АЧХ фильтра Бесселя 4.Дайте определение понятию эллиптические фильтры ЛЕКЦИЯ 11. ИЗМЕНЕНИЕ ЧАСТОТЫ ДИСКРЕТИЗАЦИИ СИГНАЛА. ПОСТАНОВКА ЗАДАЧИ ИНТЕРПОЛЯЦИИ Сочленение цифровых устройств осуществляется достаточно просто, если частоты дискретизации в коммутируемых устройствах совпадают. В этом случае достаточно согласовать законы компандирования в разных звеньях тракта. Реализуется это на основе ИКМ-сигнала с равномерным квантованием, который, как правило, является исходным для всех АЦП. По этому общему для всех цифровых устройств закону и объединяют тракты системы с разными законами компандирования. Если же частоты дискретизации не совпадают, то интерфейс должен содержать также дополнительно блок передискретизации. Эта задача, например, возникает при организации стыка аппаратуры ТФП (fд = 48 кГц) с аппаратурой ТПРП (fд = 32 кГц) или при соединении аппаратуры систем записи сигналов (fд = 44,1 кГц) с аппаратурой ТФП и каналов связи. Изменить частоту дискретизации можно двумя способами: промежуточным переходом к аналоговой форме и непосредственно путем преобразований в цифровом сигнале. Первый способ предполагает преобразование сигнала в аналоговую форму и его последующую дискретизацию с новой частотой и кодирование в новом формате. При этом способе исходный ЗС дважды подвергается квантованию, что приводит к возрастанию мощности шума квантования и не является желательным. Например, при переходе от ТФП (где fд = 48 кГц и т = 16) к ТПРП и ТВРП (где fд = 32 кГц и т = 10) уменьшение числа разрядов приводит к возрастанию шумов квантования на 20…26 дБ. Собственно же передискретизация в этом случае не сказывается на общем отношении сигнал/шум. Если передискретизация выполняется без преобразования цифрового сигнала в аналоговую форму, то для получения новой совокупности отсчетов используется процедура интерполяции. При этом реализация метода наиболее проста, если частоты дискретизации имеют общий кратный сомножитель. В качестве примера рассмотрим переход от частоты дискретизации fд = 48 кГц к новому значению, равному 32 кГц (рис. 2.5). Наименьшее значение частоты, кратное этим величинам, составляет 96 кГц. Увеличение частоты дискретизации достигается здесь интерполяцией промежуточных дополнительных значений отсчетов, а уменьшение частоты дискретизации до 32 кГц — отбрасыванием лишних значений. Итак, поскольку 32 кГц = (482):3, то процедура передискретизации состоит в следующем: 1. Вначале выполняется преобразование от fд = 48 кГц к fд = 96 кГц. Для этого методом интерполяции рассчитывают дополнительные значения отсчетов в середине исходного тактового интервала (показаны на рис. 2.5, а штриховой линией). Полученная в результате этого суммарная последовательность отсчетов соответствует fд = 96 кГц. Рис. 2.5 — К пояснению процесса передискретизации цифрового сигнала 2. Далее из этой новой последовательности выбирается каждый третий отсчет (рис. 2.5, б), что означает снижение частоты дискретизации в 3 раза по отношению к fд = 96 кГц. Для большей наглядности на рисунке показаны не кодовые слова, а соответствующие им отсчеты ЗС. В действительности же интерфейс оперирует, выполняя интерполяцию с цифровыми сигналами. Осуществляя передискретизацию сигналов, следует помнить о возможности появления искажений, вызванных так называемым наложением спектров. Для устранения этого явления при уменьшении частоты дискретизации в интерфейс в случае необходимости должны включаться цифровые фильтры для ограничения полосы частот исходного ЗС. Например, необходимость такой дополнительной фильтрации возникает при уменьшении fд с 48 до 32 кГц, когда требуется снизить полосу частот ЗС с 20 кГц (ТФП) до 15 кГц (ТПРП и ТВРП). Простейшая задача интерполяции заключается в следующем. На отрезке [a, b] заданы n + 1 точки xi = х0, х1, . . ., хn, которые называются узлами интерполяции, и значения некоторой функции f(x)  в этих точках f(x0) = y0, f(x1) =  y1,  . . ., f(xn) = yn. (1) Требуется построить функцию  (х) (интерполяционная функция), принадлежащую известному классу и принимающую в узлах интерполяции те же значения, что и f(x), т. е. такую, что (x0) = y0, (x1) =  y1,  . . ., (xn) = yn. (2) Геометрически это означает, что нужно найти кривую y =  (х) некоторого определенного типа, проходящую через заданную систему точек M(xi, yi) (i = 0, 1, ..., n) (Рисунок 1).   Рисунок 1 В такой общей постановке задача может иметь бесконечное множество решений или совсем не иметь решений. Однако эта задача становится однозначной, если вместо произвольной функции  (х) искать полином  (х) (интерполяционный полином) степени не выше n, удовлетворяющий условиям (2), т. е. такой, что  (x0) = y0,  (x1) =  y1,  . . .,  (xn) = yn. (3) Полученную интерполяционную формулу (4) обычно используют для приближенного вычисления значений данной функции  (х) для значений аргумента х, отличных от узлов интерполяции. Такая операция называется интерполяцией функций. Различают два вида интерполяции: 3) Интерполяционные формулы Ньютона До сих пор не делалось никаких предположений о заданных значениях аргумента, которые могли быть совершенно произвольными. Предположим дополнительно, что рассматриваемые значения аргумента являются равноотстоящими, т. е. образуют арифметическую прогрессию. В этом случае шаг таблицы h = хi+ - xi (i = , , , ..., n) = const яaeяется величиной постоянной. Для таких таблиц построение интерполяционных формул (как впрочем, и вычисление по этим формулам) заметно упрощается. Прежде чем перейти к рассмотрению этого вопроса, познакомимся с понятием конечных разностей. Конечные разности Пусть функция задана таблицей с постоянным шагом. Разности между значениями функции в соседних узлах интерполяции называют конечными разностями первого порядка:  yi = yi+ - yi (i = , , , ...). Из конечных разностей первого порядка образуются конечные разности второго порядка:  yi =  yi+ -  yi (i = , , , ...) Продолжая этот процесс, можно по заданной таблично функции составить таблицу конечных разностей (Таблица 1). Таблица 1 x y  y  y  y ... x y  y  y  y   x y  y  y  y   x y  y  y ...   x y  y ...     x4 y4 ...       ... ...         Конечные разности любого порядка могут быть представлены через значения функции. Действительно, для разностей первого порядка это следует из определения. Для разностей второго порядка имеем:  yi =  yi +  -  yi = (yi+ - yi+) - ( yi+ - yi) = yi+ - yi+ + yi. Для разностей 3-го порядка  yi =  yi +  -  yi = ( yi +  -  yi + ) - ( yi +  -  yi) = ( yi +  - yi +  + yi + ) - - (yi +  - yi +  + yi ) =yi +  - yi + + yi +  - yi. Методом математической индукции можно доказать, что:  kyi = yi + k - k yi +k -  + y i + k -  - ... + (-) k yi. (12) Пример 4. Построить конечные разности для функции P (x) = х3, считая шаг  х = 1.  P (x) = (х + 1)3 - х3 = 3 х2 + 3 х + 1,  2P (x) = [3 (х + 1)2 + 3 (х + 1) + 1] -(3 х2 + 3 х + 1) = 6 х + 6,  3P (x) = [6 (х + 1) + 6] - (6 х + 6) = 6,  nP (x) = 0 при n > 3. Обратите внимание, что конечная разность третьего порядка функции P (x) постоянна. Вообще, справедливо утверждение: если Pn (x) = - полином n-й стaiени, то  nPn (x) = n! an hn = const, где  x = h. Пример 5. Составить таблицу конечных разностей функции у = 2 х3 - 2 х2 + 3 х - 1 от начального значения x0 = 0, приняв шаг h = 1. Полагая x0 x 1x2, находим соответствующие значения y-1  y 2y13. Отсюда имеем:  y = yy3,  y = yy11,  2y =  yy8. Эти значения заносим в таблицу (Таблица 2). Таблица 2 x y  y  y  y -1 3 8 12 1 2 11 20 12 2 13 31 32 12 3 44 63 44 12 4 107 107 56 12 5 214 163 68 12 . . . . . . . . . . . . . . . Так как наша функция есть полином третьей степени, то третья разность ее постоянна и равна 3yi = 2  3! = 12. Поэтому дальнейшее заполнение Таблицы 2 можно производить при помощи суммирования, используя формулы  2y=  2y+ 12 (i =0, 1, 2, …),  y=  y+  2y(i = 1, 2, …), y= y+  y(i = 2, 3, …). Первая интерполяционная формула Ньютона Пусть для функции, заданной таблично с постоянным шагом, составлена таблица конечных разностей. Будем искать интерполяционный полином в виде: Pn(x) = a + a (x - x) + a(x - x) (x - x) + + a(x - x) (x - x) (x - x) + ... + an (x - x) (x - x)  ... (x - xn-). (13) Это полином степени n. Значение коэффициентов a, a, ..., an Найдем из условия совпадения значений исходной функции и полинома Pn(x) a узлах интерполяции. Полагая х = х, из (13) находим: y = Pn (x) = a a= y. Далее, придавая х значение х и х последовательно, получаем: y = Pn(x) = a + a (x - x) = a + a h y = Pn (x) = a + a (x - x) + a(x- x) (x - x) = = a+ 2 ah + 2 ah Найдем коэффициенты а, а, а: a = y, а =! =! =!, a = = ! =!. Aналогично можно найти и другие коэффициенты. Общая формула имеет вид: ак =!, k = , , ,... n Подставляя эти выражения в формулу (13), получим следующий вид интерполяционного полинома: Pn (x) = у + !(x - x) + ! (x - x) (x - x) + ... … + !(x - x) (x - x)...(x - xn-).   (14) Практически формула (14) применяется в несколько ином виде. Положим t = , т. е. х = х + ht. Тогда !=! =! t - , =! = t -  и т. д. Окончательно имеем:   (15) Полученное выражение называется первой интерполяционной формулой Ньютона для интерполирования вперед. Интерполяционную формулу (15) обычно используют для вычисления значений функций в точках левой половины рассматриваемого отрезка. Это объясняется следующим. Разности  kуi вычисляются через значение функции уi, уi + , ..., уi + k, причем i + k  n, поэтому при больших значениях i мы не можем вычислить разности высших порядков (k   n - i). Например, при i = n -  в (15) можно учесть только  у,  у,  у. Если в формуле (15) положить n = 1 , то получим формулу линейной интерполяции: P1 (x) = уt  у. При n = 2 будем иметь формулу квадратичной интерполяции: P2 (x) = уt  у. Если дана неограниченная таблица значений у, то число n в интерполяционной формуле (15) может быть любым. Практически в этом случае число n выбирают так, чтобы разность  nуi была постоянной с заданной степенью точности. За начальное значение x0 можно принять любое табличное значение аргумента х. Если таблица значений функции конечна, то число n ограничено, а именно: n не может быть больше числа значений функции у, уменьшенного на единицу. Пример 6. Приняв шаг h = 0,05, построить на отрезке [3,5; 3,6] интерполяционный полином Ньютона для функции у = ех, заданной таблицей х 3,50 3,55 3,60 3,65 3,70 у 33,115 34,813 36,598 38,475 40,447 Составляем таблицу разностей (Таблица 3). x y  y  y  y 3,50 33,115 1698 87 5 3,55 34,813 1785 92 3 3,60 36,598 1877 95   3,65 38,475 1972     3,70 40,447 Заметим, что в столбцах разностей, следуя обычной практике, мы не указываем десятичные разряды (которые ясны из столбца значений функции). Так как разности третьего порядка практически постоянны, то в формуле (15) полагаем n = 3. Приняв x0 = 3,50, y0 = 33,115, будем иметь: P3 (x) = 33,115 + 1,698 t + 0,087+ 0,005 или P3 (x) = 33,115 + 1,698 t + 0,0435+ 0,00083, где . На практике часто необходимо для функции, заданной таблично, подобрать аналитическую формулу, представляющую с некоторой точностью данные табличные значения функции. Такая формула называется эмпирической, причем задача построения ее неоднозначна. При построении эмпирической формулы следует учитывать общие свойства функции. Если из таблицы разностей будет обнаружено, что n-е разности функции для равностоящих значений аргумента постоянны, то в качестве эмпирической формулы можно взять соответствующую первую интерполяционную формулу. Пример 7. Построить эмпирическую формулу для функции у, заданной таблично х 1 2 3 4 5 у 5,2 8,0 10,4 12,4 14,0 15,2   x y  y  y 5,2 2,8 -0,4 1 8,0 2,4 -0,4 2 10,4 2,0 -0,4 3 12,4 1,6 -0,4 4 14,0 1,2   5 15,2     Таблица 4 Составляя таблицу разностей (Таблица 4), убеждаемся, что вторая разность постоянна. Используя интерполяционную формулу Ньютона в форме (14) и учитывая, что h = 1, будем иметь: у = 5,2 + 2,8 х - или у = 5,2 + 3 х - 0,2 х2.   Вторая интерполяционная формула Ньютона Для правой половины рассматриваемого отрезка разности лучше вычислять справа налево. В этом случае: t= , o. е. t  0 и интерполяционную формулу Ньютона можно получить в виде: (16) Формулу (16) называют второй интерполяционной формулой Ньютона для интерполирования назад. Пример 8. Дана таблица значений y = lg x семизначных логарифмов х у 1000 3,0000000 1010 3,0043214 1020 3,0086002 1030 3,0128372 1040 3,0170333 1050 3,0211893 Найти lg 1044. Составляем таблицу разностей (Таблица 5). Таблица 5. Таблица разностей функции у = lg x х у  y  y  y 1000 3,0000000 43214 - 426 8 1010 3,0043214 42788 - 418 9 1020 3,0086002 42370 - 409 8 1030 3,0128372 41961 - 401   1040 3,0170333 41560     1050 3,0211893       Примем xn = 1050, тогда Используя подчеркнутые разности, в силу формулы (16) будем иметь: lg 1044 = 3,0211893 + (-0,6)  0,0041560 +  0,0000401 +  0,0000008 = 3,0187005. Как первая, так и вторая интерполяционные формулы Ньютона могут быть использованы для экстраполирования функции, т. е. для нахождения значений функции y для значений аргументов х, лежащих вне пределов таблицы. Если х < х0 и х близко к х0, то выгодно применять первую интерполяционную формулу Ньютона, причем тогда < 0. Если х < хn и х близко к хn, то удобнее пользоваться второй интерполяционной формулой Ньютона, причем > 0. Таким образом, первая интерполяционная формула Ньютона обычно используется для интерполирования вперед и экстраполирования назад, а вторая интерполяционная формула Ньютона, наоборот - для интерполирования назад и экстраполирования вперед. Операция экстраполирования менее точна, чем операция интерполяции в узком смысле слова.   Точность интерполяции   График интерполяционного полинома у =  (х) проходит через заданные точки, т. е., значения полинома  (х) и данной функции у =  (х) совпадают в узлах х = хi (i = , , .., . n). Если функция  (х) сама является полиномом степени n, то имеет место тождественное равенство  (х) =  (х). В общем случае в точках, отличных от узлов интерполяции R(x) =  (х) -  (х). Эта разность и есть погрешность интерполяции, и называется остаточным членом интерполяционной формулы. Оценим его значение. Предположим, что заданные числа yi являются значениями некоторой функции y =  (х) в точках х = хi . Пусть эта функция непрерывна и имеет непрерывные производные до n + 1 порядка включительно. Можно показать, что в этом случае остаточный член интерполяционного полинома Лагранжа имеет вид: RL (x) =  (n+)( ) где  (n+) ( ) - производная (n+)-го порядка функции  (х) в некоторой точке х =  ,   [x, xn]. Если максимальное значение этой производной равно: Mn +  =  (n+) (x) , то можно записать формулу для оценки остаточного члена интерполяционного полинома Лагранжа  RL (x) Mn+ (17) Остаточный член первой интерполяционной формулы Ньютона можно записать в виде: RP (x) =  (n+) ( ) h n + , . (18) Остаточный член второй интерполяционной формулы Ньютона: RP (x) =  (n+) ( ) h n + , . (19) Если предположить, что разности  n + уn почти постоянны для функции y = f(x) и h достаточно мало, и учитывая, что  (n+) (х) = , приближенно можно положить:  (n+) ( )  . В этом случае можно записать следующую формулу остаточного члена первой интерполяционной формулы Ньютона: RP (x)   n+у. Формула остаточного члена второй интерполяционной формулы Ньютона имеет вид: RP (x)   n+уn . Точность интерполяции   График интерполяционного полинома у =  (х) проходит через заданные точки, т. е., значения полинома  (х) и данной функции у =  (х) совпадают в узлах х = хi (i = , , .., . n). Если функция  (х) сама является полиномом степени n, то имеет место тождественное равенство  (х) =  (х). В общем случае в точках, отличных от узлов интерполяции R(x) =  (х) -  (х). Эта разность и есть погрешность интерполяции, и называется остаточным членом интерполяционной формулы. Оценим его значение. Предположим, что заданные числа yi являются значениями некоторой функции y =  (х) в точках х = хi . Пусть эта функция непрерывна и имеет непрерывные производные до n + 1 порядка включительно. Можно показать, что в этом случае остаточный член интерполяционного полинома Лагранжа имеет вид: RL (x) =  (n+)( ) где  (n+) ( ) - производная (n+)-го порядка функции  (х) в некоторой точке х =  ,   [x, xn]. Если максимальное значение этой производной равно: Mn +  =  (n+) (x) , то можно записать формулу для оценки остаточного члена интерполяционного полинома Лагранжа  RL (x) Mn+ (17) Остаточный член первой интерполяционной формулы Ньютона можно записать в виде: RP (x) =  (n+) ( ) h n + , . (18) Остаточный член второй интерполяционной формулы Ньютона: RP (x) =  (n+) ( ) h n + , . (19) Если предположить, что разности  n + уn почти постоянны для функции y = f(x) и h достаточно мало, и учитывая, что  (n+) (х) = , приближенно можно положить:  (n+) ( )  . В этом случае можно записать следующую формулу остаточного члена первой интерполяционной формулы Ньютона: RP (x)   n+у. Формула остаточного члена второй интерполяционной формулы Ньютона имеет вид: RP (x)   n+уn .   Пример 9. Оценим ошибку интерполяции при отыскании lg 1044 в Примере 8. Составленная таблица разностей (Таблица 5) показывает, что третьи разности функции можно считать практически постоянными, следовательно n = 3. По формуле (19) имеем: R3(x) =  IV ( ) h 4, Так как f(x) = lg10 х, то  IV (х) = - , поэтому | IV ( ) | < . Далее h = 10 и t = - 0,6, так что окончательно R3(x)  10-10. Таким образом, интерполирование с учетом третьих разностей, т. е. с применением интерполяционного полинома третьей степени, не вносит погрешностей до десятого знака. Следует подчеркнуть, что существует один и только один интерполяционный полином при заданном наборе узлов интерполяции. Формулы Лагранжа, Ньютона и др. порождают один и тот же полином (при условии, что вычисления проводятся точно). Повышение точности интерполяции целесообразно проводить за счет уменьшения шага и специального расположения точек хi. Выбор способа интерполяции определяется различными соображениями: точностью, временем вычисления, погрешностями округления и др. Интерполяция сплайнами. Пусть отрезок [a,b] разбит точками на n частичных отрезков . Сплайном степени m называется функция , обладающая следующими свойствами: 1) функция непрерывна на отрезке [a,b] вместе со своими производными до некоторого порядка  p. 2) на каждом частичном отрезке функция совпадает с некоторым алгебраическим многочленом степени m. Разность m-p между степенью сплайна и наивысшим порядком непрерывной на отрезке [a,b] производной называют дефектом сплайна. Кусочно-линейная функция является сплайном первой степени с дефектом, равным единице. Действительно, на отрезке [a,b] сама функция (нулевая производная) непрерывна. В то же время на каждом частичном отрезке совпадает с некоторым многочленом первой степени. КОНТРОЛЬНЫЕ ВОПРОСЫ: 1. Что называют дефектом сплайна. 2.Какими свойствами обладает функция сплайном степени 3. Опишите процесс изменение частоты дискретизации сигнала ЛЕКЦИЯ 12. СПОСОБЫ ИНТЕРПОЛЯЦИЯ. МНОЖЕСТВЕННАЯ ФИЛЬТРАЦИЯ. интерполяция метод графический Интерполяция — способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений. Задачи интерполяции: научиться вычислять значение функции в любой наперед заданной точке. Интерполяция иногда делится на два вида: 1) - собственная интерполяция. 2) - экстраполяция. В данной работе я бы хотел продемонстрировать принципы работы четырех методов интерполирования: Метод Лагранжа, Метод Эйткена, Метод Ньютона и Метод кубических сплайнов. Методы интерполяции Интерполяционный многочлен Лагранжа — многочлен минимальной степени, принимающий данные значения в данном наборе точек. Для  пар чисел , где все  различны, существует единственный многочлен  степени не более , для которого . В простейшем случае () — это линейный многочлен, график которого — прямая, проходящая через две заданные точки. Лагранж предложил способ вычисления таких многочленов: где базисные полиномы определяются по формуле:  обладают следующими свойствами: являются многочленами степени   при  Отсюда следует, что , как линейная комбинация , может иметь степень не больше , и  Код метода: double p(double X,int n) // В качестве входных параметров используется { // double X – переменная относительно которой double result=0,temp=1; // Нужно найти f(x), а также int n – количество for (int i=0;i<=n;i++) // точек интерполяции. Матрицы x[], y[] заданы { // глобально. for (int j=0;j<=n;j++) if (i!=j) temp*=(X-x[j])/(x[i]-x[j]); result+=temp*y[i]; temp=1; } return result; } Метод Эйткена Итерационно-интерполяционный метод Эйткена позволяет свести вычисления коэффициентов интерполяционного полинома Лагранжа, с учетом его равенства в узлах интерполяции с исходными данными к вычислению функциональных определителей второго порядка. При этом эффективность метода повышается в тех случаях, когда нет необходимости в получении приближенного аналитического выражения функции f(х), заданной таблично, а требуется лишь определить значение в некоторой точке х*, отличной от узловых точек. Этот метод заключается в последовательной линейной интерполяции. Процесс вычисления f(x*) состоит в следующем: необходимо пронумеровать узлы интерполяции, например, в порядке убывания их от х*. Затем для каждой узловой точки интерполяции строятся соотношения, которые является интерполяционными полиномами, построенными соответственно по узлам хi, хj, хk. Полученный полином является интерполяционным полиномом, построенный по узлам хi, xj, …, хk, хm. Это утверждение верное, так как Рn-1ij…k(х) и Рn-1j…km(x) являются интерполяционными полиномами. При его реализации предполагается, что функция гладкая, а также критерием оценки погрешности определяется некоторое значение, определяемое условиями конкретной задачи. Код метода: Double P(double X,int n) // В качестве входных параметров используется { // double X – переменная относительно которой float l[n+1],result; // Нужно найти f(x), а также int n – количество for (int i=0;i<=n;i++) // точек интерполяции. Матрицы x[], y[] заданы l[i]=y[i]; // глобально. for (int i=1;i<=n;i++) for (int j=1;j+i<=n;j++) { l[j]=(l[j]*(x[i+j]-X)-l[j+1]*(x[j]-X))/(x[i+j]-x[j]); result=l[j]; } return result; Метод Ньютона Формула Ньютона: где Конечной разностью функции у=f(х) называется функция , где h – фиксированный шаг. Конечные разности иногда называются конечными разностями первого порядка. Функция обозначается: Принимаем В программе для подсчета конечных разностей мы используем следующую функцию: double Delta(double X, int n, double h) //Где n – порядок, h - шаг { if (n>1) return (Delta(X+h,n-1,h)-Delta(X,n-1,h)); else return (f(X+h)-f(X)); } Код программы: double Pn(double X, int h, int n) // В качестве входных параметров используется { // double X – переменная относительно которой double result=y[0],qq,q; // Нужно найти f(x), а также int n – количество q=(X-x[0])/h; // точек интерполяции. Матрицы x[], y[] заданы for (int i=1;i<=n;i++) // глобально. { qq=1; for (int j=1;j<=i;j++) qq*=(q-j+1); result+=Delta(x[0],i,h)*qq/fact(i); } return result+1; } { float l[n+1],result; for (int i=0;i<=n;i++) l[i]=y[i]; 2.4 Метод кубических сплайнов Кубическим сплайном на сетке x0,x1,…xn называется функция S(х), которая обладает следующими свойствами: на каждом интервале [хi-1, хi], где 1  i  n, функция S(х) является кубическим многочленом (на каждом интервале свой многочлен). на всем интервале [х0, хп] S(х) – дважды непрерывно дифференцируемая функция на краях интервала вторая производная обращается в ноль (краевое условие). S΄΄(x0)=S΄΄(xn)=0 3’. для периодических кубических сплайнов. S΄΄(x0)=S΄΄(xn)=0 ; S΄(x0)=S΄(xn)=0 CM=d Элементы матрицы С вычисляются по формуле: , () - вектор правых частей. Код программы: const int n = 5; double x[n]; double y[n]; float C[n-1][n-1]; float Cd[n-1][n]; double S(double X, int n) { int i,j; float S,k1,k2,k3,k4; float h[n],d[n-1],M[n+1],G[n-1]; for (i=0;i 1) C[i][j] = 0; } for (i=0;i=0;i--) { float summa=0; for (j=N-1;j>i;j--) summa+=y[j]*Cd[i][j]; y[i]=(Cd[i][N]-summa)/Cd[i][i]; } } 3. Графическое представление Для графического представления работы методов написана функция graphic (int n, int method) void graphic(int n, int method) { float X, Y; initwindow(1000,640); line(500,0,500,640);//y line(490,20,500,0); line(510,20,500,0); outtextxy(520,10,"y"); line(0,320,1000,320);//x line(990,310,1000,320); line(990,330,1000,320); outtextxy(985,330,"x"); X=-100; moveto(500+X*20,160+(160-(Y*3))); setcolor(4); for (int i=0;i<=n;i++) //Отмечаем на графике заданные точки интерполяции circle(500+x[i]*20,160+(160-(y[i])),5); setcolor(10); do { if (method==0) //Выбираем метод { Y=p(X,4); outtextxy(520,850,"Lagranj"); } if (method==1) { Y=P(X,4); outtextxy(520,850,"Eytkin"); } if (method==2) { Y=Pn(X,1,4); outtextxy(520,850,"Newton"); } if (method==3) { Y=S(X,n); // С помощью определенного метода находим y для outtextxy(520,850,"Splayn"); // Заданного x } lineto(500+X*20,160+(160-(Y))); //И проводим к точке (x,y) линию от предыдущей точки X=X+0.02; } while(X<=100); getch(); closegraph(); } Таким образом, например, для метода Лагранжа и значений точек интерполяции: x[0]=1; x[1]=2; x[2]=3; x[3]=4; x[4]=5; y[0]=1; y[1]=4; y[2]=9; y[3]=16; y[4]=25; Получаем график функции: КОНТРОЛЬНЫЕ ВОПРОСЫ 1. Перечислите методы интерполяции 2. Напишите формулу Лагранжа 3. Напишите формула Эйткена 4. Напишите формула Ньютона 5. Охарактеризуйте метод кубических сплайнов