Тема 2: частотный анализ цифровых фильтров (стр. 4 из 5)
sk = (-sk-2+4sk-1+4sk+1-sk+2)/6. (2.2.3)
Соответственно, оператор фильтра восстановления данных h(n) = (-1,4,0,4,-1)/6. Коэффициент усиления шумов s2 = 17/18 = 0.944.
Пример. Фактический отрезок массива данных: xk = {3,6,8,8,7,5,3,1}.
Допустим, что на отрезке был зарегистрирован явный выброс: xk = {3,6,8,208,7,5,3,1}.
Отсчет с выбросом аннулирован. Замена отсчета: x3 = (-x1+4x2+4x4-x5)/6= (-6+32+28-5)/6 » 8.17.
В массиве утрачен 5-й отсчет. Восстановление: x4 = (-x2+4x3+4x5-x6)/6 = (-8+32+20-3)/6 » 6.83.
Принимая в (2.2.3) k = 0 и подставляя сигнал sk = exp(jwk), получаем частотную характеристику, в данном случае - интерполяционного фильтра 4-го порядка:
H(w) = (4 cos w - cos 2w)/3.
Рис. 2.2.3. Разностные фильтры.
Вид частотной характеристики для фильтров восстановления пропущенных данных 4-го и 6-го порядков приведен на рис. 2.2.3. Графики наглядно показывают, что применение разностных интерполяционных фильтров восстановления данных возможно только для сигналов, высокочастотные и шумовые составляющие которых минимум в три раза меньше частоты Найквиста. Интерполяционные фильтры выше 4-го порядка применять не рекомендуется, т.к. они имеют коэффициент усиления шумов более 1.
На рис. 2.2.4 – 2.2.6 приведены примеры восстановления утраченных данных во входных сигналах оператором 3-го порядка и спектры сигналов в сопоставлении с передаточной функцией оператора восстановления данных. В сигналах утрачен каждый 10-ый отсчет (например, при передаче данных) при сохранении тактовой частоты нумерации данных. Учитывая, что все значения входных сигналов положительны, индикатором пропуска данных для работы оператора служат нулевые значения. В любых других случаях для оператора восстановления данных необходимо предусматривать специальный маркер (например, заменять аннулированные данные или выбросы определенным большим или малым значением за пределами значений отсчетов).
Рис. 2.2.4. Восстановление незашумленных данных. Рис.2.2.5. Спектры.
Рис. 2.2.6. Восстановление зашумленных данных.
Как следует из рис. 2.2.5, спектр полезного сигнала полностью находится в зоне единичного коэффициента частотной характеристики оператора, и восстановление данных выполняется практически без погрешности (рис. 2.2.4). При наложении на сигнал статистически распределенных шумов (рис. 2.2.6) погрешность восстановления данных увеличивается, но для информационной части полного сигнала она, как и во входных данных, она не превышает среднеквадратического значения (стандарта) флюктуаций шума. Об этом свидетельствует рис. 2.2.7, полученный для сигналов на рис. 2.2.6 по данным математического моделирования при разных значениях стандарта шума (выборки по 10 точкам восстановления).
Рис. 2.2.7. Погрешности восстановления сигналов.
Аппроксимация производных - вторая большая область применения разностных операторов. Оценки первой, второй и третьей производной можно производить по простейшим формулам дифференцирования:
(sn)' = (sn+1-sn-1)/2Dt. h1 = {-0.5, 0, 0.5}. (2.2.4)
(sn)'' = (sn+1-2sn+sn-1)/Dt. h2 = {1, -2, 1}.
(sn)''' = (-sn+2+2sn+1-2sn-1+sn-2)/2Dt. h3 = {0.5, -1, 0, 1, -0.5}.
Оператор первой производной является нечетной функцией и имеет мнимый спектр. Если принять s(t) = exp(jwt), то истинное значение первой производной должно быть равно: s'(t) = jw exp(jwt). Передаточная функция H(w) = jw. Оценка первой производной в точке n = 0 по разностному оператору при Dt = 1: s'(0) = (exp(jw)-exp(-jw))/2 = j sin w = H1(w). Отношение расчетного значения к истинному на той же точке: K1(w) = sin(w)/w. Графики функций в правой половине главного диапазона приведены на рис. 2.2.8.
 Рис. 2.2.8. |
Как следует из приведенных выражений и графиков, значение К(w) равно 1 только на частоте w = 0. На всех других частотах в интервале Найквиста формула дает заниженные значения производных. Однако при обработке практических данных последний фактор может играть и положительную роль, если сигнал низкочастотный (не более 1/3 главного диапазона) и зарегистрирован на уровне высокочастотных шумов. Любое дифференцирование поднимает в спектре сигнала долю его высокочастотных составляющих. Коэффициент усиления дисперсии шумов разностным оператором дифференцирования непосредственно по его спектру в главном диапазоне:
Kq = (1/p)
(sin w)2 dw = 0.5При точном дифференцировании по всему главному диапазону:
Kq = (1/p)
w2 dw = 3.29Следовательно, разностный оператор имеет практически в шесть раз меньший коэффициент усиления дисперсии шумов, чем полный по главному диапазону точный оператор дифференцирования.
На рис. 2.2.9 показан пример дифференцирования гармоники с частотой 0.1 частоты Найквиста (показана пунктиром) и этой же гармоники с наложенными шумами (сплошная тонкая кривая).
Рис. 2.2.9. Пример дифференцирования (входные сигналы – вверху, выходные – внизу).
Рис. 2.2.10. Частотные функции 2-ой производной.
Оператор второй производной относится к типу четных функций. Частотная функция оператора: H2(w) = -2(1-cos w). Собственное значение операции H(w) = -w2. Отношение фактического значения к собственному
K2(w) = [sin(w/2)/(w/2)]2
и также равно 1 только на частоте w = 0. На всех других частотах в интервале Найквиста формула дает заниженные значения производных, хотя и меньшие по относительным значениям, чем оператор первой производной. Частотные графики функций приведены на рис. 2.2.10. Коэффициент усиления дисперсии шумов оператором второй производной равен 6 при собственном значении дифференцирования, равном 19.5. Эти значения показывают, что операция двойного дифференцирования может применяться только для данных, достаточно хорошо очищенных от шумов, с основной энергией сигнала в первой трети интервала Найквиста.
В принципе, вторую производную можно получать и последовательным двойным дифференцированием данных оператором первой производной. Однако для таких простых операторов эти две операции не тождественны. Оператор последовательного двойного дифференцирования можно получить сверткой оператора первой производной с самим собой:
2h1 = h1* h1 = {0.25, 0, -0.5, 0, 0.25},
и имеет коэффициент усиления дисперсии шумов всего 0.375. Частотная характеристика оператора:
2H1(w) = 0.5[1-cos(2w)].
Графики 2H1(w) и коэффициента соответствия 2K1(w) приведены пунктиром на рис. 2.2.10. Из их сопоставления с графиками второй производной можно видеть, что последовательное двойное дифференцирование возможно только для данных, спектральный состав которых занимает не более пятой начальной части главного диапазона.
Рис. 2.2.11. Вторая производная гармоники с частотой w=0.2p при Dt=1
(пунктир – двойное последовательное дифференцирование)
Пример применения двух операторов второй производной приведен на рис. 2.2.11.
Попутно заметим, что частота Найквиста главного диапазона обратно пропорциональна интервалу Dt дискретизации данных (wN = p/Dt), а, следовательно, интервал дискретизации данных для корректного использования простых операторов дифференцирования должен быть в 3-5 раз меньше оптимального для сигналов с известными предельными частотами спектрального состава.
Частотные функции для третьей производной предлагается получить самостоятельно.
Курсовая работа 1 - Разработка простых операторов дифференцирования и методики их расчета.
Курсовая работа 2 - Разработка простых операторов второй производной и методики их расчета.
Курсовая работа 3 - Разработка простых операторов третьей производной и методики их расчета.
2.3. Интегрирование данных /24/
Интегрирование сигналов реализуется рекурсивными цифровыми фильтрами. Рассмотрим примеры анализа интегрирующих операторов.
Алгоритм интегрирования по формуле трапеций при нулевых начальных условиях:
yk+1 = yk+(sk+1+sk)/2. (2.3.1)
Принимая sk = exp(jwt) и yk = H(w)exp(jwt), подставляем сигналы в (2.3.1) при tk = kDt, Dt = 1 и решаем относительно H(w). Получаем:
H(w) = (exp(jw)+1)/[2(exp(jw)-1)].
H(w) = cos(w/2)/[2j sin(w/2)].
Истинное значение интеграла равно (1/jw)exp(jwt). Отсюда:
K(w) = H(w)exp(jwt)/[(1/jw)exp(jwt)].
K(w) = cos(w/2)[(w/2)/sin(w/2)]. (2.3.2)
Интегрирование по формуле прямоугольников (интерполяционное среднеточечное). Оператор:
yk+1 = yk+sk+1/2. (2.3.3)
После аналогичных подстановок сигнала и преобразований получаем:
K(w) = (w/2)/sin(w/2).
При численном интегрировании по формуле Симпсона уравнение фильтра имеет вид: