Виды цифровых сигналов, способы их получения. Математическое описание, операции с ЦС

Федеральное государственное бюджетное образовательное учреждение высшего образования «Национальный исследовательский университет «МЭИ» Курс лекций дисциплины вариативной части математического и естественнонаучного цикла Б.1 «Цифровая обработка сигналов» Направление подготовки 09.03.01 Информатика и вычислительная техника Профили: Автоматизированные системы обработки информации и управления, Вычислительные машины, комплексы, системы и сети, Системы автоматизированного проектирования Вычислительно-измерительные системы Авторский коллектив: Доцент кафедры ЭФИС, к.т.н. «01» февраля 2016 г. _________________ С.Н. Михалин ЗАДАЧИ КУРСА Цифровая обработка сигналов (ЦОС) – раздел науки, посвященный теоретическому обоснованию методов анализа (синтеза) сигналов и алгоритмов их обработки. Сигнал – есть функция, связывающая и характеризующая какие-либо физические параметры. Сигнал, как функцию, отождествляют с каким-либо физическим процессом и, соответственно, его параметрами (например: скорость, ускорение, температура, напряжение, сопротивление, масса, влажность, давление и т.п.). Таким образом, сигнал содержит некоторую информацию о физическом процессе, который он характеризует. Кроме полезной информации в сигнале присутствуют ненужные составляющие – помехи, происхождение которых имеет разную природу. Именно задачи выделения, преобразования, передачи полезной информации из сигнала (или их группы) являются предметом ЦОС. Обычно реальные сигналы представляют собой функцию напряжения от времени: u (t ) и содержат помехи и шум (помехи можно рассматривать как ненужные сигналы, шумы – это свойство реальных объектов). Типичные задачи ЦОС: 1) обнаружение и анализ сигнала при воздействии помех (шумов) – радиолокация, дефектоскопия, прием данных (информации) и т.п.; 2) помехоустойчивость, сжатие и кодирование сигнала – передача данных и безопасность; 3) реализация интерфейсов «человек-машина» – распознавание объектов (видеоизображений), анализ и синтез речевых сигналов; 4) распределенные системы управления и сбора информации. На сегодняшний день маловероятно найти область науки, где ЦОС еще не нашла своего применения. Поэтому ЦОС, как область знаний, никогда не потеряет своей актуальности и значимости. 2 Задачами данного курса являются: − введение терминологии свойственной предмету цифровой обработки сигналов; − ознакомление с видами сигналов, способами их получения, математического описания, операциями с ними; − изучение основных преобразований и операций, широко используемых на практике, рассмотрение их свойств; − ознакомление с основными методами обработки цифровых сигналов; − ознакомление с самыми распространенными алгоритмами цифровой обработки сигналов; − ознакомление с основными проблемами технической реализации алгоритмов на практике; − ознакомление с программной средой Matlab как средством моделирования систем обработки сигналов. 3 Раздел 1. Элементы теории сигналов Реальные физические сигналы почти всегда описываются непрерывной функцией (или кусочно-непрерывной функцией). Такие сигналы называются аналоговыми. Если функция существует только в определенные моменты времени (ее аргумента), то множество значений такой функции называется дискретным сигналом. Цифровым сигналом будем называть дискретный сигнал, который может принимать значения только из конечного множества чисел. Конечной длительности (финитные) Бесконечной длительности Детерминированные – полностью известная функция (значение функции известно для любого аргумента) Случайные Сигналы Периодические (при любом t выполняется s(t)=s(t+nT), T – период, n – целое число) Сигналы с интегрируемым квадратом: ∫s2(t)dt<∞ – значение функции для каждого аргумента есть случайная величина, устанавливающаяся с некоторой вероятностью Стационарные Непериодические Нестационарные Рис. 1.1 – Классификация сигналов Дополнительно сигналы можно подразделить на одномерные (функции одной переменной) и многомерные (функции многих переменных). Примеры детерминированных одномерных сигналов представлены на рис. 1.2. период Детерминированный, финитный, непериодический Детерминированный, бесконечный, непериодический Детерминированный, бесконечный, периодический Рис. 1.2 – Примеры детерминированных сигналов 4 Детерминированный сигнал это детерминированная функция, значения которой в любой точке области определения можно рассчитать заранее, зная математическую модель сигнала (имеется в виду формула, таблица, закон и т.п.). В частности, для переменного напряжения u (t ) = A sin (t ) можно получить значение сигнала в любой момент времени t (здесь модель определяется словами переменное синусоидальное напряжение с частотой ω и амплитудой A). Надо отметить, что детерминированные сигналы в «чистом виде» не могут существовать в природе, т.к. на реальные сигналы всегда оказывает влияние окружающая среда, поэтому реальный сигнал всегда содержит случайную составляющую. Однако понятие «детерминированный сигнал» облегчает математическое описание сигналов. На практике же реальный сигнал представляют как сумму «идеального» детерминированного сигнала и случайного сигнала. Случайный сигнал это функция, значение которой невозможно предсказать (вычислить заранее), т.е. значение функции для каждого аргумента есть случайная величина, устанавливающаяся с некоторой вероятностью. Случайность сигнала, или говорят процесса, проявляется в том, что значение функции изменяется от одного опыта к другому. Полученная в результате каждого опыта функция называется реализацией случайного процесса, т.е. это то, что мы регистрируем в ходе наблюдения (измерения) случайного сигнала. На рис. 1.3 представлен случайный сигнал. Рис. 1.3 – Реализации случайного финитного процесса 5 В каждый момент времени существует бесконечно много реализаций случайного процесса. Осуществляя измерения такого сигнала, получают конкретную, одну из многих, реализацию. Примером периодического случайного сигнала может служить сигнал, состоящий из периодически повторяющихся одинаковых фрагментов, значения которых на каждом периоде заранее не известны, но каждая реализация содержит повторяющийся сигнал. Например: электрический сигнал (с камеры наблюдения) развертки строк экрана – получая сигнал, мы не знаем наперед о цвете каждого пикселя изображения – сигнал случайный. Однако, если изображение статично, то пиксели не изменяются и получаемые строки повторяются с частотой кадров – сигнал периодичен. Предсказать какая будет «картинка» в целом (т.е. кадр) в следующий момент времени невозможно – значит рассматриваемый сигнал случайный и периодичный. Источниками сигналов являются: − генераторы сигналов (обычно зависимость напряжения от времени) это устройства, предназначенные для генерации сигналов с заданными параметрами; − датчики (преобразователи физических величин в зависимость, как правило, напряжения от времени) это устройства, дающие возможность контролировать неэлектрические процессы (температура, давление и т.п.); − искусственно созданные сигналы: табличная форма, генераторы случайных чисел, формула и т.п. (обычно применяются при математическом моделировании). Например: зависимость окружающей температуры во времени (за час, день, неделю, год и т.п.), зависимость скорости автомобиля в единицу времени, зависимость яркости свечения лампы во времени, зависимость массы объекта во времени, зависимость давления во времени и т.п. – все эти физические параметры приводятся к зависимости напряжения от времени с помощью соответствующих датчиков (сенсоров): термопара, пьезоэлемент, фотоэлемент, катушка, открытый конденсатор и т.п. 6 Параметры сигналов. Сигналы описываются математическими моделями, которые имеют ряд параметров. При этом сигналы можно представлять во временной области – в виде функций от времени или в частотной области – в виде функций от частоты (также можно представлять и в других координатах, если это будет иметь смысл). • Амплитуда – максимальное значение (по модулю) смещения или изменения переменной величины от среднего значения при колебательном или волновом движении (изменении функции). • Период повторения (обратно пропорционален частоте) сигнала – расстояние между одинаково колеблющимися точками колебательного или волнового изменения сигнала. • Максимальное (минимальное) значение – величина отклонения или изменения функции f (t ) от нулевого значения, такая, что для любой величины t из области существования функции выполняется: f (t ) ≤ Fmax ( f (t ) ≥ Fmin ). • Среднее значение – величина, определяемая в соответствии с формулой: 1T M = ∫ x(t )dt , где T – интервал наблюдения функции или ее период. T0 Среднее значение для периодических сигналов называют постоянной составляющей. • Среднеквадратическое значение – величина, определяемая в соответствии с 1T 2 выражением: rms = ∫ x (t )dt , где T – интервал наблюдения функции или T0 ее период (совпадает по смыслу с действующим значением напряжения, если x(t ) – переменное напряжение на участке электрической цепи). • Полоса сигнала – разность верхней и нижней границ частотной области (спектра), в пределах которых существует (представляется) сигнал. 7 • Эффективная полоса сигнала (эффективная ширина спектра) F=  эф – полоса частот, в пределах которой заключается основная доля энергии  эф сигнала, обычно: 90%, т.е. ∞ 2 ∫ S ()d = 0.9 ∫ S ()d , S(ω) – спектральная 2 плотность сигнала (частотное представление). S(ω) ω ωэф 90% энергии сигнала • Длительность сигнала – интервал времени, в пределах которого локализован (представлен) сигнал. • Эффективная длительность сигнала T=tэф – отрезок времени, в пределах которого заключена основная доля энергии импульса, обычно: 90%. s(t) t tэф 90% энергии сигнала • База сигнала B =  эф t эф (например: наименьшее значение базы имеет гауссов импульс B=1.353, у прямоугольного импульса B=4.59). • Отношение сигнал шум (ОСШ) равно отношению мощности полезного сигнала к мощности шума (часто приводится в логарифмических единицах). • Динамический диапазон есть отношение P мощности сигнала к наименьшей: D = 10 lg max  Pmin наибольшей мгновенной   .  8 u (t ) = 1.5 + A sin(t ) + (t ) (t ) Umax A Uд (rms) Uср -A 2  2.96  D = 10 lg  ≈ 15.4 дБ  0 .5   2.75  ОСШ = 10 lg  ≈ 23.6 дБ  0.0121  Umin ωt Рис. 1.4 – Графическое представление основных параметров сигнала (A=1 В)  1  (*) 2.75 = 1.5 +    2 2 2 9 1.1 Детерминированные сигналы 1. Дельта функция (функция Дирака) ∞, t = 0 , при этом  (t ) =  0, t ≠ 0 δ(t) ∞ ∫ (t )dt = 1 t −∞ ∞ Фильтрующее свойство дельта функции: ∫ f (t )(t − t0 )dt = f (t0 ) −∞ 2. Функция единичного скачка (функция Хевисайда) σ(t) 1, t ≥ 0  (t ) =  0, t < 0 t 3. Гармонический сигнал (Um – амплитуда, ω – круговая циклическая частота, Ψ – начальная фаза) u (t ) = U m sin (t +  ) 4. Экспоненциальный сигнал (Um – амплитуда, α – коэффициент затухания) u (t ) = U m e − t Гауссов импульс: u (t ) = U m e −t Остальные сигналы 2 / 2 . можно представить комбинацией основных сигналов. Например: 1) c помощью суперпозиции функций Хевисайда, например – прямоугольный сигнал длительностью tи: x(t ) =  (t ) −  (t − t и ) . 2) с помощью разложения сигналов в тригонометрические ряды, x(t) например: представление меандра x(t) A с амплитудой A в виде бесконечного t тригонометрического ряда x(t ) = 4A  1 1 1  sin (t ) + sin (3t ) + sin (5t ) +  + sin (kt ) +     3 5 k  3) с помощью алгебраических операций над сигналами: положим x1(t), x2(t), y(t) – непрерывные вещественные сигналы 3.1) сложение (вычитание): y(t)=x1(t)±x2(t) x1(t) y(t) + x2(t) 10 3.2) умножение на вещественную константу (c): x1(t) 3.3) Задержка сигнала x1(t) на время τ: y(t) c y(t)=сx1(t) Задержка x1(t) y(t)= x1(t-τ) y(t) 3.4) умножение сигналов: y(t)= x1(t) x2(t) (если частоты сигналов существенно различны – то имеет место амплитудная модуляция) y(t) × x1(t) x2(t) +∞ 3.5) свертка сигналов: y (t ) = x1 (t ) * x2 (t ) = ∫ x1 ( ) x2 (t −  )d , операция линейна −∞ ∞ 3.6) скалярное произведение: s(t ), v(t ) = ∫ s(t )v(t )dt −∞ u(t), В Произведение: u1 (t ) u 2 (t ) 2 2. 0V 1 1. 0V u(t), В 2 2. 0V 0V u1 -1 -1. 0V 1 1. 0V -2 -2. 0V 0s 0. 2ms V( MULT1: OUT) 0.2 u(t), В 0V u2 -1 -1. 0V 0s V( V2: +) 0. 2ms V( V1: +) 0.2 0. 4ms 0.4 0. 6ms Ti me 0.6 0. 8ms 0.8 1. 0ms t, мс 3 3. 0V 0. 4ms 0.4 0. 6ms Ti me 0.6 0. 8ms 0.8 1. 0ms t, мс Сумма: u1 (t ) + u 2 (t ) Задержка сигнала u1 на 0.1 мс 2 2. 0V 1 1. 0V 0V -1 -1. 0V 0s u(t), В V( SUM1: OUT) 0. 2ms 0. 4ms 0.2 0.4 0. 6ms Ti me 0.6 0. 8ms 0.8 1. 0ms t, мс 2 1 0.2 0.4 0.6 0.8 t, мс Рис. 1.5 – Примеры операций над сигналами 11 1.2 Пространство сигналов Для удобства анализа и обработки информации, которая заключена в сигналах, их множество «помещают» в подходящее метрическое пространство (как правило, линейное, с заранее оговоренными свойствами и единицами измерений). Это позволяет выделять из множества сигналов сигналы с определенными параметрами, сравнивать сигналы друг с другом, оценивать их изменение при их прохождении через системы обработки данных и т.п. Линейное пространство аналоговых сигналов с введенным скалярным произведением (положительно определенным) называется Гильбертовым пространством Н (второе распространенное обозначение – L2), представляет собой обобщение Евклидова пространства на бесконечномерный случай. Линейное пространство дискретных и цифровых сигналов называется пространством Эвклида (обозначение пространства – R2). В пространствах Гильберта и Эвклида определяют норму и метрику в соответствии с выражениями: s(t ) = s, s = ∞ 2 ∫ s (t ) и (s, v ) = s(t ) − v(t ) соответственно. −∞ В этих пространствах справедливо неравенство Коши-Буняковского: s(t ), v(t ) ≤ s(t ) ⋅ v(t ) – скалярное произведение векторов не превосходит произведения их норм. 1.2.1 Разложение сигнала по базисным функциям В линейном пространстве L (размерности N) всегда можно выделить множество векторов {xn ; n = 0, 1, 2,, N − 1}, для которых выполняется N −1 равенство нулю их линейной комбинации: ∑  n xn = 0 только при условии n =0 равенства нулю всех значений  n . Такое множество векторов называется линейно независимым. Ни один вектор линейно независимого множества не может быть выражен в виде какой-либо линейной комбинации других векторов этого пространства. Такое множество векторов называется базисом N-мерного пространства L. Линейная комбинация таких N линейно 12 независимых векторов образует векторное пространство, где каждый вектор U может быть выражен единственной линейной комбинацией векторов xn: U= N −1 ∑  n xn . Совокупность чисел { n } называется спектром вектора U в этом n =0 базисе. Спектр вектора в общем случае может быть комплексным. Произвольный сигнал s(t ) ∈ H (пространство Гильберта), заданный на интервале [a, b] , может быть разложен в ряд по упорядоченной системе ∞ ортонормированных базисных функций u n (t ) : s(t ) = ∑ cn u n (t ) . n =0 Для нахождения значений коэффициентов сn умножим обе части данного выражения на базисную функцию u m (t ) с произвольным номером m и проинтегрируем результаты по переменной t, при этом получим: ∞ b b ∫ s(t )u m (t )dt = k∑=0 ck ∫ u m u k dt . a a С учетом ортонормированности функций u i (t ) , в правой части этого равенства b остается только один член суммы с номером m=k при ∫ u k u k dt = 1 , который, по a левой части уравнения, представляет собой скалярное произведение сигнала и соответствующего m=k базисного вектора, т.е. проекцию сигнала на b соответствующее базисное направление: ck = ∫ s(t )u k (t )dt . a Таким образом, в геометрической интерпретации коэффициенты сk представляют собой проекции вектора (сигнала s(t ) ) на соответствующие базисные направления s(t ) u k (t ) , т.е. координаты вектора образованному системой ортогональных функций u (t ) , в базисе, в пределе – бесконечномерной. Возможность обобщенные ряды разложения по непрерывных системам сигналов ортогональных и функций функций в имеет принципиальное значение, так как позволяет вместо изучения несчетного множества точек сигнала ограничиться счетной системой коэффициентов 13 ряда. К системам базисных функций, которые используются при разложении сигналов, предъявляют следующие основные требования: - для любого сигнала ряд разложения должен сходиться; - при ограничении ряда по уровню остаточной погрешности расхождения с заданным сигналом количество членов ряда должно быть минимальным; - базисные функции должны иметь достаточно простую аналитическую форму и коэффициенты разложения в ряд должны вычисляться относительно просто. Согласно теореме Дирихле, любой сигнал s(t ) , имеющий конечное число точек нарушения непрерывности первого рода, и конечный по энергии на интервале [a, b] , может быть разложен по системе ортонормированных функций, если существуют интегралы модуля сигнала и модуля его первой b производной: ∫ s(t ) dt < ∞ , a b ∫ s′(t ) dt < ∞ . a Это означает, что любой физически полученный сигнал представим во временной области (в виде функции) и имеет спектр (в частности, при представлении сигнала по базису гармонических функций – сигнал представим также в частотной области). Именно представление сигналов в пространстве позволяет их рассматривать как множество векторов и применять методы математического анализа, статистики. 1.2.2 Примеры ортогональных систем в L2 Широкое применение в цифровой обработке сигналов нашли некоторые функциональные базисы ортогональных функций. Рассмотрим часто встречающиеся базисы. ∞   2kt   2kt  1. Тригонометрическая система функций 1, cos , sin    T   T  k =1  является полной в пространстве L2 [a, a + T ] на любом отрезке t ∈ [a, a + T ]. 14 2. Система функций Уолша {wn ( x )}n=0 является полной в пространстве L2 [0, 1] ∞ и ортонормирована на полуинтервале [0, 1) . Представим целое число n≥0 в K виде двоичного разложения: n = ∑ nk 2 k , где nk ={0, 1}. k =0 Тогда функции системы Уолша выражаются следующим образом: wn ( x ) = ∏ [rk ( x )] k , где K n k =0 число K – определяется номером функции Уолша n<2K; rk(x) – функция Радамахера. Функции Уолша n n2 1 2 3 4 5 1 1 6 7 Таблица 1. Значения функций Уолша wn ( x ) n1 n0 w0 ( x ) = 1 w1 ( x ) = r0 ( x ) 1 1 1 1 1 1 1 1 1 1 принимают w2 ( x ) = r1 ( x ) w3 ( x ) = r1 ( x ) r0 ( x ) w4 ( x ) = r2 ( x ) w5 ( x ) = r2 ( x ) r0 ( x ) w6 ( x ) = r2 ( x ) r1 ( x ) w7 ( x ) = r2 ( x ) r1 ( x ) r0 ( x ) значения ±1, что удобно при программировании. w3(x) w5(x) 1 w7(x) 1 0.5 1 x -1 1 1 -1 x 1 -1 0.125 0.5 x 0.125 0.5 Рис. 1.6 – Функции Уолша 3. Система функций Хаара {hn (x )}∞n=0 на полуинтервале ортонормированной, полной в пространстве [0, 1) является L2 [0, 1] и определяется следующим образом. Положим h0 ( x ) = 1, для n>0 представим номер базисной функции в виде: n = 2 k + m , где целые числа k≥0 и 0 ≤ m ≤ 2 k − 1 . Тогда 15  0.5 k  2m 2m + 1  2 при x ∈   2 k +1 , 2 k +1     2m + 1 2m + 2  hn ( x ) = − 2 0.5 k при x ∈  k +1 , k +1  2  2   0 в других случаях   n=1; k=0; m=0 h1(x) n=2; k=1; m=0 h2(x) 20.5 1 0.5 20.5 x 1 2 x 1 -1 n=4; k=2; m=0 h5(x) 0.5 n=5; k=2; m=1 2 0.5 x 1 0.5 n=6; k=2; m=2 h6(x) 2 0.25 x 1 0.5 h4(x) n=3; k=1; m=1 h3(x) 0.25 x 1 0.5 x 1 Рис. 1.7 – Функции Хаара (*) Функция Радамахера определяются для x ∈ [0, 1) : 1 при x ∈ [0; 0.5) , r0 ( x ) =  ) [ − 1 при x ∈ . 5 ; 1  причем r0(x) периодически продолжается на всю числовую ось. Остальные функции Радамахера определяются выражением: rk ( x ) = r0 (2 k x ), k=1, 2, … r0(x) r1(x) 1 r2(x) 1 0.5 1 -1 1 x 0.5 1 x -1 1 -1 x 0.125 0.5 Рис. 1.8 – Функции Радамахера (k=0, 1, 2) (*) Ортогональная система функций { k }∞k =0 ⊂ H является полной в гильбертовом пространстве тогда и только тогда, когда ∀x ∈ H выполняется: ∞ x = ∑ 2k  k 2 2 (равенство Парсеваля-Стеклова). Иначе говоря, существует k =0 ∞ разложение x = ∑  k  k – ряд Фурье, числа k – коэффициенты ряда Фурье. k =0 16 1.3 Корреляция сигналов Автокорреляционная функция сигнала x(t ) есть скалярное произведение +∞ сигнала и сдвинутой во времени его копии: R( ) = ∫ x(t ) x(t −  )dt . −∞ +∞ Взаимная корреляционная функция: B( ) = ∫ x1 (t ) x2 (t −  )dt . −∞ Корреляционная функция применяется в случаях, когда необходимо обнаружить сигнал в смеси нескольких сигналов разной природы, т.к. корреляционная функция характеризует степень сходства сигналов (чем два вектора ближе по длине и направлению, тем больше величина их скалярного произведения). y ā1 y1 y2 β x1 ā2 x2 a1 , a 2 = a1 a 2 = a1 a 2 cos ( ) = x1 x2 + y1 y2 x Рис. 1.9 – Геометрическая интерпретация скалярного произведения Свойства автокорреляционной функции: +∞ 1) значение функции при τ=0 есть энергия сигнала: R(0 ) = ∫ x 2 (t )dt −∞ и при τ≠0 значение R( ) ≤ R(0 ) 2) функция четна: R( ) = R(−  ) 3) функция от сигналов с конечной энергией затухает 4) если сигнал не содержит особенностей в виде дельта функции, то его корреляционная функция непрерывна 5) если функция периодична, то ее автокорреляционная функция имеет тот же период. Примеры расчета корреляционных функций. 1) Автокорреляционная функция прямоугольного сигнала с периодом T=1 мс и длительностью импульса Tи≤T 17 x(t) A t Tи T R(τ)  A 2 (Tи −  ), 0 ≤  ≤ Tи  R( ) = ∫ x(t )x(t −  )dt =  A 2 (Tи +  ), − Tи ≤  ≤ 0 −T 0,  ≥ T и  A2Tи T Tи=0.5 мс -Tи Tи τ Tи=0.4 мс Рис. 1.10 – Автокорреляционная функция прямоугольного сигнала с длительностью импульса Tи 2) Автокорреляция гармонического сигнала x(t ) = A cos( 0 t +  ) с периодом 0.5T A2 2 2 cos( 0  ) : R( ) = ∫ A cos( 0 t +  )cos( 0 (t −  ) +  )dt = T= 2 0 −0.5T Автокорреляционная функция гармонического сигнала не зависит от начальной фазы и также является гармонической функцией с тем же периодом. 3) Автокорреляционная функция экспоненциального импульса с постоянной затухания τ: x(t ) = Ae − t /  , t≥0. Автокорреляционная функция экспоненциального импульса является экспоненциальной функцией. ∞ − t  R (t1 ) = ∫ Ae Ae − t + t1  dt = A e 2 − t1 ∞  ∫e − 2t  A 2  − 1 dt = e 2 t 18 R R τ=1 мс τ=0.1 мс t1 t1 Рис. 1.11 – Автокорреляционная функция экспоненциального сигнала с постоянной времени τ 4) Корреляционная функция зашумленного сигнала с ожидаемым, имеющим форму прямоугольного импульса Рис. 1.12 – Эталонный сигнал и его автокорреляционная функция (слева), зашумленный и сдвинутый эталонный сигнал и его корреляционная функция (справа) 19 1.4 Преобразование Фурье Преобразование Фурье сигнала s(t) есть: +∞ S ( ) = ∫ s(t ) e − jt dt – фактически представляет собой разложение сигнала по −∞ системе базисных тригонометрических функций (т.е. скалярное произведение сигнала на базисные функции). Комплексную величину S(ω) называют комплексным спектром (спектральной плотностью) сигнала s(t), ω – круговая частота. Обратное преобразование Фурье: s(t ) = 1 +∞ S ( ) e jt d . ∫ 2 −∞ (*) Напомним, что для представления функции в виде интеграла Фурье ∞ необходимо, чтобы функция была абсолютно интегрируемой: ∫ s(t ) dt < ∞ и −∞ являлась кусочно-гладкой на любом конечном отрезке [a, b] ⊂ (− ∞, ∞ ) . Трудно переоценить значимость тогда еще гипотезы французского математика Фурье (1768-1830 г.) о возможности разложения любой функции в набор синусоид с кратными частотами и наоборот – синтезировать любую функцию с помощью взвешенной суммы синусоид. Поскольку сейчас это является, несомненно, фундаментальной основой многих теоретических достижений, в том числе и основ ЦОС. Важнейшим случаем преобразования Фурье является представление периодических ∞ ∞ n =1 n =1 сигналов с периодом повторения T: s(t ) = A0 + ∑ An cos(n 0 t ) + ∑ Bn sin (n 0 t ) , где 0 = 2 – круговая циклическая частота сигнала (основная гармоника); T n 0 – круговая циклическая частота n-й гармоники; 1 0.5T A0 = ∫ s(t )dt – постоянная составляющая сигнала; T −0.5T 2 0.5T 2 0.5T ( ) ( ) An = ∫ s t cos n 0t dt , Bn = T ∫ s(t )sin (n 0t )dt – коэффициенты разложения. T −0.5T −0.5T 20 К примеру: если сигнал симметричный, то все Bn=0, если антисимметричный – все An=0. Для абсолютно интегрируемой функции (сигнала) коэффициенты An и Bn стремятся к нулю при n→∞. Записанные выражения есть ряд Фурье, который можно записать через комплексные величины, ряд Фурье для действительной функции f(t) есть: f (t ) = ∑ C k e − jk t = C 0 + ∑ [C k e − jk t + C −k e jk t ] = a0 + ∑ ak cos(k1t +  k ) , +∞ +∞ ∞ 1 1 k = −∞ где 1 = 1 k =1 k =1 2 , T – период повторения функции; T C −k = C k – в противном случае сумма * +∞ ∑ C k e − jk t 1 будет комплексной, а по k = −∞ условию f (t ) – действительная функция (сигнал); a0 = 1 +0.5T ∫ f (t )dt – среднее значение функции (постоянная составляющая); T −0.5T ak = 2 C k = Ak2 + Bk2 – амплитуда гармоники (множество этих величин называется спектром амплитуд); B  k = arg(C k ) = arctg k  Ak    – фаза гармоники (множество этих величин называется спектром фаз). Уместно вспомнить аналогию представления переменных токов и напряжений посредством комплексных векторов на плоскости. Комплексные коэффициенты ряда Фурье (проекции сигнала на базисные направления) находят из формулы (скалярное произведение сигнала на 1 +0.5T 2 f (t ) e − jk1t dt , где 1 = базисные функции): C k = . ∫ T −0.5T T Каждый коэффициент представляет собой амплитуду (модуль) и начальную фазу (угол) гармоник (компонент спектра). Поэтому между спектральной плотностью одиночного импульса со +∞ спектральной плотностью S ( ) = ∫ s(t ) e − jt dt и коэффициентами ряда Фурье −∞ 21 2 − j kt 1 +0.5T T ( ) Ck = s t e dt ∫ T −0.5T для периодической последовательности таких импульсов (период повторения T, который по величине больше длительности импульса) существует связь: C k = 1  2k  S  , т.е. коэффициенты ряда Фурье T  T  есть аппроксимация огибающей непрерывного спектра одиночного импульса (сигнала с бесконечным периодом повторения). Поскольку частота спектральной компоненты k определяется величиной = 2 k , то при T→∞ (при предельном переходе от периодических T прямоугольных импульсов к одиночному импульсу) шаг по частоте 2 T (расстояние между спектральными компонентами) становится бесконечно малым и соответственно число коэффициентов Ck растет до бесконечности (поскольку ширина спектра определяется шириной импульса, которая при предельном переходе остается неизменной). Выводы: − спектр непрерывного периодического сигнала дискретен (линейчатый); − спектр непрерывного непериодического сигнала непрерывен (сплошной); − спектр вещественного сигнала симметричен (поэтому область спектра с отрицательными частотами обычно не изображают). 1.4.1 Свойства преобразования Фурье +∞ Положим: X ( ) = ∫ x(t ) e −∞ − jt +∞ dt , G ( ) = ∫ g (t ) e − jt dt и Y ( ) = −∞ +∞ ∫ y(t ) e − jt dt . −∞ 1) Линейность: если y (t ) = x(t ) +  g (t ) , то Y ( ) =  X ( ) + G ( ) , ,  ∈ ℜ 2) Задержка сигнала x(t) на время τ, т.е. y (t ) = x(t −  ) : +∞ Y ( ) = ∫ x(t −  ) e −∞ − j t dt = e − j +∞ ∫ x(t −  ) e − j ( t −  ) d (t −  ) = X ( )e − j −∞ 3) Масштабирование сигнала по времени y (t ) = x (at ) : 22 +∞ Y ( ) = ∫ x(at ) e −∞ − j t  −j t 1 +∞ 1  dt = ∫ x(at ) e a d (at ) = X   , a>0, a ∈ R a −∞ a a 4) Дифференцирование сигнала y (t ) = dx(t ) : dt x(t +  ) − x(t ) − jt X ( )e j − X ( )  Y ( ) = ∫ lim = j X ( )  e dt = lim  →0  →0   −∞  +∞ +∞ 5) Интегрирование сигнала y (t ) = ∫ x(t )dt : Y ( ) = −∞ 1 X ( ) +  X (0 ) ( ) j 6) Смещение спектра Y ( ) = X ( −  0 ) : y (t ) = +∞ +∞ −∞ −∞ j t j (− )t j t j t ∫ X ( −  0 ) e d = e 0 ∫ X ( −  0 ) e 0 d ( −  0 ) = e 0 x(t ) 7) Спектр произведения сигналов – есть свертка спектров исходных сигналов. Если y (t ) = x(t ) g (t ) , то  1 +∞  j t − j t ∫  2 ∫ X ( )e d g (t ) e dt = −∞ −∞  −∞  +∞ +∞ +∞ 1 1 = X ( ) ∫ g (t ) e − j (− )t dt d = ∫ ∫ X ( )G ( −  )d = X () * G () 2 −∞ 2 −∞ −∞ +∞ Y ( ) = ∫ x(t ) g (t ) e − jt dt = +∞ 8) Спектр свертки сигналов – есть произведение спектров исходных сигналов. Если y (t ) = x(t ) * g (t ) , то +∞ +∞  +∞  − j t  − j  Y ( ) = ∫  ∫ x( ) g (t −  )d  e dt = ∫ x( )e  ∫ g (t −  ) e − j (t − )dt  d = −∞  −∞ −∞   −∞  = X ( )G ( ) +∞ 9) Спектр действительного сигнала умноженного на гармоническую функцию (частный случай произведения сигналов): y (t ) = x(t ) cos( 0 t +  0 ) +∞ Y ( ) = ∫ x(t )cos( 0 t +  0 ) e −∞ − j t [ ] 1 +∞ dt = ∫ x(t ) e j (0t +0 ) + e − j (0t +0 ) e − jt dt = 2 −∞ +∞ +∞ 1 1 = e j0 ∫ x(t )e − j (−0 )t dt + e − j0 ∫ x(t )e − j (+0 )t dt = 2 2 −∞ −∞ 1 1 = e j0 X ( −  0 ) + e − j0 X ( +  0 ) 2 2 23 X(ω) Y(ω) 0.5X(ω) ω ω 0 ωв ω0 -ω0 ω0 ω0+ωв Рис. 1.13 – Перенос спектра сигнала при умножении на гармоническую функцию 1.4.2 Связь между корреляционными функциями и спектрами сигналов Положим сигналы s1 (t ) и s2 (t ) имеют спектры S 1 ( ) и S 2 ( ) соответственно. Тогда спектр взаимной корреляционной функций +∞ B12 ( ) = ∫ s1 (t ) s2 (t −  )dt есть произведение спектров исходных сигналов: −∞ +∞ − j ∫ B12 ( ) e d = −∞ = +∞ +∞ ∫ ∫ s1 (t ) s2 (t −  ) e − j dtd = −∞ −∞ +∞ ∫ s1 (t )e −∞ − j t +∞ ∫ s2 (t −  )e − j ( t −  ) d (t −  )dt = S 1 ( ) S 2 ( ) * −∞ * 1 +∞ S 1 ( ) S 2 ( ) e jt d . или B12 ( ) = ∫ 2 −∞ Следовательно, сигналы с неперекрывающимися спектрами являются некоррелированными. 1.4.3 Примеры расчетов спектров основных сигналов 1) Спектральная плотность неинтегрируемых сигналов 1а) Постоянный сигнал: x(t ) = A ∞ X ( ) = ∫ Ae −∞ − j t A − j t dt = e − j ∞ −∞ 1 ∞ j t = 2A ( ) , т.к.  (t ) = ∫ e d 2 −∞ Постоянный сигнал в частотной области занимает бесконечно узкую полосу, т.е. спектральная плотность на нулевой частоте бесконечно велика и равна нулю на других частотах. Наличие постоянной составляющей в периодическом сигнале говорит о наличии гармоники с нулевой частотой и амплитудой равной A. 24 0, t < 0 1б) Единичный скачок (функция Хевисайда): x(t ) =  1, t ≥ 0 ∞ X () = ∫ x(t )e − jt dt = () + −∞ 1 j 1в) Гармонический сигнал: x(t ) = A cos( 0 t +  ) ∞ X ( ) = ∫ A cos( 0 t +  )e − j t −∞ = [ ] A ∞ j (  0t +  ) dt = ∫ e + e − j (0t + ) e − jt dt = 2 −∞ ∞ A j ∞ − j (  −  0 ) t A e ∫e dt + e − j ∫ e − j (+0 ) t dt = 2 2 −∞ −∞ = Ae j  ( −  0 ) + Ae − j  ( +  0 ) Физический смысл выражения заключается в том, что на частотах ±ω0 в бесконечно узкой полосе сосредоточена энергия спектральных составляющих с конечной амплитудой, поэтому спектральная плотность на этих частотах бесконечно велика. Коэффициент ряда Фурье (комплексная амплитуда гармоники) C1 гармонического сигнала связана со спектральной плотностью X ( ) : C1 = 2 1  2   A j X   = e , где  0 = , T = 2N – период анализа (N=1). T T T  2 Таким образом, каждая синусоидальная компонента периодического сигнала отражается в спектре как линия с половинной амплитудой гармоники (масштабированная дискретная дельта-функция). 2) Спектральная плотность прямоугольного импульса (видеоимпульс) с длительностью Tи: x  A, − 0.5Tи ≤ t ≤ 0.5Tи x(t ) =  0, иначе X ( ) = 0.5Tи ∫Ae −0.5Tи = − j t A − j t dt = e − j A t 0.5Tи 0.5Tи −0.5Tи = [ ] A − j 0.5Tи e − e j 0.5Tи = − j sin (0.5Tи ) 2 A  Tи  sin   = ATи  0.5Tи  2  25 3) Спектральная плотность импульсов с высокочастотным заполнением (радиоимпульсы): x(t ) = f (t ) cos( 0 t +  0 ) Положим, спектральная плотность огибающей f (t ) известна – F ( ) , тогда с учетом п.1в и свойством преобразования Фурье о спектре произведения сигналов получим: 1 1 X ( ) = e j0 F ( −  0 ) + e − j0 F ( +  0 ) 2 2 4) Спектральная плотность Гауссова импульса: x(t ) = Ee −(t  ) 2 ∞ X ( ) = ∫ Ee − (t  ) 2 −∞ = Ee    −   2  2 ∞ ∫e e − jt dt = E ∫ e −[(t  ) + jt ]dt = E ∫ e ∞ ∞ 2 −∞  t j  − +   2  −∞ 2  t j  − +   2  2 e    −   2  2 dt = −∞     2  −  t j   d + = E   e  2   2 Гауссов импульс является самым гладким сигналом, его спектральная плотность вещественна и является гауссовой функцией (быстро затухает). 5) Спектральная плотность периодического прямоугольного сигнала с амплитудой A, периодом T и длительностью импульса τ: x(t ) = A ∞  A  k   2k  + 2∑ sin    cos t. T k =1 k T   T   k  sin      T  – спектр дискретный. Коэффициенты ряда Фурье: C k = A k T  T x(t) AQ-1 sin(x ) x A 2 T τ/2 t τ T k Скважность Q=T/τ Q 2Q 3Q 4Q Рис. 1.14 – Представление периодических прямоугольных импульсов во временной и частотной областях 26 1.5 Энергетический спектр Энергия сигнала (E), по определению, равна интегралу от мощности (w) по всему интервалу существования сигнала (x): E = ∞ ∞ −∞ −∞ 2 ∫ w(t )dt = ∫ x (t )dt . Например, энергия, выделяемая на резисторе (R) при протекании тока (i) и приложенном напряжении (u) или энергия, потребленная от источника (u) T нагрузкой (R), как известно, есть E = ∫ u (t )i (t )dt = T 1T 2 ( ) u t dt = R i 2 (t )dt , ∫ ∫ R0 где T – время существования сигнала (для периодического тока – период колебаний) и p (t ) = u (t ) i (t ) – мгновенная мощность. Для возможности сравнения полагают сопротивление R=1 Ом, т.е. речь идет об удельной энергии и мощности (или нормированная мощность). На практике работают с периодическими сигналами или с сигналами конечной длительности (и бесконечной длительности, которые искусственно делают конечными на основе их «смысла» – часто простым отбрасыванием части сигнала). Поэтому с учетом необходимости сравнения энергий различных сигналов с разной длительностью проводят нормировку: E= 1T 2 ∫ x (t )dt (т.е. говорят о средней энергии), T0 1T 2 в том числе для протяженных сигналов: E = lim ∫ x (t )dt . T →∞ T 1.5.1 Теорема Рэлея Пусть S 12 ( ) – спектр взаимной корреляционной функции двух сигналов s1(t), s2(t), взаимная корреляционная функция есть: 1 +∞ B12 ( ) = S 12 ( ) e j d (обратное преобразование Фурье), ∫ 2 −∞ +∞ B12 ( ) = ∫ s1 (t ) s2 (t −  )dt (определение взаимной корреляции). −∞ Тогда при τ=0 получим выражение известное как теорема Рэлея: 27 +∞ ∫ s1 (t ) s2 (t )dt = −∞ * 1 +∞ ( ) S  S 2 ( )d . ∫ 1 2 −∞ 1.5.2 Равенство Парсеваля +∞ Энергия сигнала x(t) во временной области: E = ∫ x 2 (t )dt , энергия −∞ сигнала со спектром X ( ) в частотной области: E= +∞ 2 ∫ X () d . −∞ Таким образом, для сигнала x(t) имеющего спектр X ( ) справедливо: +∞ +∞ 2 ∫ x (t )dt = ∫ X () d . 2 −∞ −∞ Равенство Парсеваля является частным случаем теоремы Рэлея: x(t)=s1(t)=s2(t). Функция X ( ) , представляющая собой преобразование Фурье от 2 корреляционной функции, называется энергетическим спектром. Действительно, пусть сигнал x(t) имеет спектр X(ω), т.е. набор коэффициентов ряда Фурье Ck (возможно бесконечное число коэффициентов). Тогда скалярное произведение сигнала на самого себя есть норма сигнала: x(t ) = ∞ ∑ C k e j1kt , k = −∞ ∞ ∑ C m e j1mt = m = −∞ Откуда получаем: x(t ) = ∞ ∑ Cm ∞ ∞ ∑ ∑ k = −∞ m = −∞ 2 C k C m e j1kt , e j1mt = ∞ ∞ ∑ ∑ C k C m (k − m ) k = −∞ m = −∞ . m = −∞ ∞ С другой стороны норма есть: x(t ) = ∫ x 2 (t )dt . −∞ ∞ Следовательно: 2 ∫ x (t )dt = −∞ ∞ ∑ Cm m = −∞ 2 , коэффициенты Cm есть спектр сигнала. Таким образом, мощность сигнала есть сумма квадратов модулей коэффициентов ряда Фурье. 28 1.5.3 Принцип неопределенности времячастотного представления Очевидно, что в соответствии со свойствами время-частотного представления сигнала при «сжатии функции f» по оси времени (t) в A раз, т.е. рассмотрении функции fA(t)=f(At), ее спектр растянется во столько же раз: FA(k)=const·F(k/A), поскольку частоты каждой спектральной гармоники ejkt этого разложения должны, очевидно, умножиться на A. Строго говоря, пример носит довольно частный характер, однако показывает физический смысл: когда мы сжимаем сигнал, его частотная полоса пропорционально увеличивается (верно и обратное). Это свойство подобно известному неравенству Гейзенберга: x  > h m , где Δx – неопределенность (погрешность измерения) пространственной координаты микрочастицы, Δv – неопределенность скорости частицы, m – масса частицы, а h – постоянная Планка. Поэтому нельзя одновременно измерить энергию системы и указать момент времени, которому это измерение соответствует, или нельзя одновременно локализовать сигнал во временной и в частотной областях. Теорема. Для любого дифференцируемого вещественного сигнала с энергией (E) произведение полосы (f) и длительности (t) ограничено снизу: f t ≥ E . 4 Поэтому в реальности невозможно измерить спектр ограниченного во времени сигнала (конечной длительности), т.к. такой спектр имеет бесконечно много составляющих и его приходится ограничивать искусственно. Соответственно для бесконечно протяженных сигналов (в т.ч. периодических) спектр ограничен. 29 1.6 Преобразование Гильберта Аналитическим сигналом называется сигнал s a (t ) = s (t ) + js ⊥ (t ) , где s(t ) – сигнал; s⊥ (t ) – квадратурное дополнение (сопряжение сигнала). Сопряженный сигнал связан с исходным сигналом преобразованием 1 ∞ s( ) Гильберта: s⊥ (t ) = ∫ d  −∞ t −  (преобразование – есть свертка сигнала s(t ) с функцией 1 ). t Частотная характеристика преобразования:  j,  < 0 1  K ⊥ ( ) = ∫ e − jt d = 0,  = 0 , −∞ t − j ,  > 0  ∞ т.е. преобразование Гильберта – идеальный фазовращатель. В частности, если на входе преобразователя Гильберта действительный сигнал, то и на выходе тоже действительный сигнал. Например: если s (t ) = cos ( 0 t ) , то s ⊥ (t ) = sin ( 0 t ) . Обратное преобразование: s(t ) = − 1 ∞ s⊥ ( ) d . ∫  −∞ t −  Рассчитаем спектр аналитического сигнала s a (t ) = s (t ) + js ⊥ (t ) : 0,  < 0  S a ( ) = S ( ) + j S ⊥ ( ) = S ( )[1 + j K ⊥ ( )] = S (0 ),  = 0 2 S ( ),  > 0  30 Im Re Re Re 0.5 0.5 -ω0 Im Im 1.0 -ω0 ω0 ω0 ω0 -0.5 Частота Частота s ⊥ (t ) = sin ( 0 t ) s (t ) = cos(0t ) Частота s a (t ) = s (t ) + js ⊥ (t ) Рис. 1.15 – Графическая интерпретация спектра аналитического сигнала s⊥ (t ) = sin ( 0t ) s a (t ) = s (t ) + js⊥ (t ) s (t ) = cos(0t ) Рис. 1.16 – Графическая интерпретация аналитического сигнала во временной области Таким s(t ) = cos( 0 t ) образом, симметричный становится спектр односторонним, действительного совпадающим со сигнала спектром исходного сигнала (для положительных частот, с точностью до множителя 2.0). 31 Важными свойствами аналитического сигнала является возможность определения: 1) огибающей сигнала A(t ) = s a (t ) = s 2 (t ) + s⊥2 (t )   s⊥ (t )  , если s(t ) ≥ 0 arctg   s(t )  2) полной фазы сигнала  (t ) = arg(s a (t )) =   + arctg s⊥ (t ) , если s(t ) < 0   s(t )  Например, рассмотрим амплитудно-модулированный сигнал: s (t ) = [0.4 cos( u t )]cos( 0 t ) , где ω0 – несущая частота; ωu – частота модулирующего сигнала (ω0>> ωu). Тогда квадратурное дополнение: s ⊥ (t ) = [0.4 cos( u t )]sin ( 0 t ) и огибающая сигнала: A(t ) = 0.4 cos( u t ) . Im Im Re Re -(ω0+ωu) -ω0 -(ω0-ωu) -(ω0+ωu) -ω0 -(ω0-ωu) (ω0-ωu) S ( ) ω0 (ω0+ωu) Частота (ω0-ωu) ω0 (ω0+ωu) S ⊥ ( ) Частота Рис. 1.17 – Графическая интерпретация спектра амплитудномодулированного сигнала 32 s a (t ) s⊥ (t ) Im A(t ) = s 2 (t ) + s⊥2 (t ) Re s (t ) время Рис. 1.18 – Графическая интерпретация амплитудно-модулированного сигнала во временной области (аналитический сигнал и огибающая) Области применения преобразования Гильберта: − модуляция и демодуляция (связь); − автоматическая регулировка усиления; − оценка мгновенной частоты, измерение задержки сигналов; − анализ двухмерных, трехмерных сигналов; − сжатие изображений и аудиосигналов; − обработка сигналов в радарах, сонарах, телевидении высокой четкости и т.п. 33 1.7 Понятия о модуляции сигналов Одной из важнейших проблем обработки сигналов является задача передачи сигнала по каналам связи. Здесь можно выделить два важных аспекта: 1) рациональное использование характеристик канала связи (задача передачи как можно большего количества информации за минимальное время и в узком частотном диапазоне); 2) передача данных без искажения (помехоустойчивость передачи). Обе задачи крайне важны с практической точки зрения, т.к. фактически эффективность их решения обуславливает надежность и стоимость передачи. Например, витая пара 5 класса (как среда передачи сигнала) имеет полосу частот порядка 300 МГц. Допустим необходимо передавать сигнал с полосой порядка 4 кГц. Тогда если генератор такого сигнала подключить к витой паре, то приемник естественно получит этот сигнал. Однако возможности кабеля с его полосой пропускания оказываются не востребованными. Кроме того, если бы необходимо передавать много сигналов с полосой 4 кГц (например: задача городской телефонии). Тогда очевидно, возникает задача: как «разместить» низкочастотные сигналы в полосе частот кабеля так, чтобы они не влияли (не искажали) друг на друга. В частности, для решения подобных задач необходимо уметь сдвинуть спектры сигнала в другую область частот, для чего, согласно свойствам преобразования Фурье, достаточно умножить сигнал на гармоническую функцию – рис 1.19. x1 ×cos(ω1t) ω x2 в канал передачи сигналов ×cos(ω2t) ω + y ……… x1 x2 xk xk …… ×cos(ωkt) ω1 ω2 ωk ω ω Рис. 1.19 – Схематичное представление частотного уплотнения канала 34 Фактически рассмотренный процесс и есть частный случай модуляции сигнала. Модуляция сигнала в общем случае рассматривается как изменение одного параметра сигнала другим. В самом распространенном случае модуляции подвергается гармонический сигнал с частотой ω0, который называется несущей. По характеру воздействия на несущий синусоидальный сигнал выделяют четыре вида модуляции (которые нашли практическое применение): амплитудная, частотная, фазовая и квадратурная. При этом сигнал, изменяющий параметр несущей, называется модулирующим. Сам процесс изменения – модуляцией, обратной операцией является – демодуляция. Типичными практическими задачами, в которых не обойтись без модуляции, является радиопередача информации (сотовая связь, телевидение, радары и т.п.), поскольку информационные сигналы – это широкополосные низкочастотные сигналы, которые по физическим законам неудобны для передачи антеннами с малыми габаритами и приемлемым КПД системы приема/передачи. Также при радиопередаче данных надо учитывать свойства свободного пространства – есть частоты, на которых наблюдается сильное затухание колебаний, например, в атмосфере, гидросфере и т.п. Поэтому не всякий диапазон частот пригоден для осуществления радиопередачи. 1.7.1 Общие сведения об амплитудной модуляции Наиболее простой для понимания является амплитудная модуляция, которая подразумевает изменение амплитуды несущей информационным сигналом A(t): s (t ) = A(t )cos ( 0 t ) . Одним из условий осуществления модуляции является требование, что частота несущей много больше частоты (полосы) сигнала. Рассмотрим пример. Пусть A(t ) = sin (1t ) + 0.7 sin (21t ) . Построим сигнал s (t ) = [sin (1t ) + 0.7 sin (21t )]cos ( 0 t ) – рис. 1.19. 35 Рис. 1.20 – Проблема неоднозначности при АМ модуляции Из рисунка видно, что выделение огибающей для восстановления передаваемого сигнала A(t) имеет неоднозначный характер (зеленая огибающая не соответствует передаваемому сигналу – показан красным цветом). Чтобы избежать неоднозначности (говорят перемодуляции) при амплитудной модуляции двухполярным сигналом вводят смещение (A0) и масштабирующий множитель (ξ): s (t ) = [ A0 + A(t )]cos ( 0 t ) . На рис. 1.21 представлен результат амплитудной модуляции того же сигнала (что и на рис. 1.20), но при A0=1, ξ=0.5. Соответственно, огибающая (показана красным цветом) формируется однозначным образом. Рис. 1.21 – Амплитудная модуляция без неоднозначности Если информационный сигнал A(t) есть гармоническая функция, то такую амплитудную модуляцию называют однотональной: 36 s (t ) = [ A0 + Am cos(1t + 1 )]cos( 0 t ) . При этом вводится коэффициент m = Am , называемый глубиной модуляции, A0 значение которого можно связать с минимальным и максимальным значениями модулирующего сигнала: m = Amax − Amin . Amax + Amin Можно записать (на рис.1.22 представлен сигнал во временной области): s(t ) = [ A0 + Am cos(1t + 1 )]cos( 0 t +  0 ) = = A0 cos( 0 t +  0 ) + Am cos(1t + 1 )cos( 0 t +  0 ) = = A0 cos( 0 t +  0 ) + 0.5 Am cos(( 0 − 1 )t +  0 − 1 ) + + 0.5 Am cos(( 0 + 1 )t +  0 + 1 ) = = A0 cos( 0 t +  0 ) + mA0 mA0 cos(( 0 − 1 )t +  0 − 1 ) + cos(( 0 + 1 )t +  0 + 1 ) 2 2 Таким образом, в спектре сигнала однотональной амплитудной модуляции будет присутствовать три спектральных пика – рис. 1.23 (ω0>>ω1). В общем случае в частотной области АМ сигнал содержит несущую частоту и две боковые полосы: − верхнюю с частотой ω0+ω1, представляет собой копию спектра информационного (модулирующего) сигнала; − нижнюю с частотой ω0-ω1, представляет собой зеркальную копию спектра информационного (модулирующего) сигнала. Рис. 1.22 – Однотональная АМ сигнала (A0=1, m=0.4) 37 |S| A0 0.5A0m ω0-ω1 0.5A0m ω0 ω0+ω1 arg(S) φ0+φ1 φ0 φ0-φ1 ω0-ω1 ω0 ω0+ω1 Рис. 1.23 – Амплитудный и фазовый спектры однотональной АМ сигнала Вычислим мощность однотональной АМ сигнала: мощность несущей (квадрат действующего значения) 0.5A2, мощность каждой боковой полосы 0.125m2A2. Итого мощность составляет: 0.5A2+0.25m2A2. Оценим энергетическую эффективность как отношение мощности боковых частот к 0.25m 2 A 2 m2 общей мощности сигнала: . Откуда можно видеть, = 0.5 A 2 + 0.25m 2 A 2 m 2 + 2 что в предельном случае m=1 эффективность не превосходит одной трети (в лучшем случае), т.е. более двух третей мощности тратиться на передачу бесполезной (с точки зрения информации) несущей. Это является одним из серьезных недостатков амплитудной модуляции. Другим недостатком является широкая полоса сигнала – в два раза большая, чем полоса информационного сигнала. Действительно, если принять во внимание что модулирующий сигнал содержит множество гармоник (имеет сложную форму), т.к. однотональная АМ на практике не применяется, становится очевидным, что ширина спектра АМ сигнала в два раза больше ширины спектра модулирующего сигнала. Поэтому существует несколько модификаций амплитудной модуляции: − для снижения ширины полосы частот канала одну из боковых полос подавляют (называется однополосой АМ); 38 − из энергетических соображений подавляют несущую (если приемник точно «знает» величину ω0). Однако эти усилия не могут преодолеть еще один недостаток амплитудной модуляции – плохая помехоустойчивость. Для демодуляции АМ сигнала можно имитировать работу двухполупериодного детектора – для этого необходимо вычислить модуль АМ сигнала и пропустить его через низкочастотный фильтр (сгладить косинусные импульсы – модуль кривой рис.1.22). Другим подходом является синхронное детектирование, суть которого состоит в умножении АМ сигнала на опорную функцию с несущей частотой – рис. 1.24: y (t ) = s(t )cos( 0 t +  0 ) = A(t )cos 2 ( 0 t +  0 ) = 0.5 A(t ) + 0.5 A(t )cos(2 0 t + 2 0 ) где s (t ) = [ A0 + Am cos (1t + 1 )]cos ( 0 t +  0 ) = A(t ) cos ( 0 t +  0 ) s(t) y(t) ФНЧ × A(t) cos(0t +  0 ) Рис. 1.24 – Синхронный детектор Синхронное детектирование требует точного знания частоты и фазы несущего колебания, что осложняет практическую реализацию. 1.7.2 Общие сведения об угловой модуляции Сигналами с угловой модуляцией называются частотно-модулированные (ЧМ) или фазомодулированные (ФМ) сигналы. В общем виде колебания с угловой модуляцией описываются Запишем полную фазу сигнала выражением: u (t ) = U m sin ( 0 t + (t )) . (t ) =  0 t + (t ) . Мгновенной частотой называется производная полной фазы: (t ) = d (t ) d (t ) . = 0 + dt dt При частотной модуляции частота сигнала изменяется в соответствии с информационным сигналом s(t): (t ) =  0 + ks (t ) . Поэтому ЧМ сигнал можно выразить в форме: 39 t   u (t ) = U m sin   0 t + k ∫ s(t ) , −∞   т.е. сдвиг фазы сигнала θ(t) связан с информационным сигналом через интеграл. Максимальное отклонение частоты сигнала от величины ω0 называется девиацией частоты:  = k max s(t ) . На рис. 1.25 представлены сигнал треугольной формы – информационный сигнал s(t), соответствующая ему t полная фаза (t ) =  0 t + k ∫ s(t ) и ЧМ сигнал x(t ) = sin ((t )) . −∞ s(t) t ψ(t) t x(t) t Рис. 1.25 – ЧМ сигнала В частном случае, при однотональной модуляции, когда модулирующий сигнал есть гармоническая функция, выражение для мгновенной частоты принимает вид: (t ) =  0 +  cos ( t +  ) , где Δω – девиация частоты, ω0 – частота несущей, Ω и Φ – частота и начальная фаза информационного сигнала. 40 Тогда полная фаза определяется выражением: t (t ) =  0 t + k ∫ s(t ) =  0 t + −∞ Величина m =  sin (t +  ) .   называется индексом модуляции, который характеризует  максимальное отклонение сдвига фазы от среднего значения. При фазовой модуляции сдвиг фаз изменяется пропорционально информационному сигналу: (t ) = ks(t ) , поэтому полная фаза: (t ) =  0 t + ks(t ) . На рис. 1.26 представлены информационный сигнал s(t), тот же – рис. 1.25, полная фаза (t ) =  0 t + ks(t ) и ФМ сигнал x(t ) = sin ((t )) . s(t) t ψ(t) t x(t) t Рис. 1.26 – ФМ сигнала 41 Сравнивая рис. 1.25 и рис. 1.26 можно видеть различия в ЧМ и ФМ сигналах. В общем случае, трудно по виду модулированного сигнала назвать вид угловой модуляции. В частном случае, при однотональной фазовой модуляции выражения принимают вид: (t ) = m cos(t +  ) , где m – индекс модуляции; (t ) =  0 + d (t ) =  0 − m sin (t +  ) , т.е. максимальное отклонение частоты dt равно  = m (совпадает со случаем ЧМ). В заключение подраздела, рассмотрим спектр сигналов с угловой модуляцией. Найдем Фурье коэффициенты для сигнала с угловой однотональной модуляцией: u (t ) = U max cos ( 0 t + m cos ( t +  ) +  0 ) . Однако вместо непосредственного интегрирования (расчета преобразования Фурье), e jm cos ( x ) = учтем симметричность сигнала и следующее выражение: ∞ ∑ J n (m )e jnx , n = −∞ Jn(m) – функция Бесселя n-го порядка. Тогда для сигнала справедливо: ∞   u (t ) = U max Re{e j (0t +0 )e jm cos (t + ) } = U max Re e j (0t +0 ) ∑ J n (m )e jn (t + )  = n = −∞   = U max ∞ ∑ J n (m )cos(( 0 + n )t +  0 + n ) n = −∞ Очевидно, что если удалось выразить симметричный сигнал через сумму гармонических составляющих (косинусных!), то их амплитуды есть модуль спектра, а сдвиги фаз есть начальные фазы спектральных составляющих. Анализируя поведение функций Бесселя при n>m+1, оказывается, что составляющие с частотами  0 + n становятся пренебрежимо малыми. Тогда изначально бесконечную ширину спектра сигнала с угловой модуляцией можно рассматривать как ограниченную величиной 2 + 2 (Δω – девиация частоты). 42 При узкополосой угловой модуляции (m<<1) функции Бесселя можно определить по приближенным формулам: J 0 (m ) ≈ 1, J 1 (m ) ≈ 0.5m , J 2 (m ) ≈ 0.125m 2 . Тогда спектр ЧМ сигнала содержит три составляющих (остальными можно пренебречь): несущую и две боковые (рис. 1.27) – подобно спектру АМ сигнала (разница: фаза нижней боковой составляющей противоположна). |S| Umax |-0.5mUmax| ω0-Ω 0.5mUmax ω0 ω0+Ω Рис. 1.27 – Амплитудный спектр угловой модуляции сигнала при m<<1 Преимущества сигналов с угловой модуляцией состоят в более высокой помехоустойчивости и большем динамическом диапазоне (в отличие от АМ нет эффекта перемодуляции, т.е. нелинейных искажений сигнала при росте индекса модуляции). 43 1.8 Случайные сигналы Принципиальным отличием от детерминированных сигналов является то, что мгновенные значения случайных сигналов неизвестны и могут быть лишь предсказаны с некоторой вероятностью. Характеристики случайных сигналов являются статистическими (имеют вероятностный характер). В обработке сигналов выделяют два класса случайных сигналов: шумы и информационные сигналы. Шумам присуще хаотическое изменение сигнала во времени (чаще всего связанное с беспорядочным движением носителей заряда в физических системах). Информационные случайные сигналы представляют собой осмысленные сообщения, получаемые какой-либо статистической моделью. Строго математически, случайный процесс (или случайный сигнал) x(t ) – это функция, характеризующаяся тем, что значения, принимаемые ею в любой момент времени t, являются случайными величинами. До регистрации (приема) случайного сигнала, его следует рассматривать как случайный процесс, который представляет собой ансамбль функций времени. Одна из этих функций после регистрации (приема) сигнала становится известной и называется реализацией случайного процесса. Типичными примерами случайного сигнала являются звуковые колебания (например: речь другого относительно нас человека, шум водопада), колебания температуры в какой-либо местности, сила ветра, колебания частоты в городской электросети, уровень пульсаций напряжения на выходе генератора, ошибка округления чисел при представлении в конечной разрядной сетке и т.п. 1.8.1 Вероятностные характеристики случайных сигналов Пусть X (t ) – случайный процесс, заданный ансамблем реализаций { x1 (t ) , x 2 (t ) , …, x k (t ) , …}. Выберем произвольный момент времени t1 и получим совокупность (вектор) { x1 (t1 ) , x 2 (t1 ) , …, x k (t1 ) , …}, представляющую собой случайную величину X (t1 ) . 44 Функция распределения вероятности (cumulative distribution function), обозначаемая как F ( x, t1 ) , равна вероятности того, что в момент времени t1 значение случайного процесса не превосходит x: F ( x, t1 ) = P( X (t1 ) ≤ x ) . Функция F ( x, t1 ) является неубывающей и лежит в диапазоне [0, 1]. В частности: F (− ∞, t1 ) = 0 и F (∞, t1 ) = 1 . Вероятность попадания значения случайного процесса в интервал (a, b] равна разности значений функции распределения на концах этого интервала: P(a < X (t1 ) ≤ b ) = F (b, t1 ) − F (a, t1 ) . Плотность вероятности представляет собой производную от функции распределения: p( x, t1 ) = ∂F ( x, t1 ) . ∂x Математическое ожидание случайного процесса: +∞ m x (t ) = M {X (t )} = ∫ x p( x, t )dx . −∞ { } +∞ Дисперсия Dx (t ) = M [ X (t ) − m x (t )] = ∫ x 2 p( x, t )dx − m x2 (t ) . 2 −∞ Среднее квадратичное отклонение (СКО):  x (t ) = Dx (t ) . Некоторые типовые распространенные распределения и их характеристики. 1. Равномерное распределение 0, x < a  x − a F (x ) =  , a≤ x≤b b − a  1, x > b  1 , a≤ x≤b  , p( x ) =  b − a 0, иначе x a+b dx = 2 ab−a b m x (t ) = ∫ x2  a + b  (b − a ) Dx (t ) = ∫ dx −   = 12  2  ab−a b 2 2 2. Нормальное распределение (Гауссов закон распределения) 45 p( x ) = 1 x 2     x − m x   exp − 0.5  2   x     (x − mx )   F ( x ) = Ф  2 x  1 Интеграл вероятности Ф( x ) = 2 x −0.5t ∫ e dt – табулированная функция. 2 −∞ 1.8.2 Корреляционные функции случайных процессов Ковариационная функция случайного процесса – это математическое ожидание произведения реализаций случайного процесса (или случайных процессов): K (t1 , t 2 ) = M {x(t1 ) x(t 2 )} = ∞ ∞ ∫ ∫ x1 x2 p(x1 , x2 , t1 , t 2 )dx1dx2 −∞−∞ где p( x1 , x2 , t1 , t 2 ) – двухмерная плотность вероятности, определяемая как dx dx   p( x1 , x2 , t1 , t 2 )dx1dx2 = P  X (t1 ) − x1 ≤ 1 , X (t 2 ) − x2 ≤ 2  – вероятность того, 2 2   что реализация случайного процесса X (t ) в момент времени t1 попадает в бесконечно малый интервал dx1 в окрестности x1, и в момент времени t2 – в бесконечно малый интервал dx2 в окрестности x2. Корреляционная составляющую) есть функция (характеризует статистически флуктуационную усредненное произведение центрированной случайной функции X (t ) − m x (t ) в моменты времени t1 и t2: R(t1 , t 2 ) = M {[x(t1 ) − mx (t1 )] [x(t 2 ) − mx (t 2 )]} = K (t1 , t 2 ) − mx (t1 )mx (t 2 ) . Если t1=t2=t (при совмещении сечений случайного процесса): R(t , t ) = Dx (t ) . Рассматривая два случайных процесса (или две реализации одного процесса) между ними может существовать либо отсутствовать статистическая связь. Отсутствие связи говорит о независимости одной случайной величины от другой. Тогда двумерная плотность вероятности: p( x1 , x2 ) = p( x1 ) p( x2 ) . 46 Это выражение известно как условие статистической независимости. Мерой линейной статистической связи является коэффициент корреляции: r12 = M {x1 x2 } − M {x1}M {x2 } , очевидно r12 ≤ 1 . D{x1}D{x2 } Если r12 = 0 , то между случайными величинами отсутствует статистическая связь, т.е. они некоррелированны, при этом M {x1 x2 } = M {x1} M {x2 }. 1.8.3 Стационарные процессы В общем случае вероятностные характеристики случайных процессов зависят от времени (моментов выборки, т.е. получения реализации). Однако существует обширный класс процессов, у которых вероятностные характеристики не зависят от времени, т.е. одинаковы у различных реализаций случайного процесса. Такие процессы называются стационарными. Строго говоря, случайный процесс стационарен, если его многомерная плотность вероятности p( x1 , x2 ,  , xn , t1 , t 2 ,  t n ) не изменяется при одновременном сдвиге временных сечений tk вдоль оси времени на произвольную величину τ. Таким образом, для стационарного сигнала математическое ожидание и дисперсия не зависят от времени, а корреляционная функция зависит только от интервала  = t 2 − t1 : R (t1 , t 2 ) = R ( ) , причем R( ) = R(−  ) и R( ) ≤ R(0 ) = Dx . Стационарный процесс называется эргодическим, если при определении его статистических характеристик усреднение по ансамблю реализаций эквивалентно усреднению по времени одной (любой) реализации. Тогда математическое ожидание и дисперсию эргодического процесса можно определить через бесконечно длинную реализацию: +∞ 1 0.5T m x = M {X (t )} = ∫ x p( x )dx = lim ∫ x(t )dt T →∞ T −∞ −0.5T Dx = +∞ 1 0.5T 2 2 ∫ x (t )dt − mx T →∞ T −0.5T 2 2 ∫ x p(x )dx − mx = lim −∞ 47 Таким образом, математическое ожидание эргодического процесса есть постоянная составляющая его реализации, а смысл дисперсии – мощность переменной составляющей. Достаточным условием эргодичности случайного процесса является стремление к нулю его корреляционной функции с ростом временного сдвига τ: lim R( ) = 0 .  →∞ 1.8.4 Спектральные характеристики случайных процессов Каждая реализация случайного процесса представляет собой уже известную функцию (измеренную зависимость) к которой можно применить преобразование Фурье. При этом, конечно каждая реализация будет иметь свой спектр. На практике интересуются усредненными характеристиками: ∞ S ср ( ) = ∫ xср (t )e − j t −∞ ∞ dt = ∫ m x (t )e − jt dt −∞ Таким образом, усредненный спектр случайного процесса есть спектр его математического ожидания. В частности, если m x (t ) = 0 то и S ср () = 0 (но S ср () ≠ 0 – т.е. фазы спектральных составляющих в различных реализациях процесса случайны и независимы). Случайный сигнал во временной области обладает энергией, которая равна энергии в частотном представлении сигнала – в соответствии с равенством Парсеваля. Мысленно ограничивая случайный процесс интервалом [-0.5T, 0.5T] и разделив энергию на величину T, получим среднюю мощность: PT = 1 0.5T 2 1 x (t )e − jt dt = ∫ T −0.5T 2T ∞ ∫ X () 2 d −∞ Далее совершая предельный переход при T → ∞ , получим выражение спектральной плотности средней мощности реализации x(t ) : X ( ) W ( ) = lim и соответственно T →∞ T 2 ∞ 1 ∞ ∫ x (t )dt = 2 ∫ W ()d . −∞ −∞ 2 48 Однако на практике этими формулами не пользуются, т.к. часто модель случайного процесса не позволяет произвести непосредственное вычисление предела. 1.8.5 Теорема Винера-Хинчина Как известно корреляционная функция сигнала связана с его преобразованием Фурье (рассматривая сигнал длительностью T): 0.5T ∫ x(t )x(t −  )dt = −0.5T 1 ∞ 2 X ( ) e j d ∫ 2 −∞ Разделив обе части уравнения на T и переходя в пределе при T → ∞ , получим: X ( ) j 1 0.5T 1 ∞ ( ) ( ) lim x t x t −  dt = lim e d ∫ ∫ T →∞ T 2  − ∞ T →∞ T −0.5T 2 Теорема: Корреляционная функция случайного процесса и его спектральная плотность мощности связаны преобразованием Фурье. Таким образом, для эргодического процесса R( ) = 1 ∞ W ( )e j d . ∫ 2 −∞ Учитывая четность и вещественность функций R(τ) и W(ω) можно записать: ∞ 1∞ R( ) = ∫ W ( )cos( )d и W ( ) = 2 ∫ R( )cos( )d . 0 1.8.6 Интервал корреляции Рассматривая корреляционную функцию случайного процесса можно отметить, что она убывает с ростом величины сдвига τ. Причем скорость затухания функции характеризует степень статистической связи между значениями случайного сигнала в разные моменты времени. Числовой характеристикой изменения реализаций случайного процесса является величина интервала корреляции: 1 ∞ к = ∫ R( ) d . R(0 ) 0 49 1.8.7 Эффективная ширина спектра случайного сигнала Пусть имеется случайный сигнал со спектром плотности мощности W(ω), который имеет максимальное значение Wmax. Мысленно представим себе другой случайный процесс, у которого спектральная плотность мощности постоянна и равна Wmax в некоторой (очевидно более узкой полосе) полосе частот. Ширина этой полосы называется эффективной шириной спектра случайного процесса:  эф 1 ∞ = ∫ W ()d . Wmax −∞ (*) Причем величины τк и Δωэф связаны соотношением неопределенности:  эф  к ≈ 2 . 1.8.8 Белый шум Белый шум — стационарный сигнал, спектральные составляющие которого равномерно распределены по всему диапазону частот, т.е. W(ω)= W0=const. Поэтому корреляционная R( ) = в функция соответствии с теоремой белого шума есть Винера-Хинчина дельта функция: W0 ∞ j ∫ e d = W0 ( ) , 2 0 т.е. дисперсия (средняя мощность) белого шума бесконечно велика, другими словами любые два значения сигнала некоррелированны при сколь угодно амплитуда малом τ. На рис. 1.28 представлена реализация белого шума. время Рис. 1.28 – Реализация белого шума во временной области 50 Белый шум является абстрактной моделью. В природе и технике не встречается, однако под категорию белых шумов попадают любые шумы, спектральная плотность которых одинакова (или слабо отличается, в этом случае говорят об окрашенном шуме) в рассматриваемом диапазоне частот (говорят: если эффективная полоса частот шума превышает полосу частот устройства, тогда для этого устройства шум удобно считать белым). Например, рассмотрим единичный импульс длительностью 1 мкс. Тогда, выполнив преобразование Фурье, получим огибающую спектра, описываемую модулем sinc функцией, первый ноль которой наблюдается на частоте 1 МГц – рис.1.29. Следовательно, для устройства с полосой частот не превышающей, к примеру, 20 кГц, спектр данного импульса практически постоянен, т.е. его можно считать моделью белого шума. f, кГц f, МГц Рис. 1.29 – Огибающая спектра импульса с длительностью 1 мкс 51 1.9 Программная среда Matlab Программный различного пакет назначения Применительно к Matlab систем цифровой и предназначен проведения обработке для моделирования численных сигналов пакет расчетов. представляет возможности моделирования цифровых систем обработки информации, проектирования блоков алгоритмов, вычисления характеристик и параметров алгоритмов. интерактивного Процесс моделирования визуального программирования на средства обеспечивается – матричном языке Simulink Matlab или посредством с (имеется помощью встроенный компилятор). Сам язык является высокоуровневым языком, приспособленным для выполнения матричных и векторных операций, в том числе с реализацией сложных распространенных алгоритмов: триангуляция и обращение матриц, решение систем уравнений, методы оптимизации, различные типовые преобразования и т.п., а также реализацией графического интерфейса для представления данных и результатов. Как любой язык программирования, Matlab включает операторы управления выполнением программы и собственно операторы, функции осуществляющие вычисления, преобразования данных. Вид программной строки (символ % означает комментарий, всякая запись после этого символа игнорируется): <оператор>; %комментарий <переменная>=<значение, вычисления или функция с набором параметров>; Оператор цикла. for <имя переменной>=<начальное значение>:<конечное значение> <тело цикла>; %переменная получает приращение от нач. до кон. значения, шаг 1 end; % число повторений тела цикла =(кон. значение – нач. значение +1) for <имя переменной>=< нач. значение >:<шаг>:<конеч. значение> <тело цикла>; переменная изменяется от нач. до кон. значения с заданным шагом end; Условный оператор. 52 if (<условие>) <операторы, выполняемые если условие выполняется>; else <операторы, выполняемые если условие не выполняется>; end; Пример 1. 5 Вычислим сумму ряда: S = ∑ 2n n=2 S=0; % инициализация переменной, присвоение начального значения for n=2:5 % переменная n последовательно будет равна 2, 3, 4, 5 S=S+2*n; end; % S – содержит результат Работа с массивами (векторами, матрицами) Инициализация массива: x=[1 2 3 4]; % одномерный массив – вектор (строка) с координатами (1,2,3,4) 1 2 3  a=[1 2 3; 4 5 6]; % двухмерный массив – матрица a =    4 5 6 t=[1:3:11]; % задает вектор (строка) t=[1 4 7 10] t=(1:3:11); %тоже самое: t=[1 4 7 10] t1=t’; %транспонирование, например: преобразуем строку в столбец Например, альтернативная форма программы примера 1: n=(2:1:5); % эквивалентно n=[2 3 4 5] S=2*sum(n); % операция SUM возвращает сумму элементов вектора (массива) Обращение к элементам массива (индексы нумеруются с единицы) q=x(2); % вторая координата вектора: q=2 w=a(2,1); % первый индекс – строка, второй – столбец: w=4 b=a(1, :); % первая строка матрицы a: b=[1 2 3] с=(2, 2:3); % элементы 2 и 3 столбцов во второй строке: c=[5 6] (*) команды присвоения, изменения элемента массива – аналогично. 53 Пример 2. Вычислить скалярное произведение векторов x1=(1; 2) и x2=(3; 5). По определению y = ∑ x1 (k )x2 (k ) , тогда y = 1 ⋅ 3 + 2 ⋅ 5 = 13 . k x1=[1 2]; % инициализация массивов x2=[3 5]; y=0; % процедура вычисления скалярного произведения for j=1:2 y=y+x1[j]*x2[j]; end; Альтернативные формы программы с использованием матричной алгебры: a) y=x1*x2’; % умножить строку на столбец есть число b) y=sum(x1.*x2); %поэлементно перемножить строки и сложить результаты % операции поэлементного умножения, деления, возведения в степень предваряются знаком точки: .* ./ .^ c) Программный код процедуры эквивалентен встроенной в Matlab функции: y=dot(x1,x2); % длины векторов должны быть одинаковыми Пример 3. Реализовать линейную свертку векторов (сигналов) h(k) длины K и входного сигнала x(n). y=0; for j=1:K m=n-j+1; if (m>0) y=y+h(j)*x(m); % (*) end; end; % y – содержит значение выходного отсчета y(n) (*) необходимо следить, чтобы индекс массива был положительным, в противном случае необходимо дополнить нулями вектор x (вставить в массив слева K-1 нулей). Функция, встроенная в Matlab: y=conv(h,x); %len(y)=len(h)+len(x)-1 54 Пример 4. Рассчитать спектр сигнала: x(n ) = sin (2f 0 kt s ) , f0=50 Гц, ts-1=800 Гц clear; % очистка переменных, освобождение памяти f0=50; % частота сигнала fs=800; % частота дискретизации N=128; % длина ДПФ x=sin(2*pi*f0/fs*(0:1:N-1)); % массив отсчетов заданного сигнала y=fft(x,N); % расчет комплексного спектра сигнала yy=2/N*abs(y); % расчет модуля спектра % yy – амплитудный спектр сигнала, частотный отсчет 50/800*128=8 (нумерация отсчетов сигнала – с нуля) соответствует заданной гармонике: yy(k=9)=1 (нумерация отсчетов в массиве – с единицы), все остальные отсчеты при 1≤k≤64, k≠9: yy(k)=0. Отображение результата. Команда plot. Для построения зависимостей (в виде графиков, диаграмм) применяется команда plot(x,y,) – для представления результатов на плоскости и plot3(x,y,z,) – для представления результатов в объеме. - область входных данных команд, устанавливающая атрибуты кривой, такие как цвет, тип линии, цвет и тип маркеров и т.п. В общем случае для отображения графика необходимо иметь сами данные в виде массивов x – координаты по оси абсцисс, y – координаты по оси ординат (естественно длины массивов x и y одинаковы, и каждая пара (x,y) задает точку на плоскости): % дополнение к примеру 4 figure; %показать окно области построения grid on; %включить сетку hold on; %показывать на одной области построения несколько графиков % отобразим спектр примера 4, по оси абсцисс номера отсчетов plot((0:1:N-1),yy, 'r'); %построить график красным цветом % отобразим спектр примера 4, по оси абсцисс значения частоты plot((0:1:63)/128*fs,yy(1:64), 'b'); %построить половину графика синим цветом 55 Математические функции abs(x) – модуль числа x (для комплексного значения – длина вектора) sqrt(x) – вычисляет квадратный корень числа x fix(x) – выполняет округление числа x методом отбрасывания дробной части round(x) – округление числа x до ближайшего целого sin(x) – вычисляет синус угла x (задается в радианах) cos(x) – вычисляет косинус угла x (задается в радианах) exp(x) – вычисляет значение функции ex log(x) – вычисляет значение функции ln(x) sinc(x) – вычисляет значение функции sin( x) x Функции формирования сигналов y=rectpuls(t,w) – формирует одиночный прямоугольный импульс с единичной амплитудой (t – вектор значений времени, w – ширина 1, − 0.5w ≤ t < 0.5w импульса): y =  0, t < −0.5w, t ≥ 0.5w y=tripuls(t,w,s) – формирует одиночный треугольный импульс с единичной амплитудой (t – вектор значений времени, w – длительность импульса, s – параметр асимметрии: положение максимума  2t + w  w( s + 1) , − 0.5w ≤ t < 0.5ws   2t − w 0.5ws, -1≤s≤1): y =  , 0.5ws ≤ t < 0.5w w ( s − 1 )  0, t > 0.5w   y=square(t,d) – формирует периодический сигнал прямоугольной формы (t – вектор значений времени, d – коэффициент заполнения (в процентах), равный отношению длительности импульса к периоду), y=1 в интервале импульса и y=-1 в интервале паузы. y=sawtooth(t,d) – формирует периодический сигнал треугольной формы (t – вектор значений времени, d – длительность нарастания 56 сигнала от -1 до 1 и (1-d) – длительность спадания сигнала от 1 до -1). y=randn(m,n) – формирует матрицу y – дискретный белый шум с нормальным распределением, m – число строк и n – число столбцов. y=wavread(‘имя файла’) – функция чтения данный из файла в формате wav [y,Fs,bit]= wavread(‘имя файла’) – прочитать файл и получить информацию о частоте дискретизации (Fs), количестве бит на отсчет (bit) wavwrite(y,Fs,b,’имя файла’) – записать данные y в формате wav файла, Fs – частота дискретизации данных, b=8 или 16 – число бит на отсчет. 57 Раздел 2. Дискретизация и квантование Обработка сигналов в самом общем случае классифицируется на два больших класса: обработка аналогового сигнала (сигналы – аналитические функции времени, определенные во всех моментах времени по значению) и обработка цифровых сигналов: дискретные функции, как по времени, так и по уровню (амплитуде). Однако можно встретить случаи, когда непрерывный сигнал разделяется на участки (интервалы) по времени (дискретизируется), оставаясь аналоговым по значению (не квантованным). Цифровая обработка сигналов подразумевает дискретизацию и квантование сигнала. Эти две операции являются ключевыми. Цифровой сигнал представляет собой последовательность чисел, которые являются значениями сигнала в определенные моменты времени и называются отсчетами. Как правило, отсчеты берутся через равный интервал времени, называемый периодом дискретизации ts – рис. 2.1. Аналоговый Дискретный Цифровой b Бесконечно много значений a t Бесконечно много значений тts Конечное число значений Конечное число значений тts Конечное число значений Рис. 2.1 – Процессы дискретизации и квантования Очевидно, что при представлении непрерывного сигнала дискретными отсчетами происходит потеря информации, т.к. в моменты времени между отсчетами дискретизированного сигнала о его значениях ничего неизвестно. Кроме того, для представления сигнала в «цифре», т.е. в виде двоичных чисел конечной длины, происходит квантование сигнала по уровню (округления к ближайшему уровню, например вида 2k). Соответственно, процесс преобразование отсчетов в числа называют квантованием по уровню, а округления, которые неизбежно при этом возникают, называются шумом квантования. Квантование рассматривают как случайный процесс, подчиняющийся равномерному распределению (математическое ожидание 58 равно нулю). Дисперсия ошибок округления при равномерном квантовании: 2 = , где 12 = U max 2K – «цена» младшего разряда двоичного слова разрядностью K. На практике число уровней квантования (разрядность представления чисел) выбирают так, чтобы уровень шума квантования был пренебрежимо мал (или шум квантования по уровню был ниже других видов шумов). Поэтому будем рассматривать цифровые сигналы с нулевым шумом квантования, если не оговорено иное. В общем случае, процесс дискретизации может выполняться с переменным периодом дискретизации, т.е. когда «расстояние» между соседними отсчетами не является константой, а изменяется по какому-либо закону (правилу). Например, при сигналах в виде коротких импульсов на фоне медленно изменяющихся процессов выполнять дискретизацию с малым периодом (с высокой частотой) не рационально, т.к. соседние отсчеты в пределе будут фактически отличаться только шумовой составляющей, а при быстро изменяющемся процессе – во время импульса, на длительности которого получится один отсчет или импульс возможно будет пропущен – рис.2.2. Поэтому если известны моменты появления коротких импульсов – можно увеличивать частоту дискретизации на необходимый период времени. Пропущенный импульс Рис. 2.2 – Пример дискретизации составного сигнала Процесс оцифровки (дискретизации и квантования) аналоговых сигналов выполняется с помощью аналогово-цифровых преобразователей (АЦП). В 59 первом приближении дискретизацию сигнала s(t) можно представить в виде идеального ключа – рис. 2.3, который включается и выключается с частотой дискретизации fs. Математически процесс можно представить в виде: 1, при (t − kt s ) = 0 s(k ) = s(t )(t − kts ) , где (t − kt s ) =  0, иначе s(t) s(kts) ключ fs Рис. 2.3 – Идеология процесса дискретизации Например, на рис. 2.4 показан результат равномерной дискретизации синусоидального s(t ) = 0.5 + sin (t + 0.1 ) . сигнала: Поскольку между отсчетами (обозначены красными кружочками) информации о сигнале нет, то обычно соседние отсчеты соединяют (линейная интерполяция) или считают, что за время периода дискретизации сигнал не изменяется (такая интерпретация показана зеленым цветом). s s(t) s(k) 2 4 6 8 10 12 14 16 18 20 Рис. 2.4 – Дискретизация аналогового сигнала и графического представление дискретного сигнала 60 Таким образом, дискретный сигнал s (k ) = 0.5 + sin (kt s + 0.1 ) представляется последовательностью чисел: {0.809016994375; 1.112907053653; 1.344327925502; 1.475916761939; 1.492114701314; 1.391006524188; 1.184547105929; 0.897147890635; 0.562790519529; 0.221008893961; -0.087785252292; -0.327080574275; -0.468583161129; -0.495561964603; -0.404827052466; -0.207106781187; 0.074220708435; 0.405891686681; 0.748689887165; 1.062083377852}. Поскольку равномерную дискретизацию выполнить технически существенно проще, то ее и применяют на практике. В то же время задачу дискретизации сигналов рис. 2.2 выполняют иначе – например, разделяют сигнал на составляющие: медленно меняющуюся и быстро меняющуюся. Соответственно осуществляют дискретизацию каждой составляющей с подходящей частотой дискретизации. Как и процесс дискретизации, квантование также может выполняться равномерно или неравномерно. В случае равномерного квантования сигнала (рис. 2.5), изменяющегося от минимального уровня xmin до максимального xmax, определяют шаг квантования: q= xmax − xmin , где N – число уровней N −1 квантования. xmin+q xmin xmin+2q n xmin+nq xmin+(N-2)q xmin+(n+1)q x(k) xmax Рис. 2.5 – Процесс равномерного квантования сигнала x(k) Тогда каждое значение сигнала принадлежит интервалу n, границы которого пропорциональны величине q, т.е. равны одному из элементов множества чисел xmin+nq, n – целое число из интервала [0, N-1]. Считая, что попадание значения сигнала в интервал квантования n есть случайная величина с равномерным распределением на интервале [-0.5q, 0.5q] можно записать: xкв (n ) = xmin + (n + 0.5)q , ошибка квантования x(n ) − xкв (n ) ≤ 0.5q , q2 математическое ожидание рано нулю, дисперсия (шум квантования): . 12 61 На рис. 2.6 представлена оцифровка синусоидального сигнала и, в частности, образование шума квантования (при представлении отсчетов в пяти разрядной сетке со знаком). амплитуда Аналоговый сигнал sin (ω t) время – t амплитуда Дискретный сигнал sin (ω ts k) дискретное время – индекс k амплитуда Цифровой сигнал sin (ω ts k) дискретное время – индекс k Ошибка квантования (4 бита на амплитуду) 13/16=0.8125 В Значение 0.8 В (0.8·24=12.8) 12/16=0.75 В Рис. 2.6 – Дискретизация и квантование аналогового сигнала 62 В случае неравномерного разбиения интервала значений сигнала на уровни квантования возникает вопрос об оптимальном способе разбиения. Произвольно разобьем интервал значений сигнала [xmin, xmax] на N неравномерных уровней квантования: n=[tn, tn+1). Очевидно t0= xmin и tN-1= xmax. Тогда любое значение отсчета сигнала принадлежит интервалу n, а соответствующее квантованное значение отсчета равно некоторому dn (определяемому внутри n). Поставим вопрос определения границ tn и значений dn таких, чтобы ошибка квантования была минимальна в среднеквадратичном смысле: [ ] N −1t k +1  = M (x − d k ) = ∑ 2 (x − d k )2 f (x )dx , ∫ k =0 tk где f(x) – функция плотности распределения вероятности попадания величины x в интервал tk. Экстремум функции многих переменных достигается в точке, где частные производные по каждой переменной одновременно обращаются в ноль: t j +1  ∂  N −1tk +1 2  ∑ ∫ ( x − d k ) f ( x )dx  = −2 ∫ (x − d j ) f ( x )dx = 0  ∂t j  k =0 tk tj  j = 0,, N − 1 Решение полученной системы уравнений определяет квантователь ЛлойдаМакса (P. Lloyd 1982, Joel Max 1960): dj = t j +1 t j +1 tj tj ∫ xf (x )dx ∫ f (x )dx j = 0,, N − 1 В частности, если f(x) – равномерное распределение: f ( x ) = xmax 1 , то − xmin d j = 0.5(t j +1 + t j ) – равномерный квантователь оптимален в смысле минимума среднеквадратичной ошибки. Неравномерное квантование находит применение при речевых сигналов – применяется экспоненциальное обработке квантование (рекомендации ITU-T G.711). 63 2.1 Дискретизация низкочастотных сигналов Рассмотрим два дискретных сигнала: x1 (kt s ) = sin (t s k ) и x 2 (kt s ) = sin (t s k ) + 0.2 sin (11t s k ) . На рис. 2.7 сигналы представлены: x1 – синяя кривая, x2 – красная кривая. Отметим, что значения отсчетов в моменты времени kt s абсолютно одинаковы (отмечены на рисунке зелеными крестиками). Таким образом, последовательность чисел {x(0), x(1), x(2), ….} не дает однозначного представления аналогового сигнала (можно указать сколь угодно много амплитуда различных сигналов проходящих через наперед заданные точки kt s ). дискретное время – отсчеты k Рис. 2.7 – Неопределенность представления гармонических колебаний набором дискретных отсчетов Рассмотрим причину неоднозначности. Пусть x(k ) = sin (2f 0 kt s ) – искусственно созданный сигнал единичной амплитуды с частотой колебаний f0 и периодом дискретизации ts. Определим одинаковые значения синусоиды (при равных аргументах с учетом периодичности тригонометрической   n   kt s  , функции): sin (2f 0 kt s ) = sin (2f 0 kt s + 2n ) = sin  2 f 0 +  kt s     где k, n – любые целые числа. Очевидно, можно указать кратные значения чисел n и k, тогда m = целое число и тогда sin (2f 0 kt s ) = sin (2( f 0 + mf s )kt s ) , где n – k частота дискретизации f s = t s−1 . Следовательно, синусоиды с частотами f 0 и f 0 + mf s 64 невозможно различить по дискретной последовательности чисел с частотой дискретизации fs (в этом случае частоты со значениями f 0 + mf s при m>0 называют ложными). Это фундаментальное свойство дискретной последовательности, которое приводит к более общему утверждению: поскольку каждая гармоническая составляющая любого сигнала имеет несколько копий, то спектры сложных сигналов (состоящих из множества гармоник) также являются периодичными – с частотой повторения равной частоте дискретизации (в частности, это следует из определения преобразования Фурье). Для снятия неопределенности необходимо чтобы сигнал не содержал гармонических составляющих с частотой более 0.5 f s , эту величину называют частотой Найквиста. Действительно, множество возможных значений m, определяющих величины ложных частот f 0 + mf s , при ограниченной полосе 0.5 f s есть нулевой элемент. Итак, дискретный сигнал имеет периодический спектр. Это приводит к тому, что у каждой гармонической составляющей спектра существует бесконечное множество копий отстоящих от оригинала на величину частоты дискретизации. Рассмотрим периодический сигнал x(t ) , дискретный спектр X ( ) которого ограничен частотой fв. Выполним дискретизацию сигнала x(t ) с частотой дискретизации fs – рис. 2.8. |X(f)| |X(f)| копия копия f -fв fв -0.5 fs -fв 0 fв 0.5fs fs f Рис. 2.8 – Спектры аналогового сигнала x(t) (слева) и дискретного сигнала x(kts) (справа) fs+ fв Далее, постепенно уменьшая частоту дискретизации (при неизменном сигнале и значении fв), можно видеть, что спектры исходного сигнала и его копий 65 сближаются и пересекутся при f s ≤ 2 f в – возникнет наложение спектров (aliasing) – рис. 2.9. |X(f)| копия -fв копия fв fs f Рис. 2.9 – Искажение спектра сигнала при недостаточно высокой частоте дискретизации Наглядным примером наложения может служить иллюзия, когда колесо автомобиля (велосипеда, кареты и т.п.) начинает вращаться против его движения, если между последовательными кадрами оно совершает более чем половины оборота. Пример. Сигнал x(k ) = sin (2f 0 kt s ) + 0.2 sin (8f 0 kt s ) , содержащий первую и четвертую гармоники (f0=1 кГц), дискретизирован с частотой fs. Рассмотрим два случая, определив спектр сигнала и сопоставив его с временной функцией (исходным сигналом), результаты представлены на рис. 2.10 и рис. 2.11. 66 x(k) 1) fs=12 кГц (fs>2fв), наложение отсутствует – сигнал не искажается X(f) Время (индекс k) f, Гц Полоса сигнала Копия спектра Рис. 2.10 – Дискретный сигнал x(k) при fs=12 кГц и его спектр Как и следовало ожидать, спектр периодического сигнала дискретен и содержит ровно столько частотных компонент (в полосе частот 0.5fs) сколько гармоник присутствует в исходном сигнале. 67 x(k) 2) fs=7 кГц (fs<2fв), присутствует искажение спектра Время (индекс k) X(f) Ось симметрии спектра 0.5fs Ложная частота f, Гц Полоса исх. сигнала xиск (k ) = sin (2f 0 kt s ) + 0.2 sin (6f 0 kt s ) x(k) xиск(t) Время (индекс k) Рис. 2.11 – Дискретный сигнал x(k) при fs=7 кГц, его искаженный спектр и сигнал во временной области соответствующий искаженному спектру Искажение спектра из-за эффекта наложения приводит к неправильной трактовке гармонического состава исходного сигнала. 68 Теорема отсчетов Теорема отсчетов (также известна как теорема Котельникова (1933), Найквиста (1928), Шеннона (1949)) устанавливает условия и способ однозначного представления дискретного сигнала. В 1999 году Международный научный фонд Эдуарда Рейна (Германия) признал приоритет В. А. Котельникова. Теорема отсчетов. Если аналоговый сигнал x(t ) имеет финитный (ограниченный по ширине) спектр (верхняя частота сигнала – fв), то он может быть восстановлен однозначно и без потерь по своим дискретным отсчётам, взятым с частотой дискретизации строго большей удвоенной верхней частоты сигнала: f s > 2 fв . Восстановление аналогового сигнала по дискретным отсчетам при 0 < ts < x(t ) = 1 возможно по интерполяционной формуле: 2 fв sin (f (t − kt )) ∑ x(kt s ) f (st − kt )s . ∞ k = −∞ s (*) Выражение s sin ( x ) известно как sinc функция – рис. 2.12. x sinc(x) Главный лепесток Боковые лепестки x Рис. 2.12 – График sinc функции 69 2.2 Дискретное преобразование Фурье Дискретное преобразование Фурье (ДПФ) – распространенная процедура цифровой обработки сигналов. ДПФ позволяет реализовать операции с сигналами, которые недоступны при аналоговой их обработке. ДПФ – это математическая операция, позволяющая определить спектр дискретного сигнала. Дискретный сигнал xд (t ) можно представить суммой произведения отсчетов на сдвинутые во времени дельта-функции (в виде решетчатой функции): ∞ ∞ k = −∞ k = −∞ xд (t ) = ∑ x(kt s ) (t − kt s ) = ∑ x(t ) (t − kt s ) , где x(t ) – непрерывный сигнал (до дискретизации). Спектр такого сигнала (необходимо взять преобразование Фурье): X д ( ) = ∞ ∞ x(kt s ) (t − kt s ) e ∫ k∑ = −∞ −∞ − j t dt = ∞ ∞ k = −∞ −∞ ∑ x(kt s ) ∫ (t − kt s ) e − jt dt . Применяя фильтрующее свойство дельта функции получим: X д ( ) = ∞ ∑ x(kt s ) e − jkt s – спектральная плотность дискретного сигнала. k = −∞ Если входной сигнал периодичен (или мы полагаем, что он периодичен: на периоде сигнала T взято N отсчетов с периодом дискретизации ts), то в соответствии со свойствами преобразования Фурье спектр дискретен и тоже 2 – расстояние между N ts периодичен (разрешение по частоте составляет соседними частотными отсчетами, которые часто называют бинами): 2 k −j nt s  2  N −1  2f s  Nt s   ( ) X k = x nt e = X ± 2f s  . k ∑ s   N   N t s  n =0 Для удобства записи значение ts в аргументах функций опускают, N −1 оставляя только индексы отсчетов: X (k ) = ∑ x(n )e n =0 −j 2 k n N , где k ∈ [0, N − 1] . 70 Величина X (k ) является спектральной плотностью сигнала (является комплексной величиной). Но на практике интересен спектр сигнала, характеризующий амплитудный и/или фазовый состав сигнала. Для этого спектральную плотность нормируют. В частности, придавая физический смысл нулевой гармонике – как постоянной составляющей сигнала получают следующее выражение для ДПФ (которое часто можно встретить в 1 литературе): X (k ) = N N −1 ∑ x(n ) e −j 2 k n N n =0 , k ∈ [0, N − 1] . N −1 Обратное дискретное преобразование Фурье (ОДПФ): x(n ) = ∑ X (k ) e j 2 n k N . k =0 Таким образом, ДПФ – взаимное однозначное преобразование N отсчетов одной области (пространства) в N отсчетов в другой области. Основные свойства ДПФ: 1. Сопряженная симметрия – для любого вещественного сигнала x(n): X (k ) = X ( N − k ) * 2. Линейность: x(n ) = ax1 (n ) + bx2 (n ) ↔ X (k ) = a X 1 (k ) + b X 2 (k ) 3. Циклический сдвиг сигнала на m позиций вправо X = {x0 , x1 , x2 ,, x N −2 , x N −1 } → X′ = {x N −m ,, x N −2 , x N −1 , x0 , x1 ,, x N −m−1 } Y = F {X} Y′ = F {X′} y′k = y k e −j 2 k m N , т.е. спектр циклически сдвинутого сигнала есть исходный спектр с измененными фазами спектральных составляющих. (*) сдвиг влево: y ′ k = y k e j 2 k m N . Следствие: амплитудный спектр сигнала инвариантен к сдвигу. 4. ДПФ свертки Свертка дискретных сигналов A = {a0 , a1 ,, a N −1 }; B = {b0 , b1 ,, bN −1 } N −1 N −1 j =0 j =0 ck = ∑ a j bk − j = ∑ ak − j b j am = 0 и bm = 0 если m < 0 или m ≥ N k = 0,1,,2 N − 1 71 C = A∗B ↔ F {C} = 0.5 N F {A}F {B}: c k = 0 .5 N a k b k Пример расчета спектра. Рассмотрим дискретный периодический сигнал с частотой первой гармоники f0 и периодом дискретизации ts: x(k ) = 2 sin (2f 0 kt s ) + sin (6f 0 kt s + 45°) . Вычислим ДПФ сигнала с длиной выборки N=16 отсчетов при fs=400 Гц для двух случаев: f0=50 Гц и f0=60 Гц (разрешение преобразования Фурье по частоте в этом случае составляет fs =25 Гц, условия теоремы Котельникова N выполняются f s > 2 f в ): 2 n −j k 1 15 1 15 1 15  n   n  X (n ) = ∑ x(k )e 16 = ∑ x(k )cos k  − j ∑ x(k )sin  k  , n ∈ [0, 15] . N k =0 N k =0 N k =0  8   8  Теоретически ожидается, что в результате расчета дискретного спектра в частотной полосе 200 Гц в каждом случае будет две ненулевые составляющие соответствующие первой и третьей гармоникам. 1) f0=50 Гц X (0 ) = 1 15 ∑ x(k ) = 0 16 k =0 X (1) = 1 15 1 15     ( ) x k cos k − j x(k )sin  k  = 0   ∑ ∑ 16 k =0 16 k =0 8  8  X (2 ) = 1 15 1 15     x(k )cos k  − j ∑ x(k )sin  k  = − j1 ∑ 16 k =0 16 k =0 4  4  1 15 1 15  3   3  X (3) = ∑ x(k )cos k  − j ∑ x(k )sin  k  = 0 16 k =0 16 k =0  8   8  X (4 ) = 1 15 1 15     ( ) x k cos k − j x(k )sin  k  = 0   ∑ ∑ 16 k =0 16 k =0 2  2  X (5) = 1 15 1 15  5   5  ( ) x k cos k − j x(k )sin  k  = 0   ∑ ∑ 16 k =0 16 k =0  8   8  X (6 ) = 1 15 1 15  3   3  x(k )cos k  − j ∑ x(k )sin  k  = 0.3535 − j 0.3535 = 0.5∠ − 45° ∑ 16 k =0 16 k =0  4   4  72 X (7 ) = 1 15 1 15  7   7  ( ) x k cos k − j x(k )sin  k  = 0   ∑ ∑ 16 k =0 16 k =0  8   8  X (8) = 1 15 1 15 1 15 ( ) ( ) ( ) ( ) (− 1)k x(k ) = 0 x k cos  k − j x k sin  k = ∑ ∑ ∑ 16 k =0 16 k =0 16 k =0 X (9 ) = 1 15 1 15  9   9  1 15   k x(k )cos k  − j ∑ x(k )sin  k  = ∑ (− 1) x(k )cos k  − ∑ 16 k =0 16 k =0  8   8  16 k =0 8  − j 1 15 (− 1)k x(k )sin  k  = 0 ∑ 16 k =0 8  1 15 1 15  5   5  1 15   k X (10 ) = ∑ x(k )cos k  − j ∑ x(k )sin  k  = ∑ (− 1) x(k )cos k  − 16 k =0 16 k =0  4   4  16 k =0 4  1 15   k − j ∑ (− 1) x(k )sin  k  = 0.3535 + j 0.3535 = 0.5∠45° 16 k =0 4  1 15 1 15  11   11  1 15  3  k ( ) ( ) x k cos k − j x k sin k  = ∑ (− 1) x(k )cos k  −    ∑ ∑ 16 k =0 16 k =0  8   8  16 k =0  8  X (11) = 1 15  3  k − j ∑ (− 1) x(k )sin  k  = 0 16 k =0  8  X (12 ) = 1 15 1 15  3   3  1 15   k ( ) ( ) x k cos k − j x k sin    k  = ∑ (− 1) x(k )cos k  − ∑ ∑ 16 k =0 16 k =0  2   2  16 k =0 2  − j X (13) = 1 15 (− 1)k x(k )sin  k  = 0 ∑ 16 k =0 2  1 15 1 15  13   13  1 15  5  k x(k )cos k  − j ∑ x(k )sin  k  = ∑ (− 1) x(k )cos k  − ∑ 16 k =0 16 k =0  8   8  16 k =0  8  − j 1 15 (− 1)k x(k )sin 5 k  = 0 ∑ 16 k =0  8  1 15 1 15  7   7   1 15  3  k X (14 ) = ∑ x(k )cos k  − j ∑ x(k )sin  k  = ∑ (− 1) x(k )cos k  − 16 k =0 16 k =0  4   4  16 k =0  4  1 15  3  k − j ∑ (− 1) x(k )sin  k  = j1 16 k =0  4  X (15) = 1 15 1 15  15   15  1 15  7  k ( ) ( ) x k cos k − j x k sin k  = ∑ (− 1) x(k )cos k  −    ∑ ∑ 16 k =0 16 k =0  8   8  16 k =0  8  1 15  7  k − j ∑ (− 1) x(k )sin  k  = 0 16 k =0  8  73 x(t) 2π ωt X(n) n 50 Гц 150 Гц 0.5fs Рис. 2.13 – Временное и частотное представление сигнала x(k) при f0=50 Гц Вывод: спектр соответствует ожиданиям и гармоническому составу сигнала. Отметим, плотность. что график Поэтому для показывает определения нормированную амплитуд сигнала спектральную необходимо нормированную спектральную плотность умножить на 2, т.к. каждая гармоника представляется суммой экспонент половинной амплитуды (Am): Am sin ( x ) = Am jx (e − e − jx ) 2 и каждой экспоненте во временной области соответствует дискретная дельта функция в частотной области с половинной амплитудой. Таким образом, амплитуда первой гармоники (частота 50 Гц) имеет величину 2, а амплитуда третьей гармоники (частота 150 Гц) – 1. Гармоники с номерами 10 и 14 – являются симметричными (строго говоря, это «проявление» гармоник с отрицательной частотой). 74 2) f0=60 Гц (расчет производится по тем же формулам, что и в п.1) x(t) 2π ωt X(n) n 50 Гц 150 Гц 0.5fs Рис. 2.14 – Временное и частотное представление сигнала x(k) при f0=60 Гц Вывод: спектр несоответствует ожиданиям, но, если выполнить ОДПФ, то получим исходную временную функцию. Явление, наблюдаемое на рис. 2.14, известно как утечка спектра (размытие спектра). Эффект связан с тем, что частота гармоники (сигнала) не кратна разрешению преобразования Фурье по частоте (в данном примере величина 60 Гц некратна 25 Гц). Поэтому такая гармоника не может отображаться в дискретном спектре, что приводит к аддитивному распределению ее энергии по соседним частотным отсчетам (бинам). Другими словами, некратность частоты сигнала и разрешения ДПФ по частоте приводит к тому, что в выборке сигнала во временной области содержится не целое число периодов сигнала (в данном примере: 2.4 периода). При вычислении ДПФ предполагаются бесконечные пределы интегрирования 75 – для периодического сигнала – аналитическое продолжение. В результате образуется скачек – рис. 2.15, что эквивалентно умножению функции во времени на сигнал прямоугольной формы. Это, в свою очередь, приводит к свертке спектра исходного сигнала и спектра прямоугольного сигнала, что и наблюдается в частотной области в виде эффекта размытия. x(t) исходная выборка сигнала длиной N аналитическое продолжение ωt Рис. 2.15 – Интерпретация конечной выборки сигнала как периодической последовательности 76 2.3 Методы преодоления эффекта утечки спектра Для устранения утечки спектра необходимо, чтобы в выборке сигнала (подлежащего к вычислению спектра) содержалось целое число периодов сигнала, так чтобы при аналитическом продолжении выборки не возникало скачка (разрыва функции). Таким образом, если частота сигнала известна, то возможно выбрать такую частоту дискретизации и/или длину выборки, чтобы обеспечить отсутствие эффекта утечки спектра. Однако на практике, частота сигнала, подлежащего анализу, часто неизвестна (или известна с погрешностью) или требуется анализировать сигналы с разной частотой. В этом случае возникает задача минимизации эффекта. Для этого применяют один из двух методов:  взвешивание сигнала окном;  передискретизация сигнала. 2.3.1 Окна анализа Для предотвращения образования разрыва сигнала при аналитическом продолжении его выборку взвешивают окном анализа (умножают отсчеты во временной области на специально подобранную функцию). Этим добиваются, чтобы на краях выборки отсчеты имели нулевое значение (близкое к нулю). На самом деле, всякий раз, при рассмотрении конечной выборки сигнала получается, что к сигналу применяется прямоугольное окно (все отсчеты сигнала по умолчанию умножаются на 1 и на 0 – за пределами длительности выборки). При этом в случае некратности частоты сигнала и разрешения по частоте, осцилляции спектральных составляющих описываются функцией sinc, которая представляет собой преобразование Фурье от прямоугольной функции. Таким образом, с учетом применения окна анализа w(k), ДПФ для сигнала длительностью N отсчетов имеет вид: 1 X w (m ) = N N −1 ∑ w(k ) x(k ) e −j 2 m k N . k =0 77 В нижеследующей таблице приводятся наиболее распространенные окна анализа. Параметрами окна считают: - ширину главного лепестка; - уровень бокового лепестка по отношению к главному; - положение первого минимума. № Тип окна w(n ) Таблица 2. Окна анализа График (временная и частотная области) 1 Прямоугольное 1, 0 ≤ n < N  0, иначе • уровень бокового лепестка: -13.3 дБ • положение первого минимума: fs/N 2 Треугольное  2n  N , 0 ≤ n < 0 .5 N  2n  2 − , 0.5 N ≤ n < N N  0, иначе   • уровень бокового лепестка: -26.5 дБ • положение первого минимума: 2fs/N 78 3 фон Ханна (Хеннинга)   2n  , 0 ≤ n < N 0.5 − 0.5 cos  N − 1  0, иначе  • уровень бокового лепестка: -31.5 дБ • положение первого минимума: 2fs/N 4 Хемминга   2n  , 0 ≤ n < N 0.54 − 0.46 cos  N − 1  0, иначе  • уровень бокового лепестка: -42.6 дБ • положение первого минимума: 2fs/N 6 Блэкмана   2n   4n  0.42 − 0.5 cos N − 1  + 0.08 cos N − 1 ,      0≤n< N  0, иначе   • уровень бокового лепестка: -58.1 дБ • положение первого минимума: 3fs/N 79 7 Кайзера 2    2 n − N + 1     I0  1 −       N − 1   , 0≤n< N  ( ) I      0, иначе β=7 где I0 – модифицированная функция Бесселя первого рода нулевого порядка, β – параметр, определяющий долю энергии в главном лепестке (чем больше величина - тем шире главный лепесток и меньше уровень боковых лепестков, обычно β=2÷9 • уровень бокового лепестка: -29.9 дБ • положение первого минимума: 1.63fs/N Известны также окна Чебышева, Барлетта, Тьюки, гауссово окно и различные линейные комбинации рассмотренных выше окон (Барлетта-Ханна, БлэкманаХарриса и т.п.). Каждое окно по своим параметрам является компромиссом между шириной главного лепестка (разрешающей способностью по частоте) и уровнем боковых лепестков. Сравнение окон по параметрам осуществляется относительно прямоугольного окна. Прямоугольные окна обеспечивают хорошее разрешение по частоте и полезны, таким образом, для оценки типа гармоник, присутствующих в сигнале. Поскольку затухание прямоугольного окна в частотной области описывается sinc функцией, они вводят некоторое ослабление. Альтернативные функции с меньшим ослаблением (с плоской вершиной и 80 Блэкмана-Харриса) дают максимальную амплитуду, жертвуя разрешением по частоте. Окна Хемминга и фон Ханна (Хеннинга) наиболее приемлемы для общего применения на непрерывных сигналах. В частности, поэтому они по умолчанию применяются во всех цифровых осциллографах. Например, применив к сигналу x(k ) = 2 sin (2f 0 kt s ) + sin (6f 0 kt s + 45°) (с частотой f0=50 и 60 Гц, ts-1=400 Гц) окно Хеннинга можно получить следующую аппроксимацию спектра сигнала (длина выборки N=64) – рис.2.16 и рис. 2.17 – тот же сигнал с прямоугольным окном. Рис. 2.16 – Спектры сигнала с окном Хеннинга при f0=50 Гц (слева) и f0=60 Гц (справа) Рис. 2.17 – Спектры сигнала с прямоугольным окном при f0=50 Гц (слева) и f0=60 Гц (справа) 81 В идеальном случае основной лепесток окна должен быть как можно более узким и плоским, чтобы эффективно дискриминировать все частотные компоненты, а боковые лепестки должны иметь бесконечное ослабление. 2.3.2 Передискретизация сигнала Для преодоления эффекта утечки спектра выполняют передискретизацию сигнала, т.е. изменяют число дискретных отсчетов на периоде сигнала. Этот подход требует знания частоты сигнала, которая априори неизвестна, а измеряется с погрешностью. Поэтому, на практике, удается лишь эффективно минимизировать утечку. В общем случае для реализации передискретизации сигнала (повышения частоты дискретизации или ее понижения) применяется интерполирование (обычно на практике: линейное или квадратичное). В частных случаях можно обойтись дискретными операциями децимации и/или интерполяции. Децимация цифрового сигнала Децимация дискретной последовательности (цифрового сигнала) на величину M есть прореживание последовательности (понижение частоты дискретизации в M раз), при котором из исходного множества отсчетов берется только каждый M-ый отсчет – рис. 2.18. Исходная последовательность, частота дискретизации fs x(0) x(1) x(2) ... x(M-1) x(M) x(M+1) ... x(n) y(0) y(1) y(2) ... x(2M) ↓M ... y(k) Прореженная последовательность, частота дискретизации fs/M Рис. 2.18 – Децимация дискретного сигнала При выполнении децимации следует помнить об уменьшении частоты дискретизации сигнала, что может привести к наложению спектра. Для предотвращения наложения необходимо чтобы выполнялись условия теоремы Котельникова: fs > 2 f в , где fs – частота дискретизации исходного сигнала, fв – M 82 верхняя частота исходного сигнала, M – порядок децимации. Рис. 2.19 наглядно иллюстрирует сближение спектральных компонент при децимации. Случай: fs>2Mfв (M=3) fв fs/3 fв fs/3 fs 2fs/3 3fs/3 Случай: fs<2Mfв (M=4) fв fs/4 fs fв fs/4 fs/2 fs Рис. 2.19 – Спектры сигнала до и после децимации Операция децимации является неинвариантной во времени и не влияет на амплитуду сигнала. Зададим сигнал: x(n ) = {x(0 ), x(1), x(2 ), x(3), x(4 ), x(5), x(6 ),}, выполним прореживание на M=3: y (n ) = {x(0 ), x(3), x(6 ), x(9 ), x(12 ),}. Задержим исходный сигнал (последовательность) на один отсчет: x1 (n ) = {x(1), x(2 ), x(3), x(4 ), x(5), x(6 ),}, выполнив прореживание сигнала x1 на M=3, получим последовательность: y ′(n ) = {x(1), x(4 ), x(7 ), x(10 ),}, которая не является задержанной на один отсчет копией последовательности y (n ) , т.е. ожидаемым множеством отсчетов {x(3), x(6 ), x(9 ), x(12 ),} . Интерполяция цифрового сигнала Интерполирование – это процесс повышения частоты дискретизации. Для увеличения последней в M раз, необходимо между каждой парой отсчетов исходного сигнала x(n ) добавить M-1 отсчет. На x(n) практике в исходную последовательность вводят нули. ↑M y(k) Рассмотрим сигнал: x(n ) = {x(0 ), x(1), x(2 ), x(3), x(4 ), x(5), x(6 ),}. Интерполируем его в M=3 раз, получим новую последовательность: y (k ) = {x(0 ),0,0, x(1),0,0, x(2 ),0,0, x(3),0,0, x(4 ),0,0, x(5),0,0, x(6 ),0,0,}. 83 Соответственно частота дискретизации нового дискретного сигнала y (k ) в 3 раза больше частоты дискретизации (fs) исходного сигнала x(n ) . После введения дополнительных отсчетов и соответственно увеличения частоты дискретизации для получения правильного спектра необходимо удалить ненужные копии спектра исходного сигнала, которые называются фантомами или изображениями – рис.2.20. Спектр исходного сигнала с частотой дискретизации fs fв fs 2fs 3fs фантомы Спектр интерполированного сигнала (М=3), частота дискретизации Fs=3fs fв fs 2fs 3fs=Fs Спектр исходного сигнала после интерполяции (М=3), частота дискретизации Fs=3fs fв fs 2fs Fs Рис. 2.20 – Спектры сигнала до и после интерполяции Пример. Изменить частоту дискретизации сигнала в 1.5 раза. x(t) ↑3 y(t) ↓2 v1 x v2 fв fs fв fs fв fs 1.5fs 3fs fв fs 1.5fs 3fs v1 ФНЧ с линейной фазой 2fs 3fs v2 АЧХ ФНЧ fв fs y Рис. 2.21 – Процесс изменения частоты дискретизации сигнала в 1.5 раза 84 Раздел 3. Системы обработки сигналов В самом общем случае под системой обработки сигналов понимается какой-либо преобразователь входных сигналов в выходные сигналы. В простейшем случае система представляет собой устройство (физическое или виртуальное – в виде математической модели) и содержит один входной сигнал и один выходной сигнал. Классическим примером такой системы служат аналоговые частотно-избирательные фильтры. Системы делят на линейные и нелинейные. В линейных системах, в отличие от нелинейных, справедлив принцип суперпозиции (реакция системы на сумму сигналов эквивалентна сумме реакций системы на каждый входной сигнал в отдельности). Система стационарна, если произвольная задержка входного сигнала приводит к такой же задержке выходного сигнала (система с постоянными коэффициентами), Стационарность в противном системы случае означает – система независимость нестационарна. от времени и инвариантность к сдвигу во времени. Фактически все системы цифровой обработки сигналов являются линейными инвариантными во времени системами (сокращенное название ЛИС или ЛИВ). Дискретные ЛИС системы описываются:  во временной области – импульсной характеристикой (реакция системы на единичный импульс);  в частотной области – функцией передачи (передаточной функцией). Рассмотрим простейшую ЛИС систему, на вход которой подан входной дискретный (или цифровой) сигнал x(n ) и на выходе наблюдается выходной дискретный (или цифровой) сигнал y (n ) . x(n) ЛИС система y(n) 85 Тогда, зная импульсную характеристику системы h(n ) , можно связать выходной сигнал с входным посредством разностного уравнения: ∑ an y(k − n ) = ∑ bm x(k − m ) , n m где коэффициенты an, bm полностью характеризуют импульсную характеристику системы h(n ) и определяют ее передаточную функцию. В общем случае (теоретически) пределы индексов могут быть и отрицательными. Это означает, что текущий отсчет системы зависит не только от предыстории, но и от будущего (т.е. требует знания отсчетов, которые будут получены). Поэтому, если отклик системы равен нулю h(k ) = 0 при k<0 (отклик системы зависит только от предыстории), то систему называют каузальной (причинной). В такой системе реакция на входной сигнал появляется только после поступления сигнала на ее вход. Некаузальные системы физически невозможно реализовать в реальном масштабе времени. Если требуется реализовать свертку сигналов с двусторонними операторами (например, при дифференцировании, преобразовании Гильберта, и т.п.), то это выполняется задержкой (сдвигом) входного сигнала (обычно на количество отрицательных индексов в выражении свертки). Пример каузальной ЛИС системы. Пусть имеются отсчеты импульсной характеристики {h(0), h(1), h(2)}, тогда y (n ) = x (n )h(0 ) + x (n − 1)h(1) + x (n − 2 )h(2 ) x(n) задержка задержка x(n-1) x(n-2) h(2) h(1) Σ y(n) h(0) Рис. 3.1 – Блок схема ЛИС системы 86 Таким образом, зная импульсную характеристику системы h(k ) (как множество отсчетов длины K) всегда можно рассчитать отклик y (n ) на любую K −1 входную последовательность x(n ) : y (n ) = ∑ x(n − k )h(k ) . k =0 Как известно, свертке во временной области соответствует произведение Фурье образов функций в частотной области. Если известно h(n ) , то возможно рассчитать (с помощью ДПФ) спектр импульсной характеристики H (k ) – передаточная функция системы. Аналогично для входного сигнала x(n ) рассчитаем спектр X (k ) . Тогда спектр выходного сигнала y (n ) равен: Y (k ) = H (k ) X (k ) . Выполнив ОДПФ последовательности Y (k ) можно получить временную зависимость y (n ) . Это один из путей практической реализации цифровых систем. 87 3.1 Z-преобразование В общем случае, дискретные системы, как и аналоговые, описываются с помощью дифференциальных уравнений. Отличия состоят в представлении решения. В непрерывных функциях применялось преобразование Лапласа F (s ) = ∫ f (t )e − st dt , позволяющее свести дифференциальное уравнение к алгебраическому. В случае дискретных функций (сигналов) применяется zпреобразование, применяемое к разностному уравнению. Смысл z-преобразования состоит в том, что последовательности h(k ) ставится в соответствие функция комплексной переменной, определяемая следующим образом: H ( z ) = ∞ ∑ h(k ) z −k , k = −∞ функция H ( z ) определена для тех z, при которых ряд сходится. 3.1.1 Z-преобразование основных сигналов 1. Единичная импульсная функция 1, n = 0 x(n ) =  0, n ≠ 0 X (z ) = ∞ ∑ x(k )z −k = z 0 = 1 , k = −∞ функция определена на всей комплексной плоскости 2. Единичный скачок 1, n ≥ 0 x(n ) =  0, n < 0 X (z ) = ∞ ∑ x(k )z −k k = −∞ ∞ = ∑ z −k , k =0 ряд (геометрическая прогрессия) сходится при z > 1 : X ( z ) = 1 1 − z −1 3. Степенная функция a n , n ≥ 0 x(n ) =  0, n < 0 88 X (z ) = ∑ x(k )z −k = ∑ a k z −k = ∑ (a −1 z ) , ∞ ∞ ∞ k = −∞ k =0 k =0 −k ряд (геометрическая прогрессия) сходится при z > a : X ( z ) = 1 1 − az −1 4. Затухающая синусоида x(n ) = a n cos(nt s +  ) e j (kts + ) + e − j (kts + ) k −k X ( z ) = ∑ x(k )z = ∑ a z = 2 k = −∞ k =0 k e j ∞ e − j ∞ j t s k − k ( ( = ae ) z + ae − jts ) z −k ∑ ∑ 2 k =0 2 k =0 ∞ −k ∞ cos( ) − z −1a cos(t s −  ) X (z ) = 1 − z −1 2a cos(t s ) + a 2 z −2 3.1.2 Связь z-преобразования с Фурье преобразованием Дискретное z-преобразование связано с преобразованием Лапласа и соответственно с преобразованием Фурье. ∞ Рассмотрим сигнал x(t ) = ∑ x(k ) (t − kt s ) , применим преобразование k =0 ∞ ∞ ∞ k =0 Лапласа: X (s ) = ∫ x(t )e − st dt = ∑ x(k )∫  (t − kt s )e − st dt . Далее, используя фильтрующее свойство дельта функции, получим: ∞ X (s ) = ∑ x(k ) e − skts . k =0 ∞ Делая подстановку z = e sts , запишем выражение: X ( z ) = ∑ x(k ) z −k . k =0 Таким образом, связь между преобразованием Лапласа и z- преобразованием выражается следующим образом:   1 X ( z ) = X  s = ln z  и X (s ) = X z = e sts . ts   ( ) Аналогично для Фурье преобразования:  1  X ( z ) = X  ln z  и X ( j ) = X (e jts ).  jt s  89 3.1.3 Свойства z-преобразования Пусть сигналы x1 (n ) , x2 (n ) имеют z-преобразования X 1 ( z ) и X 2 ( z ) соответственно. 1. Линейность Z-преобразование линейной комбинации сигналов ax1 (n ) + bx2 (n ) есть линейная комбинация их z-преобразований aX 1 ( z ) + bX 2 ( z ) , где a, b ∈ ℜ 2. Задержка Z-преобразование последовательности x1 (n ) задержанной на n0 отсчетов: y (n ) = x1 (n − n0 ) есть Y (z ) = ∞ ∑ x1 (n − n0 )z −n = z −n0 n = −∞ ∞ ∑ x1 (n − n0 )z −(n−n ) = z −n n = −∞ X 1 (z ) Множитель z − n0 называется оператором задержки на n0 отсчетов (тактов). 3. Свертка Z-преобразование выражения свертки сигналов y (n ) = ∞ ∑ x1 (k )x2 (n − k ) есть k = −∞ произведение их z-преобразований: Y (z ) = =  ∞ n = −∞  k = −∞   ∞  ∞ n = −∞  k = −∞   ∞ ∑  ∑ x1 (k )x2 (n − k ) z −n = ∑  ∑ x1 (k )x2 (n − k )z −(n−k ) z −k  = ∞ ∞ ∑ x1 (k )z ∑ x2 (n − k )z −(n−k ) = X 1 (z )X 2 (z ) k = −∞ −k n = −∞ 3.1.4 Обратное z-преобразование Соответствие между дискретной последовательностью (цифровым сигналом) и ее z-преобразованием является взаимно-однозначным. Поэтому переход из z-области к последовательности во временной области осуществляется по формуле: x(k ) = 1 X ( z )z k −1dz ∫ j 2 где интегрирование ведется по контуру в области сходимости функции X ( z ) , охватывающему все ее полюса. 90 На практике обратное z-преобразование находят разложением функции X ( z ) на простые дроби и последующим сопоставлением в соответствие каждому слагаемому (простой дроби) одно из типовых разложений (по аналогии с нахождением оригинала функции в операторном методе). 91 3.2 Описание цифровой системы в z-области ЛИС система цифровой обработки сигналов аналогично системе аналоговой обработки сигналов характеризуются во временной области: импульсной характеристикой h(k ) и в частотной области – передаточной функцией H ( z ) . Причем связь между ними в соответствии с определением zпреобразования осуществляется следующим образом: H (z ) = Y (z ) ∞ = ∑ h(k ) z −k , X ( z ) k =0 Система {h(k)} x(k) y(k) где Y ( z ), X ( z ) – z-преобразования выходного и входного сигналов. В общем случае, система взвешивает отсчеты входного сигнала x(k ) и отсчеты выходного сигнала y (k ) во временной области, реализуя, таким образом, дробно-рациональную функцию H ( z ) , определяющую частотные характеристики (делая подстановку z = e jts и разлагая комплексную функцию на модуль и аргумент, получаем АЧХ и ФЧХ системы). Поэтому во временной области цифровые системы описываются разностным уравнением (с вещественными коэффициентами bm и an): ∞ ∞ m =0 n =1 y (k ) = ∑ bm x(k − m ) − ∑ an y (k − n ) . Применив к разностному уравнению z-преобразование, получим ∞ ∞  −n  Y ( z )1 + ∑ an z  = X ( z )∑ bm z −m . m =0  n=1  ∞ Следовательно H ( z ) = ∑ bm z −m m =0 ∞ 1 + ∑ an z – передаточная функция в z-области. −n n =1 Характеристики системы зависят от количества и величин коэффициентов bm и an. Разложив числитель и знаменатель функции передачи на множители получим выражение: 92 H (z ) = K 0 ∏ (1 − z m z −1 ) m ∏ (1 − pn z ) −1 , n где zm – нули функции, pn – полюса функции, K0 – коэффициент усиления. Как и в аналоговых системах полюса могут быть комплексно-сопряженными или действительными, в том числе кратными. 3.2.1 Устойчивость дискретных систем Система называется устойчивой, если при любых начальных условиях свободные колебания (собственные) являются затухающими, т.е. lim y (k ) = 0 k →∞ при x(k ) = 0 . Поскольку любой сигнал на выходе ЛИС системы есть линейная комбинация задержанных во времени импульсных характеристик, то для обеспечения устойчивости необходимо, чтобы импульсная характеристика системы была затухающей. Как и в случае дифференциальных уравнений, так и в случае разностных уравнений, описывающих дискретную систему, – их решения могут иметь несколько видов (или их сочетаний) в зависимости от корней характеристического уравнения (полюсов передаточной функции): - действительные отрицательные корни; - комплексно-сопряженные корни; - кратные корни. Выражение H ( z ) можно разложить на простые дроби: N Rn H ( z ) = ∑∑ n =1 r =1 qnr (1 − p z ) −1 r n + M ∑ k m z −m , m =0 где pn – полюса функции H ( z ) (корни знаменателя), qnr – соответствующие им вычеты, Rn – кратность корня, коэффициенты km характеризуют целую часть функции H ( z ) . 93 Тогда в импульсной характеристике h(k ) слагаемым вида Rn ∑ r =1 qnr (1 − p z ) −1 r n будут соответствовать слагаемые типа: qn pnk – для отрицательных некратных корней; Rn ∑ Ai−1k i−1 pnk – для кратных корней; i =1 qn pnk + qn* ( pn* ) = 2 qn pn cos( n k + arg(qn )) – для комплексно-сопряженных k k корней ( pn =  n + j n ). Следовательно, для затухания импульсной характеристики необходимо чтобы модули полюсов системы были меньше единицы: pn < 1 . Другими словами: дискретная система устойчива, если ее полюса лежат на комплексной плоскости внутри круга единичного радиуса. 3.2.2 Фазовая и групповая задержка Задержка распространения сигналов во времени относится к характерной особенности каузальных систем. Фазовая задержка – это прямая характеристика временной задержки системой гармонических колебаний. При подаче на вход системы гармоники sin (t ) , сигнал на выходе каузальной системы, без учета изменения его амплитуды, равен sin (t −  ) , при этом: sin ( t −  ) = sin [ (t − t d )] , т.е. задержка td на частоте ω равна: t d = Распространяя результат на каждую частоту: t d () =  .  ( ) .  Постоянство значения t d ( ) в определенном частотном диапазоне обеспечивает для всех гармоник сигнала такое же соотношение их фазовых характеристик, какое было на входе системы, т.е. не изменяет формы сигнала (при постоянстве АЧХ в этом же частотном диапазоне). Каузальные системы, как правило, имеют в рабочем диапазоне определенную зависимость значения td от частоты, которая характеризуется групповым временем запаздывания (ГВЗ). ГВЗ характеризует среднюю 94 временную задержку составного сигнала (содержащего набор гармонических колебаний). Допустим, что сигнал на входе системы представляет собой сумму двух гармоник с близкими частотами: x(t ) = cos(1t ) + cos( 2 t ) , что тождественно тригонометрической записи: x(t ) = 2 cos[0.5(1 +  2 )t ]cos[0.5(1 −  2 )t ] . Это означает, что сумму двух гармоник с частотами ω1 и ω2 можно рассматривать, как амплитудную модуляцию гармоники с частотой 0.5(ω1+ω2) гармоникой с частотой 0.5(ω1-ω2). При прохождении через систему каждая из гармоник ω1 и ω2 может получить различную задержку, при этом сигнал на выходе фильтра, без учета амплитудных изменений: y (t ) = cos(1t −  1 ) + cos( 2 t −  2 ) , что тождественно: y (t ) = 2 cos[0.5(1 +  2 )t − 0.5( 1 +  2 )] ⋅ cos[0.5(1 −  2 )t − 0.5( 1 −  2 )] . Пульсацию колебаний (определяется меньшей частотой) выразим через групповую временную задержку tg: cos[0.5(1 −  2 )t − 0.5(1 −  2 )] = cos[0.5(1 −  2 )(t − t g )] . Отсюда: (1 −  2 )t g = (1 −  2 ) и t g =  1 −  2  – наклон фазочастотной = 1 −  2  характеристики системы. При распространении выражения на непрерывную частотную характеристику получим: t g ( ) = d . d Таким образом, для того чтобы система обладала линейной фазовой характеристикой необходимо и достаточно обеспечить постоянную величину ГВЗ. Это выполняется, если импульсная характеристика системы h(n ) симметрична: h(n ) = h( N − n − 1) или антисимметрична: h(n ) = −h( N − n − 1) , n=0, 1, …, (N-1)/2, если N – нечетное; n=0, 1, …, (N/2)-1, если N – четное. (*) Примем этот факт без доказательства. 95 3.3 Цифровые фильтры Цифровая фильтрация является ключевой операцией ЦОС. Как и аналоговые фильтры, цифровые фильтры предназначены для избирательного изменения спектрального состава сигнала (задачи подавления помех, сглаживания, частотного разделения, выделения сигналов и т.п.). Фильтр является классическим примером ЛИС системы. Поэтому характеристики фильтров характеризуются его импульсной характеристикой, которая полностью определяет его частотные свойства (АЧХ и ФЧХ). При этом порядок фильтра определяется числом коэффициентов в разностном ∞ ∞ m =0 n =1 уравнении: y (k ) = ∑ bm x(k − m ) − ∑ a n y (k − n ) . На рис. 3.2 представлена классификация фильтров по основным признакам. Фильтры с конечной импульсной характеристикой (КИХ фильтры) – нерекурсивные фильтры, коэффициенты an=0 Частотной селекции: ФНЧ, ФВЧ, ППФ и т.п. Фильтры с бесконечной импульсной характеристикой (БИХ фильтры) – рекурсивные фильтры Эвристические - медианный - авторские Цифровые фильтры Оптимальные (квазиоптимальные) Адаптивные На основе МНК Линейные: - фильтры Калмана - максимального ОСШ - фильтры Винера Нелинейные - во временной области - в частотной области - с обратной связью Рекурсивный МНК - прямые - решетчатые Робастные (с коррекцией импульсной или частотной характеристик) На основе нелинейных методов Рис. 3.2 – Обобщенная классификация цифровых фильтров 96 К основным параметрам фильтров относят: - порядок фильтра (число коэффициентов фильтра); - полоса пропускания; - полоса задерживания; - неравномерность АЧХ в полосе пропускания (δ1); - уровень пульсаций в полосе задерживания (δ2); - групповое время запаздывания. A переходная полоса δ1 полоса пропускания полоса задерживания δ2 f Рис. 3.3 – Обобщенный вид АЧХ фильтра Цифровые фильтры (ЦФ) по сравнению с аналоговыми обладают рядом преимуществ: • ЦФ могут иметь параметры, реализация которых невозможна в аналоговых фильтрах (например, линейную фазовую характеристику); • Работоспособность ЦФ не зависит от дестабилизирующих факторов внешней среды (например, температуры); • Один ЦФ может обрабатывать несколько входных каналов или сигналов; • Входные и выходные данные можно сохранять для последующего использования; • Точность цифровых фильтров ограничена только разрядностью отсчетов; • Фильтры могут использоваться при очень низких частотах и в большом диапазоне частот, для чего достаточно только изменять частоту дискретизации данных. 97 3.3.1 КИХ-фильтры Нерекурсивные фильтры описываются импульсной характеристикой: M y (k ) = ∑ bm x(k − m ) m =0 Импульсная характеристика КИХ фильтров имеет конечную длину (M+1 отсчет), при этом количество используемых предыдущих отсчетов M называют порядком фильтра. Набор коэффициентов bm представляет собой импульсную характеристику h(k ) . x(n) z-1 b0 x(n-1) x(n-i) z-1 … b1 z-1 … z-1 bi x(n-M) bM … … Σ y(n) Рис. 3.4 – Прямая форма реализации КИХ фильтра Вследствие отсутствия обратных связей (текущий отсчет y (k ) зависит только от набора входных отсчетов), КИХ фильтры всегда устойчивы. Это обстоятельство и простота их реализации служит широкому распространению КИХ фильтров. Однако реализация на их основе фильтров с хорошей прямоугольностью АЧХ требует высокого порядка фильтров (M – до нескольких сотен, тысяч). Несомненным достоинством КИХ фильтров является возможность реализации абсолютно линейной ФЧХ и, следовательно, постоянной величины ГВЗ, что обеспечивает условия отсутствия искажений формы сигнала. Для обеспечения линейности фазы необходимо чтобы импульсная характеристика была симметричной bm = bM −m или антисимметричной bm = −bM −m . ГВЗ таких фильтров составляет 0.5M отсчетов (или 0.5Mts секунд, ts – период дискретизации). 98 КИХ фильтры могут быть и некаузальными: y (k ) = M ∑ bm x(k − m ) , т.е. m=− M ′ реализуются сдвигом входных данных на величину M` отсчетов. При M = M ′ фильтр называется двусторонним симметричным. Симметричные фильтры, в отличие от односторонних фильтров, не изменяют фазы обрабатываемого сигнала. Интервал суммирования по m получил название "окна" фильтра. Пример 1. Пусть дана импульсная характеристика фильтра h(k ) = {0.25,0.5,0.25}. y (k ) = 0.25 x(k ) + 0.5 x(k − 1) + 0.25 x(k − 2 ) Определим частотные характеристики фильтра. H ( z ) = 0.25 z 0 + 0.5 z −1 + 0.25 z −2 H ( j ) = 0.25 + 0.5e − jts + 0.25e − j 2 ts = 0.25 + 0.5 cos(t s ) − j 0.5 sin (t s ) + + 0.25 cos(2t s ) − j 0.25 sin (2t s ) АЧХ: H ( ) = (0.25 + 0.5 cos(t s ) + 0.25 cos(2t s ))2 + (0.5 sin (t s ) + 0.25 sin (2t s ))2 ФЧХ: ( ) = − arctg 0.5 sin (t s ) + 0.25 sin (2t s ) 0.25 + 0.5 cos(t s ) + 0.25 cos(2t s ) Рис. 3.5 – АЧХ и ФЧХ фильтра 2-го порядка 99 Пример 2. Пусть дана импульсная характеристика фильтра h(k ) = {0.25,0.25,0.25,0.25} – усредняющий фильтр: y (k ) = 1 [x(k ) + x(k − 1) + x(k − 2) + x(k − 3)]. 4 Определим частотные характеристики фильтра: H ( z ) = 0.25(1 + z −1 + z −2 + z −3 ) . H ( j ) = H ( ) = 1 [1 + cos(t s ) + cos(2t s ) + cos(3t s ) − j (sin (t s ) + sin (2t s ) + sin (3t s ))] 4 (1 + cos(t s ) + cos(2t s ) + cos(3t s ))2 + (sin (t s ) + sin (2t s ) + sin (3t s ))2 sin (t s ) + sin (2t s ) + sin (3t s ) ( ) = − arctg 1 + cos(t s ) + cos(2t s ) + cos(3t s ) 1 2 Рис. 3.6 – АЧХ и ФЧХ усредняющего фильтра 3-го порядка 100 3.3.2 БИХ-фильтры Рекурсивные фильтры описываются импульсной характеристикой: M N m =0 n =1 y (k ) = ∑ bm x(k − m ) − ∑ an y (k − n ) Наличие обратной связи (за счет слагаемых с коэффициентами an) позволяет получить фильтры с бесконечной импульсной характеристикой – рис. 3.7. По этой же причине БИХ фильтры могут быть неустойчивыми. x(n) b0 y(n) Σ z-1 z-1 b1 -a1 z-1 z-1 b2 … -a2 … … bM-1 … -aN-1 z-1 z-1 bM -aN Рис. 3.7 – Прямая форма реализации БИХ фильтра Мысленно разделим слагаемые на две части: рекурсивную и нерекурсивную, а также объединим слагаемые относительно задержек – получим еще одну форму реализации БИХ фильтра – каноническую (или прямая форма 2) – рис. 3.8. Каноническая форма требует минимального объема памяти для запоминания отсчетов. Если в прямой реализации поменять местами блоки умножения и задержки, а также сгруппировать слагаемые относительно задержек можно получить транспонированную реализацию БИХ фильтра – рис. 3.9. Транспонированная форма позволяет распараллеливать процесс вычислений. 101 x(n) Σ b0 Σ y(n) z-1 b1 -a1 … … -aM … bM … … z-1 -aN Рис. 3.8 – Каноническая форма реализации БИХ фильтра Наконец, из представления передаточной функции в виде разложения на множители: H ( z ) = K 0 ∏ (1 − z m z −1 ) m ∏ (1 − pn z −1 ) следует возможность реализации БИХ n фильтров в виде каскадного соединения фильтров 1-го или 2-го порядков. Эта форма часто применяется на практике, т.к. позволяет ослабить влияние эффектов квантования (округлений). Параллельную форму можно получить, воспользовавшись разложением H ( z ) на простые дроби. Эта форма также содержит звенья 1-го и 2-го порядков и допускает параллельные вычисления. 102 x(n) b0 y(n) + z-1 b1 + -a1 z-1 b2 + -a2 … … … bM + -aM … … z-1 -aN Рис. 3.9 – Транспонированная форма реализации БИХ фильтра Пример 1. Пусть фильтр описывается уравнением y (k ) = b0 x(k ) − a1 y (k − 1) . Определим его частотные характеристики: H ( z ) = b0 , очевидно для 1 + a1 z −1 устойчивости фильтра необходимо чтобы a1 < 1 . H ( j ) = b0 b0 = 1 + a1e − jts 1 + a1 cos(t s ) − ja1 sin (t s ) АЧХ: H ( ) = b0 (1 + a1 cos(t s ))2 + (a1 sin (t s ))2 ФЧХ: ( ) = − arctg − a1 sin (t s ) 1 + a1 cos(t s ) Зададим коэффициенты фильтра: b0=0.25, a1=-0.75. 103 Рис. 3.10 – АЧХ и ФЧХ фильтра 1-го порядка x(n) y(n) 0.25 + 0.25 + x(n) y(n) z-1 z-1 0.75 0.75 Рис. 3.11 – Прямая (слева) и каноническая (справа) формы реализации БИХ фильтра Пример 2. y (k ) = b0 x(k ) + b1 x(k − 1) − a1 y (k − 1) − a2 y (k − 2 ) b0 + b1 z −1 H (z ) = 1 + a1 z −1 + a2 z −2 H ( j ) = b0 + b1 cos(t s ) − jb1 sin (t s ) 1 + a1 cos(t s ) − ja1 sin (t s ) + a2 cos(2t s ) − ja 2 sin (2t s ) Зададим коэффициенты фильтра: b0=0.2, b1=0.2, a1=-0.96, a2=0.36 и рассчитаем частотные характеристики. Определим полюса функции H ( z ) : z1 = 0.48 + j 0.36 , z 2 = 0.48 − j 0.36 и ее нули: z 0 = −1 . Так как полюса комплексные, то реализовать фильтр второго порядка двумя фильтрами первого порядка с действительными коэффициентами нельзя. 104 f/fs f/fs Рис. 3.12 – АЧХ (слева) и ФЧХ (справа) БИХ фильтра 2-го порядка Пример 3. Дана передаточная функция: 3 z −2 + 1 . Записать разностное H ( z ) = −2 2 z − 3 z −1 + 1 уравнение и сделать выводы об устойчивости фильтра. Разностное уравнение: y (k ) = x(k ) + 3 x(k − 2 ) − [− 3 y (k − 1) + 2 y (k − 2 )] . Корни знаменателя H(z): 2 z −2 − 3 z −1 + 1 = 0 , откуда z1=1, z2=2. Так как корни лежат вне круга единичного радиуса, то фильтр неустойчив. 3.3.3 Соединение фильтров Последовательное предшествующего соединение фильтра фильтров: является входным выходной для сигнал последующего. Передаточная функция общей системы равна произведению передаточных функций фильтров в нее входящих: H ( z ) = ∏ H n ( z ) . n x(k) H1(z) H2(z) H3(z) … HN(z) y(k) H(z) Рис. 3.13 – Последовательное соединение фильтров 105 Параллельное соединение фильтров: сигнал подается на входы всех параллельно соединенных фильтров одновременно, выходные сигналы фильтров суммируются. Передаточная функция общей системы равна сумме передаточных функций фильтров в нее входящих: H ( z ) = ∑ H n ( z ) . n x(k) H1(z) H2(z) H3(z) … HN(z) y(k) + H(z) Рис. 3.14 – Параллельное соединение фильтров Соединение обратной связи. Сигнал первого фильтра подается на выход системы и одновременно на вход фильтра обратной связи, выходной сигнал которого суммируется (со знаком плюс или минус в зависимости от вида связи: отрицательной или положительной), с входным сигналом системы. Эквивалентная передаточная функция системы: H ( z ) = H1 (z ) . 1 ± H 1 ( z )H 2 ( z ) x(k) + H1(z) y(k) H2(z) H(z) Рис. 3.15 – Соединение с обратной связью 106 3.4 Обобщенное описание дискретной свертки Выражение линейной свертки y (k ) = M −1 ∑ h(m )x(k − m ) применялось m =0 неоднократно и кажется простым. Однако, рассмотрим дискретный сигнал x(n ) (N отсчетов) и импульсную характеристику h(m ) (M отсчетов), тогда результат фильтрации сигнала, строго говоря, будет содержать K = N + M − 1 отсчетов. Рассмотрим пример: h(m ) = {0.25,0.5,0.25} (M=3, фильтр 2-го порядка), x(m ) = {1,1,1,1,1} (N=5) при 0≤m≤4 и x(m ) = 0 при других значениях m. Тогда свертка определяется следующим образом (7 отсчетов): y (0 ) = h(0 )x(0 ) + h(1)x(− 1) + h(2 )x(− 2 ) = 0.25 y (1) = h(0 )x(1) + h(1)x(0 ) + h(2 )x(− 1) = 0.75 y (2 ) = h(0 )x(2 ) + h(1)x(1) + h(2 )x(0 ) = 1.0 y (3) = h(0 )x(3) + h(1)x(2 ) + h(2 )x(1) = 1.0 y (4 ) = h(0 )x(4 ) + h(1)x(3) + h(2 )x(2 ) = 1.0 y (5) = h(0 )x(5) + h(1)x(4 ) + h(2 )x(3) = 0.75 y (6 ) = h(0 )x(6 ) + h(1)x(5) + h(2 )x(4 ) = 0.25 k=0 x(5) x(4) 1 x(3) 1 x(2) 1 x(1) 1 k=1 x(4) x(3) x(2) k=2 x(4) x(3) h(0) 0.25 h(1) 0.5 h(2) 0.25 x(0) 1 x(-1) x(-2) h(0) h(1) h(2) x(1) x(0) h(0) h(1) h(2) x(2) x(1) x(0) y (0 ) = x(-3) x(-5) x(-6) M −1 ∑ h(m )x(1 − m ) m =0 y (2 ) = ∑ h(m )x(− m ) m =0 x(-4) y (1) = M −1 M −1 ∑ h(m )x(2 − m ) m =0 107 h(0) k=3 x(4) k=4 k=5 k=6 h(1) h(2) x(3) x(2) x(1) h(0) h(1) h(2) x(4) x(3) x(2) h(0) h(1) h(2) x(4) x(3) h(0) h(1) h(2) x(4) y (3) = x(0) M −1 ∑ h(m )x(4 − m ) m =0 M −1 ∑ h(m )x(5 − m ) m =0 x(1) y (6 ) = x(3) x(0) y (5) = x(2) ∑ h(m )x(3 − m ) m =0 y (4 ) = x(1) M −1 x(0) M −1 ∑ h(m )x(6 − m ) m =0 x(2) x(1) x(0) Рассматривая процесс вычисления свертки и обращая внимание на индексы, можно отметить: − в общем случае, при свертке двух последовательностей неважно какую мы называем сигналом и какую импульсной характеристикой, поэтому неважно какую последовательность дополним нулями; − отрицательные индексы соответствуют «будущему», всегда можно переиндексировать последовательность дополненную нулями так чтобы не было отрицательных индексов: x′(m ) = {0,0,1,1,1,1,1} (в данном примере на M-1 нулевой отсчет); − суммирование произведений соответствует умножению вектора на строку матрицы. Поэтому операцию линейной свертки можно представить в виде произведения матриц: 108 0 h(2 ) h(1)   x(0 )  y (0 )  h(0 ) 0      0 h(2 )  x(1)   y (1)   h(1) h(0 ) 0       y (2 ) h(2 ) h(1) h(0 ) 0 0   x(2 )       y (3) =  0 h(2 ) h(1) h(0 ) 0 0   x(3)       y (4 )  0 0 h(2 ) h(1) h(0 ) 0 0   x(4 )       y (5)  0   0 h(2 ) h(1) h(0 ) 0 0        y (6 )  0   0 h(2 ) h(1) h(0 )  0     Или в общем случае: Y = TX (для определенности положим N>M) Y = [ y (0 ),..., y ( N + M − 1)] , X = [x(0 ),..., x( N − 1),0,,0] T Квадратная матрица T T является теплицевой матрицей размерностью K = N + M − 1.  h0  h h0  1  h2 h1  ...  ... hM −2 hM −3   hM −1 hM −2  0 hM −1 T=  0  ... ...   0  0   ... ...   0  0 h0 ... hM −4 hM −3 hM −2 hM −1 ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... h0 h1 h2 h3 ... ... h0 h1 h2 ... ... h0 h1 ... ... h0 ... ... 0 ... 0 ... 0 ... ... ... ... ... ... ... ... ... ... ... h0 ... h1 h0 ... ... ... ... ... ... ... 0 0 0 0 ... hM −1 hM −2 0 0 0 0 ... 0 hM −1 ... ... ... ... ... ... ... ... ... ... ... ... ... ... h2 h1  0 h3 h2  1  h4 h3  2  ... ...  … 0 hM −1  M-2  0  M-1 0  M  0  M+1 ... ...  …  0  N-1 0  N  ... ...  …  h0 0  M+N-2 h1 h0  M+N-1 109 3.5 Циклическая свертка Рассмотрим произведение двух дискретных сигналов {x1 (n )}, {x2 (n )} одинаковой длины N: y (n ) = x1 (n )x2 (n ) . Определим спектр этого произведения: Y (k ) = 1 N N −1 ∑ X 1 (n )X 2 (k − n ) , где X1 и X2 – спектры сигналов n =0 x1 и x2 . При вычислениях по этой формуле потребуются отсчеты с номерами, [0, N − 1] . выходящими за пределы интервала В этом случае учтем периодичность спектра X 2 (k ) = X 2 (k ± N ) . Тогда индексы в слагаемых суммы будут «зациклены». Поэтому это выражение называют круговой или циклической сверткой. На практике, в отличие от линейной свертки, циклическую свертку реализуют не дополнением последовательности нулями, N −1 а модулярной арифметикой индексов: y (k ) = ∑ x1 (n )x2 ((k − n ) mod N ) , n =0 где операция (k − n ) mod N означает вычисление остатка от деления числа (k − n ) на число N. Поэтому в матричной форме, в отличие от матрицы линейной свертки, матрица ДПФ будет симметричной, а не теплицевой, и представлять собой циклическую свертку степеней комплексной экспоненты на сигнал. A ДПФ 1 1 1 2  2  −j −j 2 N e e N 1 2 2  −j 2 −j 4 N 1 e e N = ... ... ...  2 2 − j ( N −2 ) − j 2( N −2 ) 1 e N e N  2 2 − j ( N −1) − j 2 ( N −1)  1 e N e N ... ... ... 1 e e −j ... ... ... e 2 − j ( N −1) N 2 2( N −2 ) N ... e 2 − j ( N − 2 )2 N −j   e  2 − j 2 ( N −1)   e N   ...  2 − j ( N − 2 )( N −1) N  e  2 − j ( N −1)2 N  e 1 2 − j ( N −2 ) N 2 ( N − 2 )( N −1) N Y = A ДПФ X – ДПФ как преобразование вектора X в вектор Y. При правильном выборе длины преобразования вектор Y будет отождествляться со спектральной плотностью сигнала X. Порядок вычислительной сложности определяется N2 операций умножения и N(N-1) операций сложения. 110 3.6 ДПФ как фильтрация сигнала Рассматривая формулу прямого вычисления ДПФ (спектральной N −1 плотности сигнала) X (n ) = ∑ x(k ) e −j 2 n k N можно отметить, что выражение k =0 идентично выражению свертки, с учетом e N −1 X (n ) = ∑ x(k ) e j 2 n ( N −k ) N −j 2 n k N =e −j 2 n N N e −j 2 n (k − N ) N j =e 2 n ( N −k ) N : , т.е. коэффициенты фильтра есть комплексные значения k =0 экспоненты hn (k ) = e j 2 n k N (для каждого частотного отсчета n будет свой набор коэффициентов фильтра). Поэтому коэффициент передачи: H n ( z ) = N −1 ∑e j 2 n k N z −k k =0 представляет собой ряд – геометрическую прогрессию: H n (z ) = 1 − z−N 1− e j 2 n N z и H n ( j ) = −1 1 − e − j t s N 1− e j 2 n N e (ts – период дискретизации). − j t s  f  sin   N  sin (0.5t s N )  fs  АЧХ фильтра имеет вид: H n ( ) = . = 2n      f  n sin  0.5 t s −   sin   −   N  N    fs H2(ω) H6(ω) Hn(ω) ~0.632 H3(ω) f/fs Рис. 3.16 – АЧХ каналов ДПФ 111 Таким образом, ДПФ представляет собой N фильтров полосового типа (длиной N отсчетов), настроенных на свои центральные частоты f n = образующих  каналы ДПФ. Если заметить,  f − f n  f n   f − fn  f n  +  = N   + N N = N  fs N  fs   f s  f s  что n fs и N f n f − fn − = fs N fs и (πn – поворот фазы, на модуль не влияет), то можно записать АЧХ канала ДПФ в виде:  f − fn   sin  N f s   , H n ( ) =  f − fn   sin   f s   в математике выражение такого вида известно как ядро Дирихле: 1 n jkx 1 n sin (nx + 0.5 x ) Dn ( x ) = ∑ e = + ∑ cos(kx ) = (свертка с ядром Дирихле дает 2 k =− n 2 k =1 2 sin (0.5 x ) частичную сумму тригонометрического ряда Фурье). Совокупность АЧХ фильтров образует флуктуацию амплитудного спектра. Поэтому амплитуда частотной составляющей на частоте между центральными частотами соседних каналов ДПФ будет ослаблена, в данном примере на частоте f=0.25fs почти на 4 дБ (0.632 раз). Этот эффект известен как гребешковые искажения ДПФ (или эффект частокола). Этому эффекту подвержены все окна анализа (у каждого окна уровень гребешковых искажений свой, минимален он у плоского окна). 112 3.7 Эффекты квантования В системах цифровой обработки сигналов при совершении операций над сигналами можно выделить несколько источников возникновения погрешности от эффектов представления чисел в разрядной сетке конечной длины: − шум квантования, возникающий при оцифровке аналогового сигнала (посредством аналого-цифрового преобразования); − искажение характеристик системы при квантовании коэффициентов цифровых фильтров, нормирующих множителей и т.п.; − переполнение разрядной сетки при вычислениях (внутри системы); − округления промежуточных и/или выходных значений. 3.7.1 Шум квантования Обработка сигналов в цифровой форме возможна в виде: − чисел с фиксированной точкой (например: целочисленное представление) − чисел с плавающей точкой (представление чисел в виде m ⋅ 2 p , m – мантисса числа, p – порядок числа) При представлении чисел в целочисленном формате (r - разрядов), что характерно для работы аналого-цифровых преобразователей (АЦП), определяют: − значение младшего значащего разряда (МЗР) как отношение диапазона входного напряжения к 2r (величина означает, что без погрешности будут представляться уровни напряжения кратные величине МЗР, остальные значения будут округляться – шум квантования); − динамический диапазон, как отношение наибольшего числа (2r-1) к наименьшему (1): D = 20 lg(2 r − 1) ≈ 6.02r (дБ) при условии 2 r >> 1. Таким образом, шум квантования (ошибка округления) для идеального АЦП не превышает величины 0.5МЗР, а по своему p МЗР-1 характеру является случайным процессом с равномерной функцией плотности вероятности. 0.5МЗР Значение 0.5МЗР ошибки 113 2 1  V pp   =  r  , Vpp – размах 12  2  2 Тогда дисперсия квантователя составит: напряжения (peak to peak), для синусоидального сигнала равен двум амплитудам. Для оценки качества процессов применяют понятие SNR (SNR – signal noise ratio) – отношение сигнал шум (ОСШ) – как отношение дисперсии 2 сигнала  сигнал к дисперсии шума квантования  2 . Тогда для гармонического сигнала с амплитудой Am ОСШ имеет вид: ( ) 2  A  2 m  ≈ 6.02r + 1.761 (дБ). SNR = 10 lg  (2 A 2 r )2 12  m   На практике из-за неидеальности АЦП полагают, что ОСШ на 4-6 дБ хуже (усредненные показатели). Например, для 8-разрядного АЦП ОСШ составляет 49.9 дБ (на практике: не более 45 дБ). В частности, для типового 12-разрядного АЦП max1421 (Maxim) SNR по документации составляет 66 дБ, а по формуле получаем 74 дБ. 3.7.2 Квантование коэффициентов Каждая цифровая система характеризуется импульсной характеристикой. В частности, цифровые фильтры задаются набором коэффициентов, которые квантуются для представления внутри вычислительного устройства. При этом округление коэффициентов вызывает изменение характеристик фильтра относительно характеристики с неквантованными коэффициентами. Особенно ярко это проявляется в рекурсивных, адаптивных системах. Масштабы искажений характеристики зависят не только от разрядности коэффициентов, но и от положения полюсов функции передачи системы (искажения носят нелинейный характер: малому изменению коэффициентов может соответствовать большое изменение функции передачи). Пример. 114 Пусть фильтр описывается уравнением y (k ) = 0.25 x(k ) + 0.7612 y (k − 1) . Тогда при квантовании коэффициентов фильтра в 3-разрядной сетке получим: yкв (k ) = 0.25 x(k ) + 0.75 y (k − 1) , т.е. ошибка округления коэффициента составляет менее 1.48%. Построим функцию ошибки задания АЧХ фильтра: err = H ( ) − H кв ( ) 100% . H ( ) err, % f/fs Рис. 3.17 – Функция ошибки задания АЧХ при квантовании коэффициентов фильтра Как видно, ошибка задания АЧХ низкочастотного БИХ фильтра достигает максимального значения (3.5-4.5)% как раз в полосе пропускания фильтра. Кроме того, чем выше порядок фильтра, тем «чувствительнее» его АЧХ к квантованию коэффициентов, уже не говоря о случаях, когда квантование коэффициентов приводит к смещению полюсов функции за пределы единичной окружности (неустойчивость фильтра). 3.7.3 Переполнение разрядной сетки Другим источником искажений сигнала при прохождении через систему обработки сигналов является переполнение разрядной сетки. Дело в том, что, проектируя систему обработки сигналов, следует помнить, что на вход устройства кроме полезного сигнала поступают помехи. В результате 115 амплитуда сигнала может оказаться выше предполагаемой величины, что при операциях накопительного сложения (например: свертка), взвешивания сигнала, масштабирования может привести к результату, который не может быть представлен в разрядной сетке устройства. Например: при 8-разрядном представлении чисел без знака, сложение чисел 200 и 70 даст 14, что может привести к непредсказуемому результату работы устройства. 3.7.4 Округления В процессе вычислений (реализации алгоритма ЦОС) неизбежно возникают округления результата для представления в разрядной сетке вычислительного устройства. Особенно это проявляется в целочисленных устройствах, где надо иметь в виду, что при операциях сложения необходимо увеличивать разрядность слова (представляющего результат) на 1 разряд, а при умножениях – фактически удваивать разрядность представления результата. В случаях вычислений с плавающей точкой проблема округлений не так критична, т.к. представления в таком формате имеет большой динамический диапазон (здесь опасность потери точности заключается в операциях с числами имеющими разный порядок). Кроме того, при реализации рекурсивных систем (реализуются каскадно – сначала рекурсивная часть потом нерекурсивная) промежуточные значения могут превышать по модулю значение входного сигнала в несколько раз, что может быть источником переполнений и грубых округлений – потери точности. Чтобы избежать потери точности рекомендуется складывать сначала малые числа, затем большие. Однако на практике, как правило, неизвестно как распределены числа в массиве данных по величине, а сортировка требует вычислительных мощностей и времени. С другой стороны ошибки округления в цифровых фильтрах могут приводить к катастрофическим последствиям. Рассмотрим БИХ фильтр y (k ) = x(k ) + 0.93 y (k − 1) , который имеет единственный полюс внутри единичной окружности и значит – фильтр устойчив (т.е. если на входе конечный сигнал, то на выходе также конечный сигнал). Подадим на вход 116 фильтра скачок с амплитудой 1 В, а точность вычислений обеспечим в формате с фиксированной точкой (ddd.d). Результат (импульсная характеристика фильтра) представлен нижеследующей таблице (y(k) – выходной сигнал при точном расчете, yd(k) – выходной сигнал с учетом округлений в разрядной сетке). Видно, что при отсутствии округлений выходной сигнал экспоненциально затухает, а при наличии округлений – выходной сигнал остается постоянным – на уровне 0.7 и измениться не может, т.к. 0.93·0.7=0.651, что при округлении дает 0.7. Таблица 3. k 1 2 3 4 5 6 7 8 9 сигнал 1 y(k) 1.0 0.93 0.8649 0.804357 0.74805201 0.6956883693 0.646990183449 0.60170087060757 0.55958180966504 0.520411082988487 yd(k) 1.0 0.9 0.8 0.7 0.7 0.7 0.7 0.7 0.7 0.7 Наблюдаемое явление известно как предельный цикл. Для анализа возможности возникновения предельного цикла рассматривают выражение для эффективных значений коэффициентов знаменателя – как отношение округленного результата умножения к внутреннему состоянию фильтра (соответствующего множителя). Для нашего примера эффективные значения коэффициентов имеют вид: a = round (0.93 y (k − 1)) . y (k − 1) Видно (рис. 3.19), что модуль полюса фильтра становится большим 1, что обуславливает цикл. 117 a y(k) Рис. 3.18 – Эффективное значение коэффициента фильтра с учетом округлений 118 Раздел 4. Практические вопросы ЦОС Теоретические знания важны для анализа уже имеющихся систем и необходимы для построения алгоритмов обработки сигналов. Практическая реализация алгоритма еще на этапе выбора метода решения задачи требует рассмотрения вопроса эффективности решения, его точности, устойчивости к помехам, шумам и т.п. Поэтому на практике часто оказывается, что решение задачи «в лоб» является неэффективным, либо может приводить к существенной погрешности, т.е. недостоверному решению. Для выявления подобных проблем и определения путей решения рассмотрим гипотетическую систему обработки сигналов. Любая система обработки данных состоит из модулей (блоков) ввода сигналов и их вывода, а также собственно блока обработки – рис. 4.1. Входные сигналы Блок входных сигналов Блок обработки сигналов Блок выходных сигналов Выходные сигналы Рис. 4.1 – Обобщенная структура системы обработки сигналов Техническая реализация блоков зависит от видов сигналов и целей обработки. Блок входных сигналов обычно представляет собой (рис. 4.2): - датчик физической величины (или их множество), осуществляющий преобразование физической величины в электрический сигнал; - усилитель сигнала с датчика, т.к. часто аналоговый сигнал не соответствует по параметрам требованиям последующих цепей обработки сигнала; - фильтр низких частот, необходимый для ограничения полосы сигнала; - аналогово-цифровой преобразователь. Или же вместо цепей аналоговой обработки сигнала имеется какой-либо приемник цифровых данных (если датчик представляет собой готовое цифровое устройство или данные поступают от другой цифровой системы обработки). 119 Входные аналоговые сигналы Усилитель + ФНЧ Входные цифровые сигналы Преобразователь цифрового интерфейса К блоку цифровой обработки сигналов АЦП К блоку цифровой обработки сигналов Рис. 4.2 – Варианты источника входных данных для системы ЦОС Блок обработки сигналов представляет собой какой-либо вычислительный модуль. Состав модуля, как правило, одинаков: процессор (контроллер), память, система тактирования, цифровые интерфейсы ввода/выводы, блок питания и т.п. Блок характеризуется производительностью (как правило, числом операций в единицу времени) и потреблением энергии. Его структура может быть очень сложной и включать, кроме процессора, множество устройств обработки цифрового сигнала: аппаратные фильтры, аппаратная реализация ДПФ, различного назначения аппаратные кодеры/декодеры, вспомогательные контроллеры, процессоры и т.п. Блок выходных сигналов может представлять собой (рис. 4.3): - цифро-аналоговый преобразователь (ЦАП) с низкочастотным фильтром (ФНЧ), необходимым для сглаживания ступенчатого вида выходного сигнала ЦАП; - усилитель мощности, предназначенный для передачи сигнала по каналам аналоговой связи или для управления аналоговым устройством. Или вместо преобразования сигнала в аналоговую форму имеется какой-либо трансивер для передачи данных в цифровом виде или просто другое цифровое устройство (имеющее такой же цифровой интерфейс). 120 От блока цифровой обработки сигналов От блока цифровой обработки сигналов ФНЧ+ усилитель мощности ЦАП Выходные аналоговые сигналы Выходные цифровые сигналы Преобразователь цифрового интерфейса Рис. 4.3 – Варианты реализации блока выходных данных в системе ЦОС Пример. Рассмотрим задачу измерения сетевого синусоидального напряжения промышленной частоты (220 В, 50 Гц) рис. 4.4. Учитывая высокий уровень разности потенциалов сети, блок входных сигналов должен содержать преобразователь уровня напряжения – это может быть делитель напряжения (резистивный или емкостной) или просто понижающий трансформатор. Для ограничения полосы сигнала (для последующей оцифровки) необходимо применение аналогового ФНЧ с полосой пропускания включающей величину частоты сигнала – 50 Гц. Частота отсечки фильтра определяет минимально возможную частоту дискретизации, а уровень подавления в полосе задерживания – допустимую разрядность АЦП. Если принять, что кроме гармоники 50 Гц в сигнале нет других гармонических составляющих (а уровень шумов пренебрежимо мал), то ФНЧ можно опустить и частота дискретизации должна превышать величину 100 Гц. Задачей блока цифровой обработки сигнала (БЦОС) является определение амплитуды сигнала. Выходным сигналом является бинарный код амплитуды сигнала в 8-разрядной сетке. Сетевое напряжение uвх(t) ФНЧ uд(t) АЦП БЦОС Бинарный код амплитуды Блок входных сигналов Рис. 4.4 – Обобщенная структура измерителя сетевого напряжения 121 Рассмотрим сигнал на выходе АЦП (оцифрованный сигнал uд(t)) – рис. 4.5. Отметим, что в общем случае полученные отсчеты могут не содержать максимального значения напряжения, которое можно отождествить с Амлитуда амплитудой сигнала. T1 T2 ωt Рис. 4.5 – Период синусоиды, дискретизированный с разной частотой дискретизации Поэтому измерения амплитуды необходимо проводить косвенно – например, можно вычислить амплитуду (A) исходя из известной связи амплитуды с 1 T 2 2 A sin (t )dt = 0.5 A 2 . В случае дискретного ∫ T действующим значением: сигнала интеграл можно заменить суммой: 1 N N −1  2  nt s  = 0.5 A 2 ±  .   ∑ A2 sin 2  T n =0 При этом точность «замены» будет определяться количеством отсчетов (N) на периоде сигнала (T) и методом численного интегрирования. При применении формулы прямоугольников (т.е. можно складывать квадраты отсчетов) для достижения точности 5% необходимо, чтобы на периоде сигнала (T=20 мс) было N≥126 отсчетов, т.к. абсолютная погрешность численного интегрирования методом прямоугольников функции f ( x ) на интервале [a, b] определяется формулой: 2 ( b − a) = max f ′( x ) . 2N [a ,b ] T2 Фактически  = 2N  2 4   A T 2  , т.е. относительная погрешность  = 2N −1 . 122 Это означает, что для решения данной задачи с точностью не хуже 5% (δ<0.05) необходимо дискретизировать сигнал с частотой не менее 50 ⋅ 126 = 6300 Гц (надо помнить о других составляющих погрешности – например: погрешность квантования АЦП, поэтому реальная частота дискретизации должна быть еще больше). Тогда алгоритм, выполняемый блоком обработки (БЦОС), состоит в суммировании квадратов 126 отсчетов (или более), поступающих от АЦП, и вычислении амплитуды как корня квадратного из деленного на 252 результата суммирования: A = 0.5 125 [u д (n )]2 (с позиций эффективности вычислений ∑ 126 n =0 удобно выбрать N=128). В результате такой цифровой обработки определяется параметр сигнала – амплитуда, которая характеризует физический процесс – уровень напряжения в сети. При этом в виду независимости обрабатываемого сигнала от системы обработки необходимо учитывать следующие факторы: − начальная фаза сигнала неизвестна (относительно фронта тактирования) поэтому алгоритм обработки не должен зависеть от начальной фазы сигнала, в т.ч. от момента его поступления; − сигнал всегда содержит помехи и/или шумы, поэтому всегда надо учитывать возможность представления такого сигнала в разрядной сетке (в том числе при накопительном взвешивании и сложении отсчетов), а при наличии операций сравнения не использовать команд точного равенства (заменять командами попадания значения в интервал); − выбор частоты дискретизации сигнала должен быть основан не только на основе теоремы отсчетов (т.е. полосе сигнала), но в соответствии с требованиями к результату обработки (т.е. с учетом параметров алгоритма обработки). 123 4.1 Вычислительная сложность алгоритма Преимущества ЦОС достигается не только за счет возможности применения развитого математического аппарата (в отличие от возможностей физики при аналоговой обработке). Главным преимуществом является возможность реализации алгоритмов ЦОС на электронно-вычислительных устройствах (для простоты процессорах), которые в отличие от аналоговых устройств не подвержены влиянию внешних факторов на работу алгоритма. Это могут быть обыкновенные ПЭВМ, серверы и т.п., а также специализированные микросхемы, представляющие собой, как правило, «систему на кристалле» (микроконтроллеры, цифровые процессоры обработки сигналов и прочие). Поскольку каждый процессор, как цифровое устройство, функционирует на основе тактирующей системы, т.е. каждая команда выполняется за конечное число тактов (или циклов процессора) – за конечное время. Поэтому команды сложения, умножения выполняются также за конечное число тактов. С другой стороны разнообразие процессоров и их реализации делает неудобным использование числа тактов для сравнения производительности и сложности алгоритмов. Поэтому сложность алгоритмов сравнивают по числу операций сложения и умножения. Операции сложения, как правило, выполняются быстрее. Поэтому, исторически так сложилось, алгоритмы с преобладанием операций сложения считаются лучше, чем их же реализации на основе операций умножения (при одинаковом числе всех операций). Пример 1. Рассмотрим реализацию КИХ фильтра с линейной фазой, содержащего M коэффициентов (пусть M-четно): y (k ) = M −1 ∑ bm x(k − m ) m =0 Для вычисления одного выходного отсчета необходимо ровно M умножений и M-1 сложение на каждый отсчет выходного сигнала. Это и определяет вычислительную сложность этого простейшего алгоритма. 124 Однако, учитывая симметричность коэффициентов фильтра, количество умножений можно уменьшить в два раза, если сгруппировать слагаемые: y (k ) = 0.5 M −1 ∑ bm [x(k − m ) + x(k − m + M − 1)], при таком же числе сложений. m =0 Пример 2. Рассмотрим реализацию ДПФ (N точек): 1 X (n ) = N N −1 ∑ x(k ) e k =0 −j 2 n k N = 1 N N −1  k =0   2n   2n  k  + j sin  k   N   N  ∑ x(k )cos Очевидно, аргумент тригонометрических функций не зависит от отсчетов сигнала, как и результат вычисления функций cos и sin – поэтому их можно вычислить заранее и хранить в виде таблицы коэффициентов, в эти же коэффициенты можно включить деление на N (это возможно, если длина ДПФ известна и фиксирована). Тогда расчет действительной и мнимой частей спектра при действительном входном сигнале осуществляется по формулам: X re (n ) = 1 N  2n  ∑ x(k )cos N k  , k =0   N −1 X im (n ) = 1 N N −1  2n  k  N  ∑ x(k )sin k =0 Это требует 2N умножений и 2N-2 сложений на каждый частотный отсчет (бин) ДПФ. Всего отсчетов N, поэтому для полного вычисления N точечного ДПФ необходимо выполнить 2N2 умножений и 2N(N-1) сложений. Однако, учитывая симметрию спектра действительного сигнала, на практике вычисляют только 0.5N отсчетов. В результате вычислительная сложность алгоритма расчета спектра сигнала из N отсчетов определяется как 0.5 N ⋅ 2 N = N 2 умножений и 0.5 N ⋅ 2( N − 1) = N 2 − N сложений. Поэтому, например, в случаях, когда необходимо рассчитать не весь спектр сигнала, а только его малую часть, вместо реализации алгоритма расчета ДПФ можно применить более простые в вычислительном контексте алгоритмы. 125 4.1.1 Алгоритм Герцеля Рассмотрим импульсную характеристику отдельного канала ДПФ hn (k ) = e j 2 n k N , обозначим W = e −j 2 N . Тогда спектральную плотность сигнала можно представить в виде: N −1 X (n ) = ∑ x(k ) e j 2 n ( N −k ) N k =0 = x(0 )(W −n ) + x(1)(W −n ) N −1 N + ... + + x( N − 2 )(W −n ) + x( N − 1)(W −n ) 2 1 Для фиксированного частотного отсчета n можно вынести множитель W − n : [ ] X (n ) = x(0 )(W −n ) + x(1)(W −n ) + ... + x( N − 2 )(W −n ) + x( N − 1) (W −n ) ,  N −1 N −2 1 y N −1 затем в выражении yN-1 повторить операцию и т.д.:           X (n ) =  ...(x(0)W −n + x(1))W −n + ... W −n + x( N − 2)W −n + x( N − 1)W −n        y1    y  N − 2   y N −1 Таким образом, процесс вычисления представим в виде рекуррентного соотношения: yi = yi −1W − n + x(i ) , при y −1 = 0 . После вычисления N-1 итерации значение рекурсивной функции совпадет с частотным отсчетом: X (n ) = y N −1W − n . Рекурсивное выражение можно интерпретировать как БИХ фильтр первого порядка (но с комплексным коэффициентом) – рис. 4.6. Функция передачи этого фильтра имеет вид: Pn ( z ) = сравнить ее с функцией передачи канала 1 . Интересно 1 − W −n z −1 ДПФ: H n (z ) = 1 − z−N 1− e j 2 n N z −1 (отличаются числителем – при рассмотрении выборки N входных отсчетов задержка на N тактов не влияет на результат и формулы совпадают). 126 x(i) y(i) ( ) Y ( z ) 1 − W − n z −1 = X ( z ) + z-1 W-n Рис. 4.6 – Реализация рекурсивной функции в виде БИХ фильтра Чтобы избежать комплексных операций на каждой итерации знаменатель функции передачи умножают на комплексно-сопряженное число, тогда можно получить: −j 2 n N 1− e z −1 Pn =  2n  −1 −2 1 − 2 cos z + z  N  Этот фильтр представим каскадным соединением КИХ фильтра (рис.4.7) Pn′ = 1 − e −j 2 n N z −1 и БИХ фильтра Pn′′ = x(i) 1 −1 1 − z + z −2  2n  , где  = 2 cos .  N  vi + z-1 + y(i) z-1 -W-n α z-1 -1 БИХ часть КИХ часть Рис. 4.7 – Реализация рекурсивной функции в виде БИХ фильтра 2-го порядка Соответственно получается БИХ фильтр второго порядка, но уже с действительными коэффициентами {α, -1}. Комплексное значение числителя W =e n −j 2 n N (реализация КИХ части) участвует в процессе вычислений один раз – после последней итерации (при i=N-1), т.к. промежуточные значения yi интереса не представляют. Описанный процесс фильтрации известен как алгоритм Герцеля (был предложен Джеральдом Гёрцелем (Goertzel) в 1958 127 году). Вычислительная сложность алгоритма состоит в выполнении 2N сложений (одно из них комплексное) и N умножений (одно из них комплексное). Несомненным достоинством алгоритма является минимальные требования к памяти (для хранения промежуточных результатов) – для расчета текущего значения vi одновременно надо «помнить» только два отсчета vi-1 и vi-2: vi = xi + vi −1 − vi −2 . Пример применения алгоритма Герцеля. Пусть необходимо в дискретном сигнале произвольной формы обнаружить тон (гармоническую составляющую) с известной частотой f0=1 кГц (подобная задача возникает при обнаружении сигнала в системах декодирования, распознавания, управления, передачи данных и т.п.). Полагаем, что сигнал дискретизирован с частотой fs=10 кГц и имеется N=200 отсчетов сигнала (в памяти воображаемого устройства обработки). Тогда интересующий нас спектральный компонент имеет номер n = j  40  поэтому  = 2 cos  ≈ 1.618034 , e  200  2 n N f0 N = 20 , fs = e j 36° . Тогда необходимо вычислить vi = xi + 1.618034vi −1 − vi −2 для каждого i = 0,...,199 и X (20) = v199 − v198 e j 36° – значение спектральной плотности, амплитуда A20 = 2 X (20) . N В качестве тестовых сигналов рассмотрим сигналы: x1 (n ) = sin (100nt s ) + 5 sin (2000nt s ) содержащий тон 1 кГц и x2 (n ) = sin (100nt s ) + 5 sin (1000nt s ) не содержащий такого тона. N=200; Fs=10000; %частота дискретизации %создаем сигналы t=2*3.1415/Fs*(0:N-1); x1=sin(50*t)+5*sin(1000*t); x2=sin(50*t)+5*sin(500*t); v(1)=0; %инициализация промежуточных переменных v(2)=0; 128 for i=1:N v(3)=x1(i)+1.618034*v(2)-v(1); %рекурсивная функция for j=1:2 %осуществление задержки v(j)=v(j+1); end; end; y=v(2)-v(1)*(0.80902+0.58779i); %комплексное умножение на ej36º A20=abs(y)/N*2; %вычисление амплитуды тона Таблица 4.1 x2 Сигнал x1 X(20)=y -293.83 +404.55i 0.0006 - 0.0015i 5.0 1.6e-5 Тон обнаружен Тон отсутствует A20 Выводы Эту же задачу можно было бы решить, применяя свойства корреляционной функции. Однако в этом случае, число операций сложения и умножения было бы порядка 2N 2. N=200; Fs=10000; %создаем сигналы t=2*3.1415/Fs*(0:N-1); x1=sin(50*t)+5*sin(1000*t); x2=sin(50*t)+5*sin(500*t); x=sin(1000*t); %сигнал – эталон %расчет корреляционной функции for i=-N+1:N-1 q1=0; q2=0; for j=-N+1:N-1 m=i+j; if (m>=1)&&(m<=N)&&(j>=1)&&(j<=N) q1=q1+x(j)*x1(m); q2=q2+x(j)*x2(m); end; end; y1(i+N)=q1; %отсчеты корреляционной функции эталона и x1 y2(i+N)=q2; %отсчеты корреляционной функции эталона и x2 end; 129 y1(n) отсчеты n Рис. 4.8 – Корреляционная функция эталонного сигнала с сигналом x1 Анализируя отсчет y(0) по уровню и наблюдая симметричную функцию можно достоверно сделать вывод о наличии тона с амплитудой 2 y (0 ) =5. N y2(n) отсчеты n Рис. 4.9 – Корреляционная функция эталонного сигнала с сигналом x2 Анализируя отсчет y(0) по уровню и наблюдая несимметричную функцию можно достоверно сделать вывод об отсутствии тона. 130 4.2 Быстрое преобразование Фурье Многообразие задач требующих вычисления ДПФ требуют эффективной реализации алгоритма расчета спектральных коэффициентов, т.к. с ростом длины преобразования N количество умножений и сложений растет квадратично (при прямой реализации ДПФ). В 1965 году была опубликована статья Дж. В. Кули (Cooley) и Дж. Тьюки (Tukey), в которой описывался эффективный алгоритм вычисления ДПФ (уменьшает число умножений с N 2 до N log 2 N , но требует, чтобы N было степенью двойки). Кроме этого алгоритма известны работы Винограда, гнездовой алгоритм, алгоритм Блюстейна и т.п. Каждый из методов использует свойства ДПФ и накладывает ограничения на выбор длины преобразования. Например, эффективные алгоритмы БПФ требуют длин выборки равной простому числу или его степени. Известный гнездовой алгоритм предполагает, что длина преобразования представляется произведением простых чисел. Широкое распространение алгоритма по основанию два (Кули-Тьюки) объясняется простотой реализации в вычислительных системах. Рассмотрим ДПФ длины N=2k (k – целое число). N −1 X (k ) = ∑ x(n )e −j 2 k n N n =0 Последовательность отсчетов можно разбить на две части – с четными индексами и нечетными: X (k ) = = 0.5 N −1 ∑ x(2n ) e −j 2 k 2n N n =0 0.5 N −1 ∑ x(2n ) e n =0 −j 2 k n 0.5 N + 0.5 N −1 ∑ x(2n + 1) e −j 2 k ( 2 n +1) N n =0 +e −j 2 k 0.5 N −1 N ∑ x(2n + 1) e −j = 2 k n 0.5 N n =0 Обе суммы (по четным и нечетным индексам исходной последовательности) представляют собой ДПФ длины 0.5N. Однако, число отсчетов в частотной области также равно 0.5N: X (k ) = X чет (k ) + e −j 2 k N X нечет (k ) , k=0, …, 0.5N-1. Для получения отсчетов из диапазона 0.5N ÷ N необходимо воспользоваться периодичностью дискретного спектра: 131 X чет (k + 0.5 N ) = X чет (k ) , X нечет (k + 0.5 N ) = X нечет (k ) X (k ) = X чет (k − 0.5 N ) + e −j 2 k N −j 2 ( k − 0 .5 N ) N = X чет (k − 0.5 N ) − e X нечет (k − 0.5 N ) = X нечет (k − 0.5 N ) Далее, рассматривая полученные две последовательности (отсчеты с четными и нечетными индексами) как независимые, разбиваем их по аналогии, получая ДПФ вдвое короткой последовательности, и так далее до тех пор, пока не останется набор последовательностей ДПФ длины 2 – рис. 4.10. x(0) x1(0) x3(0)=x(0) x(1) x1(1) x3(1)=x(4) x(2) x1(2) x4(0)=x(2) x(3) x1(3) x4(1)=x(6) x(4) x2(0) x5(0)=x(1) x(5) x2(1) x5(1)=x(5) x(6) x2(2) x6(0)=x(3) x(7) x2(3) x6(1)=x(7) N=8 0.5N=4 0.25N=2 № двоичное представление 1 2 3 4 5 6 7 двоичноинвертированное представление № 000 100 010 110 001 101 011 111 4 2 6 1 5 3 7 000 001 010 011 100 101 110 111 Рис. 4.10 – Процесс разбиения последовательности в алгоритме Кули-Тьюки Итеративный процесс разбиения последовательности отсчетов на четные и нечетные индексы приводит к перестановке данных. При этом связь индексов исходной и конечной последовательностей осуществляется так называемой двоичной инверсией (пример – таблица на рис. 4.10). В процессе разбиения каждая сумма отсчетов с нечетными поворачивающий множитель типа W = e k N −j индексами 2 k N умножается на (каждая пунктирная стрелка на рисунке): X (k ) = X (k1 ) = N −1 ∑ x(n ) e −j 2 k n N n =0 0.5 N −1 ∑ x(2n ) e n =0 , k = 0 ÷ N −1 −j 2 k1 n 0 .5 N +e −j 2 k1 0.5 N −1 N ∑ x(2n + 1) e n =0 −j 2 k1 n 0 .5 N = A1 (k1 ) + WNk1 B1 (k1 ) N N  X  k1 +  = A1 (k1 ) − WNk1 B1 (k1 ), k1 = 0 ÷ − 1 2 2  132 A1 (k 2 ) = 0.25 N −1 ∑ x(4n )e −j 2 k 2 n 0.25 N n =0 +e −j 2 k 2 0.25 N −1 0 .5 N ∑ x(4n + 1)e −j 2 k 2 n 0.25 N n =0 = = А 2 (k 2 ) + W0k.52 N B 2 (k 2 ) N N  A1  k 2 +  = A 2 (k 2 ) − W0k.52 N B 2 (k 2 ), k 2 = 0 ÷ − 1 4 4  B1 (k 2 ) = 0.25 N −1 ∑ x(4n + 2)e −j 2 k 2 n 0.25 N n =0 +e −j 2 k 2 0.25 N −1 0 .5 N ∑ x(4n + 3)e n =0 −j 2 k 2 n 0.25 N = = А3 (k 2 ) + W0k.52 N B 3 (k 2 ) N N  B1  k 2 +  = A3 (k 2 ) − W0k.52 N B 3 (k 2 ), k 2 = 0 ÷ − 1 4 4  Далее каждая из последовательностей A2, B2, A3, B3 снова делятся на две. Процесс следует прекратить, когда в каждой последовательности останется всего два отсчета. Поэтому на произвольном шаге разбиения получается унифицированная операция, связывающая отсчеты текущего шага с предыдущим, которая известна как бабочка Фурье. Алгоритм с разбиением отсчетов во временной области получил название быстрого преобразования Фурье (БПФ) с прореживанием по времени. X (k i ) = X чет + WNkii X нечет Xчет Xнечет WNk i ( X k i + 0.5 N i )= X чет − WNkii X нечет Рис. 4.11 – Бабочка Фурье при прореживании по времени Пример. Даны отсчеты синусоидального сигнала (8 отсчетов на период) x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7) 0.5 2 1 0.5 2 − 0.5 2 −1 − 0.5 2 Спектр сигнала X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7) -j0.5 j0.5 133 7 X (k ) = ∑ x(n )e −j 2 k n 8 n =0 3 X (k1 ) = ∑ x(2n )e −j , k =0÷7 2 k1 n 4 +e n =0 −j 2 k1 3 8 ∑ x(2n + 1)e −j 2 k1 n 4 n =0 = A1 (k1 ) + W8k1 B1 (k1 ) X (k1 + 4 ) = A1 (k1 ) − W B1 (k1 ), k1 = 0 ÷ 3 k1 8 1 A1 (k 2 ) = ∑ x(4n )e −j 2 k 2 n 2 n =0 +e −j 2 k 2 1 4 ∑ x(4n + 1)e − jk 2 +e  − j k2 2 A1 (k 2 + 2 ) = A 2 (k 2 ) − W4k2 B 2 (k 2 ), k 2 = 0 ÷ 1 B1 (k 2 ) = ∑ x(4n + 2 )e −j 2 k 2 n 2 n =0 +e −j = x 2 + x6 e [x 1 ∑ x(4n + 3)e +e  − j k2 2 A2(0) A2(1) x(4) e − j k 2 A3(0) x(2) A3(1) x(6) e − j k 2 [x x(5) B2(1) e − j k 2 B3(0) x(3) B3(1) x(7) e − j k 2 3 −j 2 k 2 n 2 = ] + x7 e − jk2 = А3 (k 2 ) + W4k2 B 3 (k 2 ) A1(0)  − j k2 e 2  − j k2 e 2 B2(0) x(1) ] n =0 − jk 2 = + x5 e − jk2 = А 2 (k 2 ) + W4k2 B 2 (k 2 ) 2 k 2 1 4 B1 (k 2 + 2 ) = A3 (k 2 ) − W0k.52 N B 3 (k 2 ), k 2 = 0 ÷ 1 x(0) 2 k 2 n 2 n =0 = x0 + x 4 e 1 −j A1(2) A1(1) A1(3) X(0)  − j k2 e 4 X(1)  − j k2 e 4 B1(0)  − j k2 e 2  − j k2 e 2 B1(2) B1(1) B1(3) X(4) X(5) X(2)  − j k2 e 4 X(6) X(3)  − j k2 e 4 X(7) Рис. 4.12 – Структурная схема алгоритма Кули-Тьюки (N=8, прореживание по времени) Аналогичным образом можно построить алгоритм с прореживанием по частоте. 134 X (k ) = 0.5 N −1 ∑ x(n ) e 2 k n N n =0 С учетом, что e X (2k ) = −j −j 2 k 0.5 N N + 0.5 N −1 ∑ x(n + 0.5 N ) e −j 2 k ( n + 0 .5 N ) N n =0 = (− 1) можно получить: k 0.5 N −1 ∑ [x(n ) + x(n + 0.5 N )]e −j 2 k n 0 .5 N n =0 X (2k + 1) = e −j 2 k 0.5 N −1 N ∑ [x(n ) − x(n + 0.5 N )]e −j 2 k n 0 .5 N n =0 x(m) x(m+0.5N) xчет Wk xнечет Рис. 4.13 – Бабочка Фурье при прореживании по частоте Сравнение эффективности вычислений ДПФ по методу Кули-Тьюки по порядку числа операций умножения относительно прямого вычисления по формуле представлено на рис. 4.14. Число операций N2 N log 2 N N Рис. 4.14 – Эффективность БПФ по сравнению с прямым ДПФ 135 ДПФ действительного сигнала, невзирая на алгоритмы БПФ, все равно содержит избыточность, т.к. спектр сигнала симметричен (это означает, что фактически расчет одной и той же величины выполняется дважды). В связи с этим появились модификации быстрых алгоритмов под задачи, которые требуют расчета ДПФ двух разных сигналов. Известны два подхода: метод двойного преобразования и метод ДПФ сигнала удвоенной длины. 4.2.1 Двойное преобразование Если два независимых действительных входных сигнала x1 (n ) и x2 (n ) со спектрами X 1 (k ) и X 2 (k ) соответственно представить одним комплексным сигналом: x(n ) = x1 (n ) + jx2 (n ) , то выполнение одного БПФ позволяет рассчитать спектры обоих сигналов одновременно. N −1 X (k ) = ∑ ( x1 (n ) + jx2 (n ))e n =0 −j 2 k n N N −1 = ∑ x1n e −j 2 k n N n =0 N −1 + j ∑ x2 n e −j 2 k n N n =0 = N −1 N −1   2k    2k   2k    2k   = ∑ x1n  cos n  − j sin  n   + j ∑ x2 n  cos n  − j sin  n = N N N N n =0 n =            N −1 N −1    2k   2k    2k   2k   = ∑  x1n cos n  + x2 n sin  n   + j ∑  x2 n cos n  − x1n sin  n = n =0  n =0   N   N   N   N  = X re (k ) + jX im (k ) Поскольку сигнал на входе был комплексным, то его спектр, очевидно, не обладает симметрией. Однако в силу линейности свойств ДПФ он представляет собой линейную комбинацию спектров сигналов x1 и x2. N −1   2k   2k   X ( N − k ) = ∑  x1n cos n  − x2 n sin  n + n =0   N   N  N −1   2k   2k   + j ∑  x2 n cos n  + x1n sin  n   = X re ( N − k ) + jX im ( N − k ) N N n =0      Поэтому возможно по данным действительным X re (k ) , X re ( N − k ) и мнимым X im (k ) , X im ( N − k ) частям спектра комплексного сигнала x(n ) = x1 (n ) + jx2 (n ) определить спектры сигналов x1 (n ) и x2 (n ) : X 1 (k ) = 0.5[ X re (k ) + X re ( N − k )] + j 0.5[ X im ( N − k ) − X im (k )] X 2 (k ) = 0.5[ X im ( N − k ) + X im (k )] + j 0.5[ X re (k ) − X re ( N − k )] 136 Эффективность вычислений теоретически увеличивается в два раза за вычетом операций для вычисления спектров исходных сигналов из спектра комплексного сигнала. В таблице 4.2 представлено аналитическое сравнение метода двойного преобразования с расчетом двух независимых БПФ для каждого сигнала в отдельности. Таблица 4.2 – Эффективность метода двойного преобразования Вид алгоритма Операции Сложения Умножения 6 N log 2 N 4 N log 2 N 3 N log 2 N + 2 N 2 N log 2 N + N Два БПФ «Двойное» БПФ Тогда можно вычислить коэффициент эффективности в зависимости от длины преобразования (считая, что операции умножения и сложения имеют одинаковую трудоемкость вычисления):  = 5 log 2 N − 3 100% . На рис. 4.15 10 log 2 N построена полученная зависимость. η, % 48 46 44 42 40 38 36 34 0 10 N 1 10 2 10 3 10 4 10 5 10 Рис. 4.15 – Эффективность метода двойного преобразования по сравнению с расчетом двух независимых БПФ 137 4.2.2 БПФ действительного сигнала удвоенной длины Допустим, дана выборка действительного сигнала x(n ) четной длины: 2N отсчетов, спектр которого обозначим X (k ) . Сформируем комплексный сигнал, состоящий из четных и нечетных отсчетов сигнала x(n ) : v(n ) = x(2n ) + jx(2n + 1), n = 0,, N − 1 (этот сигнал в два раза короче). Определим спектр комплексного сигнала v(n ) : −j 1 N −1 V (k ) = ∑ [x(2n ) + jx(2n + 1)]e N n =0 k = 0,  , N − 1 2 k n N 1 = N N −1 ∑ x(2n )e n =0 −j 2 k n N j + N N −1 ∑ x(2n + 1)e −j 2 k n N n =0 Первое слагаемое это ДПФ масштабированного сигнала: N −1 ∑ x(2n )e −j 2 k n N n =0 = 1  X чет   2 2 Второе слагаемое это ДПФ масштабированного и сдвинутого сигнала: N −1 ∑ x(2n + 1)e −j 2 k n N n =0 = 1  X нечет  e j 0.5 2 2 Делая замену переменных: x1 (n ) = 1 1 x(2n ) и x2 (n ) = x(2n + 1) сводим к N N задаче метода двойного преобразования: N −1 V (k ) = ∑ x1 (n )e n =0 −j 2 k n N N −1 + j ∑ x2 (n )e n =0 −j 2 k n N = N −1 N −1   2k    2k   2k    2k   = ∑ x1 (n ) cos n  − j sin  n   + j ∑ x2 (n ) cos n  − j sin  n = n =0 n =0  N   N    N    N  N −1    2k   2k    2k   2k   = ∑  x1 cos n  + x2 sin  n   + j ∑  x2 cos n  − x1 sin  n = n =0  n =0   N   N   N   N  = Vre (k ) + jVim (k ) N −1 Аналогично: N −1   2k   2k   X ( N − k ) = ∑  x1 cos n  − x2 sin  n + n =0   N   N  N −1   2k   2k   + j ∑  x2 cos n  + x1 sin  n   = Vre ( N − k ) + jVim ( N − k ) n =0   N   N  138 Тогда с учетом свойств масштабирования и сдвига можно получить:  k  X re (k ) = 0.5[Vre (k ) + Vre ( N − k )] + 0.5[Vim ( N − k ) + Vim (k )]cos  − N  k  0.5[Vre (k ) − Vre ( N − k )]sin   N  k  X im (k ) = 0.5[Vim (k ) − Vim ( N − k )] − 0.5[Vim (k ) + Vim ( N − k )]sin   − N  k  0.5[Vre (k ) − Vre ( N − k )]cos  N Эффективность вычислений теоретически меньше чем у метода двойного преобразования. В таблице 4.3 представлено аналитическое сравнение метода с расчетом двух независимых БПФ для каждого сигнала. Таблица 4.3 – Эффективность метода «удвоенной длины сигнала» Вид алгоритма Операции Сложения Умножения 6 N log 2 N 4 N log 2 N 3 N log 2 N + 8 N 2 N log 2 N + 8 N Два БПФ Удвоенная длина Тогда можно вычислить коэффициент эффективности в зависимости от длины преобразования (считая, что операции умножения и сложения имеют одинаковую трудоемкость вычисления):  = log 2 N − 2 100% . На рис. 4.16 2 log 2 N + 2 построена полученная зависимость. η, % 50 40 30 20 10 1 10 2 10 3 10 4 N 10 Рис. 4.16 – Эффективность метода БПФ удвоенной длины сигнала 139 4.3 Преобразование случайного сигнала в дискретной системе Рассмотрим дискретную систему с импульсной характеристикой h(k ) . Пусть на ее вход поступает стационарный случайный процесс {x(k )} с корреляционной функцией Rx. {x(k )} y (k ) Блок обработки сигналов h(k ) Рис.4.17 – Цифровая система обработки сигналов Тогда по определению корреляционной функции: R y (n ) = ∞ ∑ y(k ) y(k + n ) . k = −∞ Выходной сигнал является результатом свертки входного сигнала и импульсной характеристики системы: y (k ) = ∞ ∑ x(k − m1 )h(m1 ) и y(k + n ) = m1= −∞ ∞ ∑ x(k + n − m2 )h(m2 ) m 2 = −∞ В результате, преобразовав произведение сумм в сумму произведений и поменяв порядок суммирования, получим: R y (n ) = = = ∞  ∞  ∞  ( ) ( ) x k − m h m ∑  ∑ 1 1   ∑ x (k + n − m2 )h(m2 )  = k = −∞   m1= −∞  m 2=−∞  ∞ ∞ ∞ ∑ ∑ ∑ x(k − m1 )x(k + n − m2 )h(m1 )h(m2 ) = m1= −∞ m 2 = −∞ k = −∞ ∞ ∞ ∑ ∑ Rx (n − m2 + m1 )h(m1 )h(m2 ) m1= −∞ m 2 = −∞ Введем новый индекс m = m2 − m1 : R y (n ) = = ∞ ∑ m = −∞ ∞ ∞ m1= −∞ m = −∞ Rx (n − m ) ∑ h(m1 )h(m2 ) = ∑ Rx (n − m ) ∞ ∑ h(m1 )h(m + m1 ) = m1= −∞ ∞ ∑ Rx (n − m )Bh (m ) m = −∞ где Bh(m) – корреляционная функция импульсной характеристики. Таким образом, корреляционная функция случайного сигнала на выходе системы представляет собой свертку корреляционной функции входного шума и корреляционной функции системы. Дисперсия выходного случайного процесса есть 140 D y = R y (0 ) = ∞ ∑ Rx (m )Bh (m ) . m = −∞ Если входной случайный процесс является белым шумом, тогда: R y (m ) = Dx Bh (m ) и D y = R y (0 ) = Dx ∞ ∑ h 2 (m ) m = −∞ т.е. при воздействии на вход системы белого шума дисперсия выходного сигнала пропорциональна сумме квадратов импульсной характеристики. 141 4.4 Усреднение сигналов Реальные данные (получаемые в результате измерений) содержат шум. Одним из простых способов снижения его влияния является усреднение. В цифровой обработке сигналов усреднение предполагает суммирование ряда отсчетов и деление суммы на их количество. Фактически усреднение по выборке сигнала, длина которой определяется количеством отсчетов в сумме – порядком усреднения, является математическим ожиданием: xср = 1 Количественной мерой флуктуации является дисперсия:  = N 2 1 N N −1 ∑ x(k ) . k =0 ∑ [x(k ) − xср ] N −1 2 , k =0 где σ – среднеквадратическое отклонение (СКО). Различают когерентное и некогерентное усреднение. Когерентное усреднение заключается в том, что усредняемые значения выборок сигнала синхронизированы по фазе, в то время как шум некоррелирован с сигналом. Поэтому когерентное усреднение приводит к уменьшению дисперсии шума и как следствие – к улучшению ОСШ. Рассмотрим выборку сигнала x(k), включающую N отсчетов. Получим M таких выборок в предположении, что все они синхронизированы относительно фазы сигнала (когерентны). Таким образом, получены данные: Набор 1: x1(0), x1(1), x1(2), …, x1(N-1) Набор 2: x2(0), x2(1), x2(2), …, x2(N-1) … Набор M: xM(0), xM (1), xM (2), …, xM (N-1) Рассматривая i-ый отсчет в каждом наборе можно утверждать, что он является суммой некоторого точного значения сигнала и шумовой составляющей. Поэтому величина xср (i ) = 1 M M ∑ xk (i ) = k =1 1 M ∑ [xточное (i ) + k ] = xточное (i ) + M k =1 1 M M ∑ k k =1 при случайном процессе (шуме) ηk с нулевым средним величина xср(i) стремится к точному значению сигнала xточное(i), где i∈[0; N-1]. 142 Для оценки эффективности усреднения сравнивают дисперсию усредненного сигнала σср2 с дисперсией исходного сигнала σх2 : 1  = М 2 ср ∑ [xср (i ) − xточное (i )] M 2 k =1 ∑ [ xk (i ) − xточное (i )] M 1 = М 2 1 M    = ∑ k  M 2 ∑  k  k =1  k =1   M M 1 М где i – номер отсчета в наборе (выборке) с номером k; k =1 = 2 ∑ [k ]  2x = 2 1 М 1 ∑M k =1  M 2 k =1 xточное(i) – по смыслу есть математическое ожидание, вокруг которого флуктуирует сигнал xk(i) и к которому стремится среднее значение xср(i). M  k  ∑ k =1  Величину 2 в соответствии с биномом Ньютона можно 2 преобразовать в выражение: M M −1 M M  2  =  + 2 ∑ k ∑ ∑ i  j . Рассматривая k ∑ k =1 k =1 i =1 j =i +1  второе слагаемое можно переобозначить переменную суммирования j: M −1 M M −1M −i +1 ∑ ∑ i  j = ∑ ∑ i i+ j , где i =1 j =i +1 i =1 j =1 отсчеты шума k с индексами k>M равны нулю (требование формальное). Поэтому можно продолжить равенство и записать: M −1 M ∑ M M ∑ i  j = ∑∑ i i+ j . i =1 j =i +1 Полученное выражение есть сумма отсчетов j =1 i =1 автокорреляционной функции белого шума (с исключенной точкой «нулевого сдвига», соответствующей энергии сигнала): M M M j =1 i =1 m =1 ∑∑ i i+ j = ∑ R (m ) . Поэтому с учетом свойств отсчетов белого шума (любая пара некоррелирована) можно 2 M M  записать: ∑ k  = ∑ 2k . k =1  k =1  2 Тогда M  k  2 ∑  ср k =1  = 1 , = M 2 x M ∑ 2k M т.е.  2x  = M 2 ср – дисперсия шума в k =1 усредненном сигнале в M раз меньше дисперсии шума в исходном сигнале. Таким образом, коэффициент улучшения ОСШ сигнала при когерентном 143 усреднении равен M (как отношение ОСШ сигнала на входе фильтра к ОСШ сигнала на его выходе). На практике когерентное усреднение активно применяется в современных цифровых осциллографах. Усреднение некогерентных данных (известно как среднеквадратическое или скалярное усреднение) применяется в частотной области (к модулю спектра), т.к. во временной области его применение малоэффективно (подобно ФНЧ с посредственными характеристиками). Например, измеряется спектр периодического сигнала, содержащего некоррелированный шум – рис. 4.18. X(f) Спектр сигнала без шума f, Гц k=1 k=2 k=3 k=4 Рис.4.18 – Спектры смеси сигнала с шумом при уровне некогерентного усреднения k 144 Рассчитали ДПФ сигнала – получили оценку его спектра X 1 (k ) . Затем, получили другую выборку того же сигнала (может пересекаться с предыдущей выборкой) и снова рассчитали ДПФ – получили оценку спектра X 2 (k ) и так далее, пока не будет получено X M (k ) . Поскольку модуль спектра не зависит от фазы сигнала, то усреднение в частной области на каждой частоте (бине с номером k) становится когерентным: X ср (k ) = во временной 1 M M ∑ X m (k ) . Но в отличие от когерентного усреднения m =1 области среднее значение шумовых составляющих (присутствующих в X m (k ) ) стремится не к нулю, а к некоторому значению (математическое ожидание положительных величин отлично от нуля). При этом среднее значение точной составляющей бина остается без изменений, т.е. данный метод уменьшает дисперсию шума, но не улучшает отношение мощности сигнала к мощности шума. Выигрыш по критерию ОСШ при усреднении в частотной области составляет M раз, где M – порядок усреднения. Реализация усредняющего фильтра возможна в виде как КИХ фильтра – рис. 4.19, так и в виде БИХ фильтра – рис. 4.20 (эти схемы эквивалентны при нулевых начальных условиях). Усредняющие фильтры также называются однородными (когда все коэффициенты фильтра одинаковы h(n) = 1 N – импульсная характеристика фильтра), что обеспечивает минимум величины ∑ h 2 (n ), где величина (N-1) – порядок усредняющего фильтра. n 145 x(n) z-1 x(n-1) z-1 … z-1 x(n-M) … 1 H (z ) = N N −1 ∑z Σ −k N-1 y(n) k =0 Рис. 4.19 – Усредняющий КИХ фильтр x(n) z-1 N-1 1 1 − z−N H (z ) = N 1 − z −1 x(n-1) -1 Σ y(n-1) z-1 y(n) Рис. 4.20 – Усредняющий БИХ фильтр 146 4.4.1 Уменьшение шума квантования АЦП Наличие шума квантования АЦП в реальном цифровом сигнале принципиально, а потому и методы его уменьшения актуальны. Как известно, дисперсия шума квантования определяется числом разрядов АЦП, т.е. МЗР 2 величиной младшего значащего разряда (МЗР):  = . При этом 12 2 полагаем, что квантование – случайный процесс с равномерным распределением и отсчеты шума некоррелированы (спектр шума квантования равномерен). Тогда спектральная плотность мощности (СПМ) оказывается равномерно распределенной функцией в диапазоне частот от -Fs до +Fs (Fs – 2 частота дискретизации), т.е. СПМ = . Тогда для уменьшения СПМ шума Fs квантования возможно: − увеличить разрядность АЦП (тривиальное решение, требующее материальных затрат или вовсе является невозможным); − увеличить частоту дискретизации Fs (этот подход называется сверхдискретизацией). Суть метода сверхдискретизации состоит в «размазывании» шума квантования в большей полосе частот. Поэтому при последующей НЧ фильтрации сигнала (ограничении спектра, например, под условия теоремы отсчетов) вклад шумовой составляющей (шума квантования) в сигнал уменьшится, т.к. часть ее энергии шума будет подавлена фильтром. Графически процесс распределения шума квантования при сверхдискретизации представлен на рис. 4.21. Спектр сигнала Спектр шума -0.5Fs -0.5Fsn 0.5Fs частота 0.5Fsn Рис. 4.21 – Шум квантования при сверхдискретизации (FsM частотных отсчетов аппроксимирующих тот же спектр. Следует отметить, что существует неправильное мнение о том, что дополнение нулями улучшает разрешение ДПФ, поскольку оно увеличивает длину преобразования. На деле дополнение нулями позволяет получить лишь интерполированное преобразование более сглаженной формы. В частности, оно устраняет неопределенности, обусловленные наличием узкополосных компонент сигнала, частоты которых лежат между соседними бинами спектра. Рассмотрим пример вычисления 32-точечного ДПФ синусоидального сигнала, в том числе при дополнении его нулями. Рассматривая спектр сигнала рис.4.27А невозможно сказать, сколько гармоник присутствует в спектре – две (48 Гц и 72 Гц) или одна (два ненулевых бина могут соответствовать одной гармонике расположенной на частотной оси между ними). Однако с ростом длины преобразования аппроксимация спектра (на рис.4.27 показана красной штриховой линией) становится более подробной – по рис. 4.27В и рис.4.27Г можно уверенно говорить именно о двух гармонических составляющих спектра. Однако значения амплитуд и частот гармонических составляющих маскируются влиянием добавленных нулей, т.е. фактически искажением исходного сигнала. 152 без дополнения нулями (N=32) А f n дополнение нулями до N`=64 Б f n дополнение нулями до N`=128 В f n дополнение нулями до N`=256 Г n f Рис.4.27 – Сигналы (слева) и их спектры (справа) при дополнении нулями 153 4.5.2 Растяжение участка спектра Метод растяжения участка спектра сигнала позволяет увеличить разрешающую способность ДПФ, но требует выделения части частотной полосы сигнала (т.е. игнорирования спектральных составляющих за пределами выделенной полосы). Соответственно метод оправдан, когда необходимо высокое разрешение ДПФ в небольшой области частотной полосы сигнала. Допустим, в сигнале x(n) с полосой частот [0, fв] необходим подробный спектральный анализ в области частот [f1, f2], причем 016 кГц. Полоса пропускания ФНЧ должна быть 16 кГц, положим ширину переходной полосы фильтра менее 16 кГц. Тогда получим: fд=32 кГц и M=16. Для реализации ДПФ с разрешением 1 кГц потребуется выборка N`=32 отсчетов, т.е. 160 комплексных умножений. Если ФНЧ реализовать в виде КИХ фильтра с линейной фазой порядка K, то фильтрация потребует Pf=0.5(K+1)N` комплексных умножений. Для реализации ФНЧ с полосой 16 кГц и переходной полосой 16 кГц и частоте дискретизации сигнала 512 кГц необходим фильтр порядка K≥27. 155 Итого алгоритм требует 512+160+0.5·28·32=1120 комплексных умножений. Поэтому выигрыш (в смысле вычислительной эффективности по количеству умножений) по сравнению с прямым вычислением ДПФ требуемого разрешения составляет более 75.6%, т.е. решение данной задачи с помощью рассмотренного метода (структура представлена на рис. 4.32) эффективнее более чем в 4.1 раза. x(t) АЦП fs=512 кГц x(n) x1(n) × БПФ ФНЧ ↓16 x2(n) x3(n) fд=32 кГц X(k) e − j 2 f c nt s Рис.4.32 – Структурная схема устройства спектрального анализа участка спектра сигнала 156 4.5.3 Алгоритм локализации спектральных пиков Расчет ДПФ сигнала с неизвестной частотой (даже с применением оконных функций или с передискретизацией сигнала) фактически всегда содержит эффект размытия. Это означает, что при попытке точного определения частоты гармоники или ее амплитуды неизбежно возникает проблема точности измерения (поскольку максимум, который ассоциируется с дельта-функцией, отображающей в спектре гармоническую компоненту сигнала, оказывается смещенным по частоте и искаженным по амплитуде) – рис. 4.33 (синусоидальный сигнал с амплитудой 1 В, частотой 100 Гц, частота дискретизации Fs=2 кГц, длина ДПФ N=128 отсчетов: разрешение по частоте = Fs = 15.625 Гц ). N амплитуда 93.75 Гц 109.375 Гц f, Гц Рис. 4.33 – Проблема оценки частоты и амплитуды синусоидальной компоненты Очевидно, частота синусоидальной компоненты лежит между спектральными бинами с частотами 93.75 Гц и 109.375 Гц (или номера частотных отсчетов 6 и 7 соответственно, нумерация с нуля). Конечно, можно применить методы по уменьшению размытия:  увеличить длину ДПФ (уменьшив разрешение по частоте) – не гарантирует решения проблемы и требует значительных вычислительных затрат;  применить метод дополнения нулями – не решает задачи, т.к. частоты спектральных пиков при этом смещаются, а амплитуды искажаются; 157  применить подходящее окно анализа – размытие уменьшается, но не достаточно – рис. 4.34. 1 амплитуда 0.8 0.6 0.4 0.2 f, Гц 50 100 150 200 250 300 Рис. 4.34 – Фрагмент спектра сигнала (рис. 4.33) при применении окна Решение (без увеличения выборки сигнала, применения окно и т.п.) может быть найдено с помощью аппарата аппроксимации функции. При этом необходимо быть уверенными, что в окрестности аппроксимируемого спектрального пика (в данном примере 100 Гц) нет других частотных составляющих сигнала. Положим, что индекс (номер отсчета), соответствующий максимуму спектрального пика (т.е. гармонической составляющей сигнала с частотой Fp = Fs m p ) имеет вид: m p = m + Re{}, N где m – индекс (номер отсчета), соответствующий максимуму в рассчитанном ДПФ (содержащим размытие спектра) – является целым числом; δ – характеризует смещение спектрального пика (гармоники сигнала) от текущего максимума. Определив положение максимума – индекс m по рассчитанному ДПФ (результат подвержен размытию спектра), можно выделить три комплексных отсчета: X(m-1), X(m), X(m+1), т.е. значения спектра в точке максимума и значения в соседних точках. По этим трем равноотстоящим точкам можно построить полином второй степени. Затем по найденным коэффициентам уравнения определить положение его единственного экстремума (максимума). 158 Положим коэффициенты a, b, c – коэффициенты полинома второй степени am 2 + bm + c = X (m ) . Тогда справедливы уравнения: am 2 + bm + c = X (m )  2 a(m + 1) + b(m + 1) + c = X (m + 1)  2 a(m − 1) + b(m − 1) + c = X (m − 1) Решая систему уравнений, определим коэффициенты a, b и вычислим положение экстремума:  = mopt = − b X (m + 1) − X (m − 1) = . 2a X (m + 1) + X (m − 1) − 2 X (m ) Тогда, для примера рис. 4.33, получаем следующие результаты: m p = m + Re{}=6+0.401825 Fp = Fs m p =100.0285 Гц N 159 4.6 Эффективная реализация КИХ фильтров Как уже говорилось, применение КИХ фильтров имеет ряд преимуществ, однако при решении практических задач порядок фильтра получается большим (десятки, сотни). Это требует решения вопросов эффективной реализации свертки. Одним из простых способов является метод, в основе которого лежит свойство ДПФ: спектр свертки сигналов – есть произведение спектров исходных сигналов. Тогда выполнив ДПФ от сигнала и импульсной характеристики фильтра (т.к. у КИХ фильтра она не изменяется, ее ДПФ вычисляют однократно и хранят в памяти устройства обработки сигнала), выполняют перемножение полученных спектров и затем вычисляют обратное ДПФ, которое будет представлять собой свертку импульсной характеристики фильтра и сигнала. Положим, дана импульсная характеристика фильтра h(k), порядок фильтра K. Требуется отфильтровать сигнал x(n), длина выборки N. x(n) x(n)*h(k) y(n) X(k) x(n) БПФ × ОБПФ y(n) H(k) Рис.4.35 – Прямая (слева) и эффективная (справа) реализация КИХ фильтрации Прямое решение задачи требует (K+1)N операций умножения. Реализация посредством алгоритма БПФ требует вычисления прямого ДПФ ( N log 2 N операций умножения), перемножения спектров (4N операций умножения), вычисления обратного ДПФ ( N log 2 N операций умножения). Итого: 2 N log 2 N + 4 N (отметим, величина не зависит от порядка фильтра). Рассчитаем коэффициент эффективности как отношение числа операций умножения при прямом решении задачи к числу операций умножения при решении задачи посредством применения алгоритма БПФ:  = (K + 1) 2 log 2 N + 4 . Зададимся вопросом, при каком порядке фильтра K имеет место улучшение, т.е. η>1. Это выполняется если K > 2 log 2 N + 3 . Например: при N=1024, K>23. 160 4.7 Вычисление свертки секционированием Было показано, что сложность арифметических операций при реализации КИХ фильтров с помощью ДПФ не зависит от порядка фильтра, в то время как сложность реализации с непосредственным вычислением свертки пропорциональна порядку фильтра. Если порядок фильтра мал, можно ожидать, что более эффективной будет реализация методом прямой свертки, однако по мере увеличения порядка фильтра реализация с ДПФ станет более эффективной. Для фильтров очень высоких порядков выигрыш может достигать десятков и сотен раз. С другой стороны, реализация с ДПФ требует значительных объемов памяти и соответственно растет время обработки блока (требуется накопить буфер отсчетов заданной длины). Методы секционирования свертки являются компромиссным решением. Суть их заключается в том, что операция свертки выполняется над секциями или блоками данных с использованием БПФ для ускорения процесса вычислений. Ограничение размера секций уменьшает объем требуемой памяти, а использование алгоритмов БПФ сохраняет вычислительную эффективность процедуры, компенсируя «усилия» на организацию процесса и дополнительную вычислительную нагрузку. 4.7.1 Метод перекрытия с суммированием Разделим сигнал x(n) (длина выборки N может быть и бесконечной) на секции длиной по M отсчетов. Мысленно выделим секцию с индексами mM≤n≤(m+1)M-1 (m – произвольный номер секции) – рис. 4.36. Отсчеты сигнала x(n) 0 … M-1 … секция секция 1 x`m(i) mM … … секция m mM+1 секция m+1 mM+2 … … N-1 секция N/M-1 mM+i … mM+M-1 Рис.4.36 – Разбиение сигнала на непересекающиеся секции длиной M отсчетов 161 Расчет линейной свертки импульсной характеристики фильтра длиной K отсчетов и секцией сигнала длиной M отсчетов дает M+K-1 отсчетов (рис.4.37): mM ≤ n ≤ (m + 1)M − 1  x(n), , xm′ (i ) =  иначе 0, K −1 y m′ (m1 ) = ∑ h(k ) xm′ (m1 − k ) , 0≤i≤M-1. n = mM + m1 , k =0 m1 ∈ [0, M − 1] – смещение в секции, m ∈ [0, N / M − 1] – номер секции. Отсчеты сигнала x(n) 0 … M-1 … 2M-1 секция 1 секция 2 … … … секция m секция m+1 … N-1 секция N/M h(k) K отсчетов Выходные отсчеты y`(n) K-1 отсчет отклик 0 отклик 1 … M+K-1 отсчет отклик m отклик m+1 сложить … отклик N/M-1 Рис.4.37 – Метод перекрытия с суммированием В силу линейности операции свертки выходной сигнал есть сумма откликов фильтра h(k) на каждую секцию: y (n ) = N / M −1 ∑ ym′ (n mod mM ) = m =0 N / M −1 K −1 ∑ ∑ h(k ) xm′ ((n m =0 k =0 mod mM ) − k ) . Для достижения вычислительной эффективности свертку в каждой секции выполняют с помощью эффективных алгоритмов, например: с помощью БПФ. А пересекающиеся результаты вычисления сверток суммируются по границам индексов исходных секций. 162 4.7.2 Метод перекрытия с накоплением При подаче возмущения на вход системы (фильтра) с конечной импульсной характеристикой h(n) длиной K отсчетов ее отклик имеет длину K-1 отсчет. Таким образом, при подаче на вход системы сигнала длиной N>>K в начале результата расчета свертки (относительно нулевого отсчета) будет иметь место переходной процесс длительностью K-1 отсчет. Для устранения влияния переходного процесса в выходном сигнале необходимы K-1 предыдущих отсчетов сигнала, обеспечивающих «запуск» расчета свертки к моменту времени, соответствующему отсчету y(n) выходного сигнала – рис.4.38. … x(n) x(n+1) x(n+2) … x(n+M-2) x(n+M-1) x(n+M) … … 0 0 0 x(n) x(n+1) x(n+2) … x(n+M-2) x(n+M-1) 0 0 … h(K-1) … h(2) h(1) h(0) K −1 y (n ) = ∑ h(k )x(n − k ) y(n) y(n+1) … y(n+K-1) … y(n+M-1) … y(n+M+K-1) k =0 Отсчеты искажены переходным процессом … x(n-K+1) … x(n-1) x(n) x(n+1) x(n+2) … x(n+M-2) x(n+M-1) x(n+M) … x(n-K+1) … x(n-1) x(n) x(n+1) x(n+2) … x(n+M-2) x(n+M-1) 0 0 … h(K-1) … h(1) h(0) y(n) y(n+1) … y(n+K-1) … y(n+M-1) … y(n+M+K-1) Рис. 4.38 – Идеология метода перекрытия с накоплением Таким образом, желая в выходном сигнале получать секции данных расчета свертки сигналов длиной N, необходимо входной сигнал x(n) разбить на неперекрывающиеся секции длиной N-K+1 отсчет и добавить перекрытие в K-1 отсчет слева – рис. 4.39. При этом для получения выходного сигнала 163 необходимо собрать результат: составляя выходные секции методом «стык в стык». Отсчеты сигнала x(n) секция 1 N-1 отсчет секция 2 h(n) … секция m K-1 отсчет секция m+1 … Выходные отсчеты y(n) отклик 1 отклик 2 … отклик m отклик m+1 … Рис.4.39 – Метод перекрытия с накоплением Кроме того, расчет линейной свертки можно производить только для N отсчетов (вместо N+K-1), в противном случае – правые K-1 отсчеты необходимо отбросить. Для достижения вычислительной эффективности свертку в каждой секции выполняют с помощью эффективных алгоритмов, например: с помощью БПФ. 164 4.7.3 Скользящее ДПФ Часто возникает задача получения спектра сигнала с той же интенсивностью, с которой поступают отсчеты сигнала. Одним из решений является применение алгоритма Герцеля. Однако в случае необходимости расчета спектра по приходу каждого нового отсчета сигнала часть расчетов (операций) повторяется, поэтому в этом случае возможно уменьшить число операций. Такой алгоритм известен как скользящее ДПФ (СДПФ). Алгоритм СДПФ вычисляет значение одного бина ДПФ длины N, используя результаты предыдущего расчета. 1 Для m-го бина известно: X (m ) = N N −1 ∑ x(n ) e −j 2 mn N . n =0 Рассмотрим процесс вычисления ДПФ от сигнала при сдвиге на один отсчет – рис. 4.40. x(0) x(1) … x(N-1) x(N) … x(p) … x(p+N-1) … X0(m) X1(m) … … … Xp(m) … … Рис. 4.40 – Соответствие выборок при скользящем ДПФ По рис. 4.40 видно, что из скользящего окна крайний левый отсчет удаляется и справа прибавляется новый отсчет. Поэтому записав выражение ДПФ для сигнала, начиная с произвольного индекса p, и для того же сигнала, начиная с индекса p+1, получим связь между спектрами двух выборок длины N: 1 X p (m ) = N N −1 ∑ x( p + n ) e −j 2 mn N n =0 165 N −1 1 X p +1 (m ) = N = e ∑ x( p + 1 + n ) e −j 2 mn N 2 − j m ( n −1) 1 N = ∑ x( p + n ) e N = N n=1 n =0 j 2 m N N N ∑ x( p + n ) e j 2 m N 2 − j mn  N −1  =  ∑ x( p + n ) e N − x( p ) + x( p + N ) N  n =0  2 − j mn N e n =1 Перегруппировав слагаемые и подставляя выражение для Xp(m), получим: 1 X p +1 (m ) =  N N −1 ∑ x( p + n ) e −j 2 mn N n =0 x( p + N ) − x( p ) j N m + = e N  2 x( p + N ) − x( p ) j N m  =  X p (m ) +  e N  2 В полученном выражении присутствует слагаемое x(p+N) – отсчеты «из будущего». Однако надо учитывать, что расчет ДПФ подразумевает накопление буфера в N отсчетов, т.е. задержку в N отсчетов. Учитывая этот факт, без потери общности можно ввести новую переменную (индексацию), учитывающую временной сдвиг: x( p ) − x( p − N ) j N m  X p (m ) =  X p −1 (m ) +  e N  2 Важным фактом применения этого выражения является: − для расчета очередного спектра сигнала достаточно двух действительных сложений, действительного умножения и комплексного умножения; − расчет очередного спектра не зависит от величины N. Полученный результат можно представить в виде соединения фильтров (гребенчатого и комплексного резонатора) – рис. 4.41. e j 2 m N x(n) + N-1 + × Xp(m) Xp-1(m) z-N z-1 гребенка Комплексный резонатор Рис. 4.41 – Структура однобинового СДПФ 166 При необходимости расчета нескольких бинов потребуются несколько комплексных резонаторов на каждую частоту m, включенных параллельно. Имея разностное уравнение можно найти передаточную функцию: 2  X ( z ) − X ( z )z − N  j N m −1 X p ( z ) =  X p ( z )z + e N   H (z ) = (1 − z )e −N 1− e j 2 m N j 2 m N z −1 Этот комплексный фильтр имеет N нулей, расположенных равномерно на единичной окружности в z-плоскости и единственный полюс, подавляющий ноль в точке z = e j 2 m N . Полученная характеристика повторяет зависимость, полученную при рассмотрении алгоритма Герцеля, но требует меньших вычислительных затрат при расчете очередного спектра. Эквивалентность рассматриваемого алгоритма расчету одного бина ДПФ приводит к необходимости рассмотрения вопроса утечки спектра. В соответствии с теоремой о свертке применение окна можно осуществить и в частотной области. Пример применения окна Хэмминга представлен на рис. 4.42. x(n) + N-1 z e -N Xp(m-1) + × j 2 ( m −1) N z-1 × -0.23 гребенка Xp(m) Xp(m) + e j 2 m N z-1 Xp(m+1) + e × × × j 2 ( m +1) N z-1 + 0.54 × -0.23 Рис. 4.42 – Структура СДПФ с окном Хемминга 167 Обращая внимание на комплексный резонатор, имеющий полюс на единичной окружности, можно говорить о проблеме устойчивости СДПФ. Действительно, при реализации алгоритма в устройствах с конечной разрядной сеткой неизбежно возникают округления, которые приводят к проблемам аналогичным при реализации БИХ фильтров. Для предотвращения возможного «выхода» полюса за пределы единичного круга вводят коэффициент затухания r (величина близка к 1): H уст ( z ) = (1 − r N z −N 1− re )r e 2 j m N j 2 m N z −1 Соответствующая выражению структура фильтра представлена на рис. 4.43. re j 2 m N x(n) + N-1 + × Xp(m) z-N Xp-1(m) × rN гребенка z-1 Комплексный резонатор Рис. 4.43 – Структура устойчивого однобинового СДПФ Другим методом решения проблемы устойчивости вычисления отдельно взятого бина является уменьшение наибольшей (действительной или мнимой) коэффициента e j 2 m N составляющей на величину младшего разряда. 168 4.8 Фильтры Фарроу Одной из основных задач анализа сигналов, как уже отмечалось, является эффективное вычисление спектров и, в частности обеспечение условий отсутствия эффекта размытия спектра. Одним из способов решения задачи является передискретизация сигнала (раздел 2.3). В общем случае передискретизация является задачей интерполяции. Этот же механизм можно применять для решения другой важной задачи – задержка дискретного сигнала на время некратное периоду дискретизации (т.е. на «доли индекса» отсчета). Задачу математической интерполяции функции удобно рассматривать в виде поиска коэффициентов полинома проходящего по отсчетам сигнала (узлам интерполяции). Положим, имеется N отсчетов сигнала x(n). Допустим необходимо выполнить передискретизацию сигнала (resampling) x(n). Для решения задачи представим сигнал (или его участок) полиномом, что позволит вычислить любые значения функции между узлами интерполяции – т.е. получить сигнал с другой частотой дискретизации или сигнал сдвинутый на «долю индекса отсчета» – рис. 4.44 (черным для наглядности показан аналоговый сигнал, красным – имеющийся дискретный сигнал с периодом дискретизации ts, синими ромбиками отмечены отсчеты дискретного сигнала сдвинутого на время v=0.4ts, зеленым пунктиром показан сигнал с другим периодом дискретизации tsn=1.2ts). ts x(n) v tsn ts n Рис. 4.44 – Процесс интерполирования дискретного сигнала Известно, что через N точек (отсчетов) проходит единственный полином степени N-1. Таким образом, выходной сигнал представим в виде: N −1 y (t ) = PN −1 (t ) = ∑ ak t k , где коэффициенты ak необходимо определить по k =0 исходному сигналу x(n). 169 Для этого необходимо решить систему уравнений (значения интерполирующего полинома в узлах интерполяции совпадает со значениями сигнала): N −1  ( ) x = ∑ a k 0 k = a0  k =0  N −1 N −1  x(1) = a 1k = a ∑ k ∑ k  k =0 k =0  N −1  k ( )  x 2 = ∑ ak 2 k =0    N −1  x( N − 1) = ∑ ak ( N − 1)k  k =0  Однако появление каждого нового отсчета сигнала требует нахождения новых коэффициентов полинома, т.е. решения новой системы уравнений, имеющей вычислительную сложность порядка N3 операций умножения за период дискретизации. Подобное решение должно считаться неэффективным в плане необходимой производительности вычислительного устройства. Например, реализация интерполяции N=4 порядка для определения коэффициентов полинома требуется 64 умножения и 3 умножения для вычисления каждой точки нового сигнала. Если же, к примеру, частота дискретизации сигнала 44.1 кГц (обработка музыкального контента), то производительность вычислительной системы должна превосходить 44100х67=2954700 умножений в секунду (на каждый канал стерео сигнала), что требует применения серьезных вычислительных процессоров и соответствующего потребления энергии, что просто нерационально. На самом деле задачу интерполяции полиномом можно упростить и, в частности, достичь решения рассмотренной задачи (N=4) при использовании только трех умножений. В основе подхода лежит представление сигнала с помощью ортогональных полиномов, например – полиномов Лагранжа: N −1 y (t ) = ∑ x(n )LNn −1 (t ) n =0 170 N −1 N −1 n Полином Лагранжа L (t ) = ∏ (t − t k ) k =0 ,k ≠ n N −1 ∏ (t n − t k ) это полином N-1 степени, равный k =0 ,k ≠ n единице в точке t=nts – момент дискретизации n-го отсчета и равен нулю в другие моменты дискретизации. На рис. 4.45 представлены кубические полиномы Лагранжа (N=4). Полином, выделенный красным цветом, равен единице в точке t=-2 и равен нулю в других точках. Полином, выделенный синим цветом, равен единице в точке t=-1 и равен нулю в других точках и т.д. L3(t) t Рис. 4.45 – Полиномы Лагранжа (N=4) Для получения эффективной вычислительной процедуры распишем полиномы Лагранжа для случая N=4 и подставим необходимые моменты времени: t0=-2; t1=-1; t2=0; t3=1: (t − t1 )(t − t 2 )(t − t3 ) = (t + 1)t (t − 1) 1 = − (t 3 − t ) (t0 − t1 )(t0 − t 2 )(t0 − t3 ) (− 2 + 1)(− 2)(− 2 − 1) 6 (t − t0 )(t − t 2 )(t − t3 ) = (t + 2)(t )(t − 1) = 0.5t 3 + 0.5t 2 − t L13 (t ) = (t1 − t0 )(t1 − t 2 )(t1 − t3 ) (− 1 + 2)(− 1)(− 1 − 1) (t − t0 )(t − t1 )(t − t3 ) = (t + 2)(t + 1)(t − 1) = −0.5t 3 − t 2 + 0.5t + 1 L32 (t ) = (t 2 − t0 )(t 2 − t1 )(t 2 − t3 ) (2)(1)(− 1) (t − t0 )(t − t1 )(t − t2 ) = (t + 2)(t + 1)(t ) = 1 t 3 + 0.5t 2 + 1 t L33 (t ) = (t3 − t0 )(t3 − t1 )(t3 − t2 ) (1 + 2)(1 + 1)(1) 6 3 L30 (t ) = 171 Произвольный отсчет интерполируемого сигнала является взвешенной суммой 3 кубических полиномов: y (t ) = ∑ x(n )L3n (t ) и при подстановке полученных n =0 выражений для полиномов Лагранжа получим: 1 3 (t − t )x(0) +  1 t 3 + 1 t 2 − t  x(1) + − 1 t 3 − t 2 + 1 t + 1 x(2) + 6 2 2 2   2  1 1  1 +  t 3 + t 2 + t  x(3) 2 3  6 y (t ) = − Раскрыв скобки и сгруппировав слагаемые относительно степеней t, получим искомый полином: y (t ) = a3t 3 + a2 t 2 + a1t + a0 , осуществляющий расчет значения сигнала в любой момент времени с погрешностью определяемой порядком интерполяции, т.е. решающей задачу передискретизации или задержки сигнала на момент времени некратный периоду дискретизации: 1 1 1 1 a3 = − x(0 ) + x(1) − x(2 ) + x(3) 6 2 2 6 1 1 a2 = x(1) − x(2 ) + x(3) 2 2 1 1 1 a1 = x(0 ) − x(1) + x(2 ) + x(3) 6 2 3 a0 = x(2 ) Каждый коэффициент полинома является взвешенной суммой отсчетов сигнала, как и процесс расчета значения полинома можно представить в виде КИХ фильтров – рис. 4.46. Полученный фильтр известен как фильтр Фарроу (аналогично рассмотренному примеру можно получить фильтр любого порядка). В представленной структурной схеме всего 9 умножений на расчет коэффициентов полинома (умножения на ноль, на плюс/минус единицу не считаем) и 3 умножения на расчет значения сигнала (по сравнению с 67 умножениями это достижение). Однако нетрудно видеть, что отсчеты исходного сигнала умножаются на одни и те же коэффициенты – значит, есть возможность группировки слагаемых перед умножением, что позволит снизить их количество. 172 x(n) -1/6 x(0) × z-1 × z-1 0.5 x(1) × z-1 + × z-1 + x(3) × z-1 0.5 + × a3 a2 × × z-1 × × + + × a1 + + 1 + 1/3 + × z-1 × z-1 + 0.5 + 0.5 z-1 × z-1 -1 × -1 + × z-1 1/6 × × z-1 -0.5 x(2) 1/6 + a0 + × + t y(t) Рис. 4.46 – Структурная схема фильтра Фарроу (N=4) Заметим, что коэффициенты: 1 1 1 1 1 1 a3 = − x(0 ) + x(1) − x(2 ) + x(3) = [x(3) − x(0 )] + [x(1) − x(2 )] 6 2 2 6 6 2 a1 + a3 = 1 [x(3) − x(1)] 2 a1 + a 2 + a3 = x(3) − x(2 ) → a1 = 1 [x(3) − x(1)] − a3 2 → a 2 = x(3) − x(2 ) − a1 − a3 a0 = x(2 ) Таким образом, модифицированный фильтр Фарроу содержит только 3 умножения для вычисления коэффициентов полинома еще 3 умножения для вычисления значения отсчета. Структурная схема модифицированного фильтра представлена на рис. 4.47. 173 x(0) x(n) Σ x(1) z-1 + Σ Σ - + Σ - - + Σ 0.5 0.5 + a0 + - Σ - a3 x(3) z-1 + 1/6 + x(2) z-1 - Σ a1 + a2 × + × + × + t y(t) Рис. 4.47 – Структурная схема модифицированного фильтра Фарроу (N=4) В заключение, необходимо рассмотреть какие значения может принимать переменная t. При синтезе фильтра рассматривались четыре момента времени: t0=-2; t1=-1; t2=0; t3=1, поэтому значение переменной t должно лежать внутри интервала [-2; 1]. Однако учитывая конечную степень интерполирующего полинома, максимальная точность интерполяции (минимальная погрешность) обеспечивается в интервале [-1; 0], т.е. в середине отрезка. Продолжая рассмотрение примера обработки сигнала с частотой дискретизации fs=44.1 кГц, изменим его частоту дискретизации до fsn=48 кГц. Поскольку требуется немедленное воспроизведение гипотетического контента, то данная задача относится к задачам реального времени, т.е. выполнение операций по обработке сигнала выполняется потоковым методом, при этом положим без потери общности, что сигнал монофонический. Обозначим: - шкала времени исходного сигнала: t x = nt s ; - шкала времени передискретизированного сигнала: t y = kt sn . 174 Введем для удобства обозначений коэффициент R = f s t sn = = 0.91875 . Тогда f sn t s для расчета отсчета выходного сигнала с произвольным индексом k понадобятся четыре отсчета исходного сигнала с индексами: n-2, n-1, n, n+1 (взяты для удобства нумерации в соответствии с моментами времени t0=-2; t1=-1; t2=0; t3=1). Рассчитать индекс n можно сопоставив шкалы времени сигнала, при этом ближайшим отсчетом к требуемому моменту времени kt sn будет являться отсчет с индексом n = kR  + 1 , где односторонние скобочки x  означают отбрасывание дробной части величины x. Таким образом, точке t0=-2 соответствует отсчет x(n-2), а точке t3=1 – отсчет x(n+1). В этом случае нетрудно пересчитать шкалу nts в интервал [-2; 1]: t = kR − n . Рассмотрим пример – запись звучания ноты «ля» (частота порядка 440 Гц) с частотой дискретизации fs=44.1 кГц (показан на рис. 4.48) требуется передискретизировать до частоты fsn=48 кГц. x(t) t, мс Рис. 4.48 – Затухающая синусоида с частотой 440 Гц (модель ноты «ля») В соответствии со спроектированным фильтром Фарроу для каждого отсчета выходного сигнала y(k) рассчитаем: − индексы отсчетов исходного сигнала n = kR  + 1 , т.е. из исходного сигнала делаем выборку из N=4 отсчетов: x(n-2), x(n-1), x(n), x(n+1); − переменную t = kR − n 175 − коэффициенты интерполятора a3 = a1 = 1 [x(3) − x(0)] + 1 [x(1) − x(2)], a0 = x(2) , 6 2 1 [x(3) − x(1)] − a3 , a2 = x(3) − x(2) − a1 − a3 , 2 − собственно интерполированные значения – исходный сигнал с другой частотой дискретизации y (k ) = [(a3t + a2 )t + a1 ]t + a0 Фрагмент результата работы фильтра Фарроу представлен на рис. 4.49 (зеленые кружочки – исходные отсчеты, красные крестики – отсчеты передискретизированного сигнала). y(k) t, мс Рис. 4.49 – Результат передискретизации (фильтром Фарроу N=4) 176 4.9 Передискретизация как задача регрессии данных С математической точки зрения задача передискретизации сигнала, как уже отмечалось, относится к задачам интерполяции функции по равноотстоящим точкам (узлам интерполяции). Однако реальные сигналы часто зашумлены, а ограниченный набор отсчетов сигнала не позволяет в полной мере говорить о статистических характеристиках шума. Поэтому интерполирование сигнала с малым чувствительности коэффициентов ОСШ приводит интерполирующего к большой многочлена от реализации шума (помех). Конечно, можно применить фильтрацию сигнала до передискретизации или после нее. Однако это требует дополнительных вычислительных затрат. Альтернативой может быть «подбор» такой функции, которая на участках максимально близка именно к полезной составляющей сигнала – это задача оптимальной аппроксимации. Это обеспечит процесс передискретизации и подавление шума одновременно. В этом контексте задачу можно сформулировать в следующем виде: минимизация ошибки между сигналом x(t) и аппроксимирующей функцией d(t):   p min ∫ [x(t ) − d (t )] dt  , где p – показатель нормы.   Если p=1 – мера наименьших модулей, например, метод Лагранжа дает хорошие результаты при обработке случайных сигналов с законом распределения близким к распределению Лапласа. Если p=2 – квадратичная мера обеспечивает максимальное правдоподобие функции приближения при нормальном законе распределения случайной составляющей зависимой переменной. При p=∞ – норма ошибки равна максимальному абсолютному отклонению сигнала от аппроксимирующей функции. Минимизация этой нормы соответствует минимаксной (равноволновой) аппроксимации (известна также как Чебышевская аппроксимация). 177 Аппроксимация данных с учетом их статистических параметров относится к задачам регрессии, которые обычно возникают при обработке экспериментальных данных, полученных в результате измерений статистических по своей природе сигналов. Как правило, при регрессионном анализе усреднение данных производится методом наименьших квадратов (p=2). Рассмотрим стационарный сигнал x(n), который содержит шум, распределенный, как правило, по нормальному закону. Требуется определить такую функцию d(n), часто в виде полинома: d (n ) = минимизирует выражение: ∑ [x(n ) − d (n )] 2 M −1 ∑ am n m , которая m =0 , обеспечивающее минимальную n погрешность приближения в среднеквадратическом смысле. Для определения коэффициентов полинома минимизируемое выражение дифференцируют относительно неизвестных величин am. Полученные частные производные приравнивают нулю и решают систему уравнений, из которой находят коэффициенты полинома (это соответствует отысканию экстремума функции). Виды регрессии обычно называются по типу аппроксимирующих функций: полиномиальная, экспоненциальная, логарифмическая и т.п. При применении регрессии большое значение имеет адекватность модели (вид, порядок) полезной составляющей сигнала. 4.9.1 Линейная регрессия Аппроксимирующий N −1 полином: d (n ) = an + b , минимизируется функционал: F = ∑ [x(n ) − an − b] , что приводит к системе уравнений: 2 n =0 N −1  ∂F = − 2 ∑ [x(n ) − an − b]n = 0  ∂a n =0  ∂F N −1  = −2 ∑ [x(n ) − an − b] = 0 n =0  ∂b Решение системы уравнений: 178 a= N −1 6  2 N −1  ( ) nx n − x(n ) ∑ ∑  N ( N + 1)  N − 1 n=0 n =0  b= 1 N N −1 ∑ x(n ) − n =0 N −1 a 2 4.9.2 Квадратичная регрессия Аппроксимирующий полином: d (n ) = an 2 + bn + c , [ минимизируется ] N −1 функционал: F = ∑ x(n ) − an 2 − bn − c , что приводит к системе уравнений: n =0 2 [ [ ] ] [ ] N −1  ∂F = − 2 ∑ x(n ) − an 2 − bn − c n 2 = 0  ∂a n =0  ∂F N −1  = − 2 ∑ x(n ) − an 2 − bn − c n = 0  ∂b n =0  N −1  ∂F = −2 x(n ) − an 2 − bn − c = 0 ∑  ∂c n =0 Решение системы уравнений: N −1 a = 30 N −1 N −1 n =0 n =0 6 ∑ n x(n ) − 6( N − 1)∑ nx(n ) − ( N − 1)( N − 2 )∑ x(n ) 2 n =0 (N − 2)(N − 1)N (N + 1)(N + 2) N −1 N −1 12 6 b= ∑ nx(n ) − N (N + 1) ∑ x(n ) − a(N − 1) N ( N − 1)( N + 1) n=0 n =0 c= 1 N N −1 ∑ x(n ) − a n =0 (N − 1)(2 N − 1) − b (N − 1) 6 2 4.9.3 Экспоненциальная регрессия Аппроксимирующий [ полином: d (n ) = Ae n , минимизируется ] N −1 функционал: F = ∑ x(n ) − Ae n , что приводит к системе уравнений: n =0 2 [ ] N −1  ∂F = − 2 ∑ x(n ) − Aen en = 0  ∂A n =0  ∂F N −1  = −2 ∑ x(n ) − Ae n Ae n n = 0  ∂ n =0 [ ] 179 Решение системы уравнений при x(n)>0: =6 N −1 N −1 n =0 n =0 2 ∑ n ln ( x(n )) − ( N − 1)∑ ln ( x(n )) 1 A = exp  N N ( N − 1)( N + 1) N −1 ∑ ln(x(n )) −  n =0 (N − 1) 2  4.9.4 Применение регрессии к зашумленному сигналу Рассмотрим несколько сигналов (при n<0 сигнал отсутствует, N=128): 0.5 + 2(1 − e n ), 0 ≤ n ≤ N − 1 x1 (n ) =  ,  = −0.0625 – сигнал имеет излом;  (n− N ) . 5 + 2 e , N ≤ n ≤ 2 N − 1  x 2 (n ) = 0.4 sin (n ) + 0.6 cos(n ),  = 2 , n ∈ [0, N − 1] – гладкий сигнал. N Оценим возможности задержки этих сигналов на половину периода дискретизации при воздействии стационарного гауссова шума (n) с параметрами: M=0, =0.4. Таким образом, по имеющемуся дискретному сигналу x(n) необходимо вычислить значение отсчета, расположенного посередине каждой пары исходного сигнала посредством следующих методов (для возможности сравнения выберем значения апертур M=4 и M=8):  регрессия линейная;  регрессия квадратичная;  регрессия экспоненциальная;  фильтрация по Фарроу (апертура M=4). Сопоставим результаты обработки сигнала с помощью различных регрессий и величин апертуры, а также результатов кубической регрессии и фильтрации по Фарроу 3-го порядка: y (k ) = [(a3t + a2 )t + a1 ]t + a0 , t=0.5 a3 = 1 [x(3) − x(0)] + 1 [x(1) − x(2)], a0 = x(2) , 6 2 180 a1 = 1 [x(3) − x(1)] − a3 , a2 = x(3) − x(2) − a1 − a3 2 Вычислим СКО разности обработанного сигнала и эталона в каждом случае. Результаты обработки сигнала x1 (n ) + (n ) :  рис. 4.50 – исходный сигнал и его зашумленная реализация ОСШ=5.49;  рис. 4.51 – обработка зашумленного сигнала линейной регрессией;  рис. 4.52 – обработка зашумленного сигнала квадратичной регрессией;  рис. 4.53 – обработка зашумленного сигнала экспоненциальной регрессией;  рис. 4.54 – обработка зашумленного сигнала фильтром Фарроу 3-го порядка. n Рис. 4.50 – Исходный сигнал и его зашумленная реализация (СКО=0.4) 181 n Рис. 4.51 – Применение линейной регрессии синий – апертура 4 (СКО=0.198) зеленый – апертура 8 (СКО=0.149) n Рис. 4.52 – Применение квадратичной регрессии синий – апертура 4 (СКО=0.307) зеленый – апертура 8 (СКО=0.209) 182 n Рис. 4.53 – Применение экспоненциальной регрессии синий – апертура 4 (СКО=0.206) зеленый – апертура 8 (СКО=0.159) n Рис. 4.54 – Применение фильтра Фарроу 3-го порядка (синий) СКО=0.416 183 Полученные данные сведены в таблицу 4.4, также в таблице отражены результаты расчета дисперсии  2spd разности спектров сигналов без шума в каждом случае и спектра исходного сигнала (для уменьшения размытия применено окно Хеннинга). Таблица 4.4 – Результаты расчета погрешностей -5 СКО  2spd ·10 № Описание 1. Исходный сигнал с шумом 2. 3. 4. 5. Линейная регрессия (M=4 | М=8) Квадратичная регрессия (M=4 | М =8) Экспоненциальная регрессия (M=4 | М =8) Фильтр Фарроу 3-го порядка 0.400 - 0.198 0.149 0.1009 0.5842 0.307 0.209 0.0683 0.3834 0.206 0.159 0.1062 0.6291 0.416 0.3552 Результаты обработки сигнала x 2 (n ) + (n ) :  рис. 4.55 – исходный сигнал и его зашумленная реализация ОСШ=1.63;  рис. 4.56 – обработка зашумленного сигнала линейной регрессией;  рис. 4.57 – обработка зашумленного сигнала квадратичной регрессией;  рис. 4.58 – обработка зашумленного сигнала экспоненциальной регрессией;  рис. 4.59 – обработка зашумленного сигнала фильтром Фарроу 3-го порядка. Полученные данные сведены в таблицу 4.5, также в таблице отражены результаты расчета дисперсии  2spd разности спектров сигналов без шума в каждом случае и спектра исходного сигнала (для уменьшения размытия применено окно Хеннинга). 184 Таблица 4.5 – Результаты расчета погрешностей -6 СКО  2spd ·10 № Описание 1. Исходный сигнал с шумом 2. 3. 4. 5. Линейная регрессия (M=4 | М=8) Квадратичная регрессия (M=4 | М =8) Экспоненциальная регрессия (M=4 | М =8) Фильтр Фарроу 3-го порядка 0.395 - 0.193 0.133 0.1218 0.7487 0.310 0.209 0.1067 0.5910 0.202 0.145 0.1182 0.7159 0.413 0.1073 n Рис. 4.55 – Исходный сигнал и его зашумленная реализация (СКО=0.4) 185 n Рис. 4.56 – Применение линейной регрессии синий – апертура 4 (СКО=0.1934) зеленый – апертура 8 (СКО=0.1331) n Рис. 4.57 – Применение квадратичной регрессии синий – апертура 4 (СКО=0.3100) зеленый – апертура 8 (СКО=0.2091) 186 n Рис. 4.58 – Применение экспоненциальной регрессии синий – апертура 4 (СКО=0.2017) зеленый – апертура 8 (СКО=0.1448) n Рис. 4.59 – Применение фильтра Фарроу 3-го порядка (синий) СКО=0.4125 187 В заключении, оценим вычислительную сложность обработки сигнала с помощью квадратичной регрессии с апертурой 4 и сравним ее со сложностью фильтрации по Фарроу 3-го порядка (полагая операции умножения и сложения одинаковыми по длительности выполнения, а операции с памятью – мгновенными). 1. Модифицированный фильтр Фарроу 3-го порядка требует 17 операций (11 – расчет коэффициентов полинома и 6 – на расчет отсчета) на каждый отсчет выходного сигнала. 2. Квадратичная регрессия с апертурой 4 требует 29 операций (25 – расчет коэффициентов полинома и 4 – на расчет отсчета) на каждый отсчет выходного сигнала (при расчетах полагаем величину апертуры известной, умножение на единицу и ноль, а также сложение с нулем не считаем). Таким образом, большая точность (по критериям СКО и  2spd ) решения задачи задержки зашумленного сигнала на половину периода дискретизации (независимо от гладкости сигнала) при применении квадратичной регрессии обеспечивается почти двукратными вычислительными затратами по сравнению с фильтрацией по Фарроу 3-го порядка. 188