ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ
Тема 2: ЧАСТОТНЫЙ АНАЛИЗ ЦИФРОВЫХ ФИЛЬТРОВ.
Не перестаю удивляться дерзкой гениальности Стефенсона и братьев Черепановых. Как они отважились построить паровоз, не располагая теорией его движения?
Архив Кифы Васильевича (Наука и жизнь, 1984).
Пока нет теории, есть возможность войти в Историю. Бог прославился созданием Евы из ребра Адама без всякого теоретического обоснования. А когда теория есть, можно только влипнуть в какую-нибудь историю.
Лариса Ратушная. Уральский геофизик (XX в.).
Содержание:
Введение.
2.1. Сглаживающие фильтры и фильтры аппроксимации. Фильтры МНК 1-го порядка. Фильтры МНК 2-го порядка. Фильтры МНК 4-го порядка.
2.2. Разностные операторы. Разностный оператор. Восстановление данных. Аппроксимация производных.
2.3. Интегрирование данных.
2.4. Расчет фильтра по частотной характеристике.
Литература.
Введение.
Основной инструмент цифровой фильтрации данных и проектирования цифровых фильтров – частотный анализ (второй распространенный термин – спектральный анализ). Частотный анализ базируется на использовании периодических функций, в отличие от численных методов анализа и математической статистики, где предпочтение отдается полиномам. В качестве периодических используются преимущественно гармонические функции – синусы и косинусы. По-существу, спектральный состав сигналов – это тонкая внутренняя структура данных, которые несет сигнал, и которая практически скрыта в динамическом (графическом) представлении больших множеств данных даже для опытных обработчиков. Точно так же частотная характеристика цифрового фильтра – это его однозначный функциональный паспорт, полностью определяющий сущность преобразования фильтром входных данных.
Однако следует отметить, что хотя сущность фильтрации сигналов состоит именно в направленном изменении частотного состава данных, которые несет сигнал, тем не менее, у начинающих специалистов существует определенное эмоциональное противодействие частотному подходу и его роли в анализе данных. Преодолеть это противодействие можно только одним путем – на опыте убедиться в эффективности частотного подхода.
Рассмотрим несколько примеров частотного анализа фильтров применительно к известным способам обработки данных.
2.1. Сглаживающие фильтры и фильтры аппроксимации /24/.
Предположим, что требуется осуществить сглаживание (регуляризацию, аппроксимацию) по методу наименьших квадратов (МНК) равномерного по аргументу массива данных.
Фильтры МНК 1-го порядка
(МНК-1). Простейший способ аппроксимации по МНК произвольной функции s(t) - с помощью полинома первой степени, т.е. функции вида y(t) = A+Bt (метод скользящих средних). В качестве примера произведем расчет симметричного фильтра на (2N+1) точек с окном от -N до N.
Для определения коэффициентов полинома найдем минимум функции приближения (функцию остаточных ошибок). С учетом дискретности данных по точкам tn
= nDt и принимая Dt = 1 для симметричного НЦФ с нумерацией отсчетов по n от центра окна фильтра (в системе координат фильтра), для функции остаточных ошибок имеем:
s(A,B) = [sn
- (A+B·n)]2
.
Дифференцируем функцию остаточных ошибок по аргументам 'А, В' и, приравнивая полученные уравнения нулю, формируем 2 нормальных уравнения:
(sn
-(A+B·n)) ºsn
- A1 - Bn = 0,
(sn
-(A+B·n))·n ºn×sn
- An - Bn2
= 0,
С учетом очевидного равенства n = 0, результат решения данных уравнений относительно значений А и В:
А = sn
, B =n×sn
/
n2
.
Подставляем значения коэффициентов в уравнение аппроксимирующего полинома, переходим в систему координат по точкам k массива y(k+t) = A+B·t, где отсчет t производится от точки k массива, против которой находится точка n = 0 фильтра, и получаем в общей форме уравнение фильтра аппроксимации:
y(k+t) = sk-n
+ tn×sk-n
/
n2
.
Для сглаживающего НЦФ вычисления производятся непосредственно для точки k в центре окна фильтра (t = 0), при этом:
yk
= sk-n
. (2.1.1)
|
Рис. 2.1.1. |
Импульсная реакция фильтра соответственно определяется (2N+1) значениями коэффициентов bn
= 1/(2N+1). Так, для 5-ти точечного НЦФ:
h(n) = {0.2, 0.2, 0.2, 0.2, 0.2}.
Передаточная функция фильтра в z-области:
H(z) = 0.2(z-2
+z-1
+1+z1
+z2
).
Коэффициент усиления дисперсии шумов:
Kq
= Sn
h2
(n) = 1/(2N+1),
т.е. обратно пропорционален ширине окна фильтра. Зависимость значения Kq
от ширины окна приведена на рис. 2.1.1.
Частотная характеристика фильтра
(передаточная функция фильтра в частотной области) находится преобразованием Фурье импульсной реакции h(n) (фильтр симметричный, начало координат в центре фильтра) или подстановкой z = exp(-jw) в выражение передаточной функции H(z). И в том, и в другом случае получаем:
H(w) = 0.2[exp(2jw)+exp(jw)+1+exp(-jw)+exp(-2jw)]. (2.1.2)
Можно использовать и непосредственно уравнение фильтра, в данном случае уравнение (2.1.1). Подадим на вход фильтра гармонический сигнал вида sk
= exp(jwk). Так как сигнальная функция относится к числу собственных, на выходе фильтра будем иметь сигнал yk
= H(w)exp(jwk). Подставляя выражения входного и выходного сигналов в уравнение (2.1.1), получаем:
H(w) exp(jwk) = 0.2exp(jw(k-n))= 0.2 exp(jwk) exp(-jwn).
Отсюда, выражение для передаточной функции:
H(w) = 0.2exp(-jwn) = 0.2[exp(2jw)+exp(jw)+1+exp(-jw)+exp(-2jw)],
что полностью идентично выражению (2.1.2).
Следует запомнить
: если оператор фильтра известен, то для получения его частотной характеристики достаточно подставить сигнал exp(jwn) непосредственно в линейное уравнение фильтра. Тем самым выполняются сразу 2 операции: производится z- преобразование h(n) и подставляется z = exp(-jwn), т.е. осуществляется трансформация h(n)→ h(z) → H(w).
Так как импульсная реакция фильтра МНК симметрична (функция h(n) четная), частотное представление передаточной функции должно быть вещественным, в чем нетрудно убедиться, объединив комплексно сопряженные члены выражения (2.1.2):
H(w) = 0.2(1+2 cos w+2 cos 2w).
Альтернативное представление передаточной функции H(w) для фильтра с произвольным количеством коэффициентов 2N+1 нам достаточно хорошо известно, как нормированный фурье-образ прямоугольной функции, каковой по существу и является селектирующее окно фильтра (2.1.1):
H(w) = sin((N+1/2)w)/[(N+1/2)w] = sinc((N+1/2)w). (2.1.3)
Рис. 2.1.2. Сглаживающие фильтры МНК.
Графики передаточных функций (2.1.3) приведены на рисунке 2.1.2. По графикам можно видеть коэффициент передачи сигнала с входа на выход фильтра на любой частоте. Без ослабления (с коэффициентом передачи 1) сглаживающим фильтром пропускается (и должен пропускаться по физическому смыслу сглаживания данных) только сигнал постоянного уровня (нулевой частоты). Этим же определяется и тот фактор (который стоит запомнить), что сумма коэффициентов сглаживающего НЦФ всегда должна быть равна 1 (отсчет ненормированного дискретного фурье-преобразования на частоте w = 0 равен сумме значений входной функции).
Чем больше число коэффициентов фильтра (шире окно фильтра), тем уже полоса пропускания низких частот. Подавление высоких частот довольно неравномерное, с осцилляциями передаточной функции относительно нуля. На рис. 2.1.3 приведен пример фильтрации случайного сигнала (шума) фильтрами с различным размером окна.
Рис. 2.1.3. Фильтрация шумов фильтрами МНК 1-го порядка.
Частотное представление передаточных функций позволяет наглядно видеть особенности фильтров и целенаправленно улучшать их характеристики. Так, если в рассмотренном нами фильтре с однородной импульсной реакцией hn
= 1/(2N+1) уменьшить два крайних члена в 2 раза и заново нормировать к сумме S hn
= 1, то частотные характеристики фильтра заметно улучшаются. Для нахождения передаточной функции модифицированного фильтра снимем в выражении (2.1.3) нормировку (умножим на 2N+1), вычтем значение 1/2 крайних членов (exp(-jwN)+exp(jwN))/2 = cos(wN) и заново пронормируем полученное выражение (разделим на 2N). Пример новой передаточной функции при N=3 также приведен на рисунке 2.1.2. Передаточные функции модифицированных таким образом фильтров приводятся к нулю на частоте Найквиста, при этом несколько расширяется полоса пропускания низких частот и уменьшается амплитуда осцилляций в области подавления высоких частот. Если смотреть на сглаживание, как на операцию подавления высокочастотных помех, то модифицированные фильтры без сомнения больше соответствует своему целевому назначению.
|
Рис. 2.1.4. |
При выборе окна фильтра следует учитывать как коэффициент подавления дисперсии шумов, так и степень искажения полезного сигнала, на который наложены шумы. Оптимальное окно фильтра может быть определено только в том случае, если спектр сигнала известен и ограничен определенной верхней частотой, а мощность шумов не превышает определенного уровня. Рассмотрим это на конкретном примере.
Допустим, что нужно обеспечить максимальное подавление дисперсии шумов при минимальном искажении верхней граничной частоты сигнала fв
, на которой мощность шумов равна мощности сигнальной гармоники fв
. Значение fв
равно 0.08 частоты Найквиста дискретизации данных, т.е. fв
= 0.04 при Dt=1. Относительные значения мощности (дисперсии) гармоники и шума принимаем равными 1. Спектр модели сигнала + шума в сопоставлении с передаточными функциями фильтров приведен на рис. 2.1.4.
Таблица 2.1.1.
| N |
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
| Ку
|
1 |
0.98 |
0.94 |
0.88 |
0.8 |
0.7 |
0.6 |
0.51 |
| Wu
|
1 |
0.96 |
0.88 |
0.77 |
0.64 |
0.51 |
0.38 |
0.26 |
| Wq
|
1 |
0.33 |
0.2 |
0.14 |
0.11 |
0.09 |
0.08 |
0.07 |
| Кс/ш
|
1 |
2.88 |
4.4 |
5.4 |
5.8 |
5.6 |
4.89 |
3.85 |
| d2
|
1 |
0.35 |
0.23 |
0.18 |
0.17 |
0.18 |
0.21 |
0.26 |
| s2
|
1 |
0.32 |
0.2 |
0.15 |
0.15 |
0.18 |
0.23 |
0.31 |
|
Рис. 2.1.5. |
По формуле (2.1.3) вычисляем коэффициенты Ку
(fв
) усиления фильтров с N от 0 до 6 на частоте fв
(см. таблицу 2.1.1). При мощности гармоники Wu
= 1 амплитудное значение гармоники на входе фильтра равно U = = 1.41. Мощности гармоник на выходе фильтров в зависимости от N:
Wu
(N)= 0.5·[U· Ку
(fв
)]2
.
Соответственно, при мощности входного шума Wq
=1 мощности шумов на выходе фильтров будут численно равны коэффициентам усиления дисперсии шумов Wq
(N) = Wq
·Kq
(N).
Максимум отношения
Кс/ш
(N) = Wq
(N)/Wu
(N)
определяет оптимальный фильтр с максимальным увеличением отношения сигнал/шум, т.е., по существу, коэффициент усиления отношения сигнал/шум при выполнении фильтрации с учетом изменения амплитудных значений полезной части сигнала.
|
Рис. 2.1.6. |
При Ку
(fв
) > 0.5 и Wu
(N) = Wq
(N) = 1 численные значения величины d2
(N) = 1/ Кс/ш
(N) в первом приближении могут служить оценкой s2
(N) квадрата среднего квадратического отклонения выходных сигналов от "чистой" гармоники fв
, заданной на входе. Свидетельством этому служат последние строки таблицы 2.1.1, где приведены результаты математического моделирования фильтрации по данным условиям на выборке 10000 точек. На рис. 2.1.6 приведены результаты сопоставления расчетных d2
(N) и модельных s2
(N) значений данных коэффициентов. Эффект фильтрации можно видеть на рис. 2.1.7, где приведен пример сигналов моделирования на ограниченном отрезке данных.
Рис. 2.1.7. Сигналы на входе и выходе фильтра МНК 1-го порядка.
Фильтры МНК 2-го порядка
(МНК-2) рассчитываются и анализируются аналогично. Рассмотрим квадратный многочлен вида y(t)=A+B·t+C·t2
. Для упрощения анализа ограничимся симметричным сглаживающим НЦФ с интервалом дискретизации данных Dt=1.
Минимум суммы квадратов остаточных ошибок:
s(A,B,C) = [sn
-(A+B·n+C·n2
)]2
. (2.1.4)
Система уравнений после дифференцирования выражения (2.1.4) по А, В, С и приравнивания полученных выражений нулю:
A1 + Bn + Сn2
=sn
.
An + Bn2
+ Сn3
=n·sn
.
An2
+ Bn3
+ Сn4
=n2
·sn
.
При вычислении значения квадратного многочлена только для центральной точки (t=0) необходимости в значениях коэффициентов В и С не имеется. Решая систему уравнений относительно А, получаем:
A = {n4
sn
-n2
n2
sn
} / {1n4
- [n2
]2
}. (2.1.5)
При развертывании выражения (2.1.5) для 5-ти точечного НЦФ:
yo
= (17sn
- 5n2
sn
) /35 = (-3·s-2
+12·s-1
+17·so
+12·s1
-3·s2
) /35. (2.1.6)
Импульсная реакция: hn
= {(-3, 12, 17, 12, -3)/35}.
Передаточная функция фильтра:
H(z)= (-3z-2
+12z-1
+17+12z1
-3z2
)/35. (2.1.7)
Рис. 2.1.8. Сглаживающие фильтры МНК.
Аналогичным образом выражение (2.1.5) позволяет получить импульсную реакцию для 7, 9, 11 и т.д. точек фильтра:
3
hn
= {(-2,3,6,7,6,3,-2)/21}.
4
hn
= {(-21,14,39,54,59,54,39,14,-21)/231}.
5
hn
={(-36,9,44,69,84,89,84,69,44,9,-21)/459}.
Подставляя значение z = exp(-jw) в (2.1.7) или непосредственно в (2.1.6) сигнал sn
= exp(jwn) и объединяя комплексно сопряженные члены, получаем частотную характеристику 5-ти точечного сглаживающего фильтра МНК второго порядка:
H(w) = (17+24 cos(w)-6 cos(2w))/35.
Вывод формул передаточных функций для 7, 9, 11-ти точечных фильтров МНК предлагается для самостоятельной работы.
Рис. 2.1.9. Рис. 2.1.10.
Вид частотных характеристик фильтров при N=3 и N=5 приводится на рис. 2.1.8. При сравнении характеристик с характеристиками фильтров МНК-1 можно видеть, что повышение степени полинома расширяет низкочастотную полосу пропускания фильтра и увеличивает крутизну ее среза. За счет расширения полосы пропускания главного частотного диапазона при тех же значениях N коэффициенты усиления дисперсии шумов фильтров МНК-2 выше, чем фильтров 1-го порядка, что можно видеть на рис. 2.1.9.
Методика выбора окна фильтра под частотные характеристики входных сигналов не отличается от фильтров МНК 1-го порядка. На рис. 2.1.10 приведены значения d2
(N) и s2
(N) фильтров МНК-2 в сопоставлении со значениями фильтров МНК-1 для частоты fв
= 0.08 Гц при Dt=1. Из сопоставления видно, что для получения примерно равных значений подавления шумов фильтры МНК-2 должны иметь в 2 раза большую ширину окна, чем фильтры МНК-1. Об этом же свидетельствует и пример моделирования фильтрации, приведенный на рис. 2.1.11.
Рис. 2.1.11.
Модификация фильтров.
Фильтры МНК второго порядка (равно как и другие фильтры подобного назначения) также можно модифицировать по условию H(w) → 0 при w → p. Один из простейших методов модификации заключается в следующем. В выражение передаточной функции (со всеми коэффициентами фильтра, вида (2.1.7)) подставляем z = exp(-jw), заменяем значения концевых коэффициентов фильтра на параметры, принимаем w = p, и, приравняв полученное выражение нулю, находим новые значения концевых коэффициентов, после чего сумму всех коэффициентов нормируем к 1 при w = 0.
Пример модификации фильтра МНК 2-го порядка.
Передаточная функция: выражение (2.1.7). Частотная характеристика (нормировку можно снять):
H(w) = -3exp(2jw)+12exp(jw)+17+12exp(-jw)-3exp(-2jw).
Замена концевых коэффициентов {значение 3} на параметр b и упрощение:
H(w) = 17+24 cos(w)+2b cos(2w).
При w = p: H(p) = 17-24+2b = 0. Отсюда: b = 3.5
Новая частотная характеристика (с приведением коэффициентов к целым числам):
H(w) = 68+96 cos(w)+14 cos(2w). Сумма коэффициентов при w = 0: H(0) = 68+96+14 = 178.
Нормированная частотная характеристика: H(w) = (68+96 cos(w)+14 cos(2w))/178.
Коэффициенты фильтра: hn
= {(7,48,68,48,7)/178}.
Пример-
задание: Модифицировать 7, 9 и 11-ти точечные сглаживающие фильтры МНК 2-го порядка.
Контроль: 7
hn
= {(1,6,12,14,12,6,1)/52}. 9
hn
= {(-1,28,78,108,118,108,78,28,-1)/548}.
11
h n
= {(-11,18,88,138,168,178,168,138,88,18,-11)/980}.
Сравнительные графики частотных характеристик модифицированных фильтров МНК второго порядка приведены на рисунке 2.1.8.
Фильтры МНК третьего порядка по своим частотным характеристикам эквивалентны фильтрам второго порядка.
Фильтры МНК 4-го порядка.
Расчет по аналогичной методике сглаживающих фильтров МНК 4-ой степени дает следующие результаты:
h0-3
= (131,75,-30,5)/231,
h0-4
= (179,135,30,-55,15)/429,
h0-5
= (143,120,60,-10,-45,18)/429,
h0-6
= (677,600,390,110,-135,-198,110)/2431.
На рис. 2.1.12 приведено сопоставление частотных характеристик одноразмерных фильтров МНК 1-го, 2-го и 4-го порядка.
Рис. 2.1.12. Сглаживающие фильтры МНК.
В целом, по сглаживающим фильтрам МНК можно сделать следующие выводы:
1. Повышение порядка фильтра увеличивает степень касания частотной характеристикой уровня коэффициента передачи Н=1 на частоте w = 0 и расширяет полосу пропускания фильтра.
2. Увеличение количества членов фильтра приводит к сужению полосы пропускания и увеличивает крутизну ее среза.
3. Модификация фильтров уменьшает осцилляции передаточной функции в полосе подавления сигналов.
Совместное изменение этих параметров позволяет подбирать для сглаживания данных такой фильтр МНК, частотная характеристика которого наилучшим образом удовлетворяет частотному спектру сигналов при минимальном количестве коэффициентов фильтра.
2.2. Разностные опера
Рассмотрим примеры частотного подхода при анализе разностных операторов.
Разностный оператор
1-го порядка имеет вид:
Dsk
= sk+1
-sk
.
Последовательное n-кратное применение оператора записывается в виде оператора n-го порядка:
Dn
(sk
) = D[Dn-1
(sk
)] = Dsk
*
Dn-1
(sk
) (2.2.1)
| k
|
sk
|
D
|
D
|
D
|
D
|
D
|
D
|
| -7 -6 -5 -4 -3 -2 -1 0 1 |
0 0 0 0 0 0 0 1 0 |
0 0 0 0 0 0 1 -1 0 |
0 0 0 0 0 1 -2 1 0 |
0 0 0 0 1 -3 3 -1 0 |
0 0 0 1 -4 6 -4 1 0 |
0 0 1 -5 10 -10 5 -1 0 |
0 1 -6 15 -20 15 -6 1 0 |
| Кq
|
2 |
6 |
20 |
70 |
252 |
924 |
Слева приводится таблица выходных значений импульсной реакции разностных операторов на единичный импульсный сигнал Кронекера. Как видно из таблицы, ряды последовательных разностей содержат знакопеременные биномиальные коэффициенты. В представленной форме разностные операторы являются каузальными фазосдвигающими (односторонними), но нетрудно заметить, что операторы четных степеней могут быть переведены в симметричную форму сдвигом вперед на половину окна оператора.
В последней строке таблицы приводятся коэффициенты усиления дисперсии шумов, значение которых резко нарастает по мере увеличения порядка оператора. Это позволяет использовать разностные операторы с порядком выше 1 для определения местоположения статистически распределенных шумов в массивах данных. Особенно наглядно эту возможность можно видеть на частотных характеристиках операторов.
Подставляя сигнал s(k) = exp(jwk) в (2.2.1) и упрощая, получаем:
Dn
s(k) = (jn
) exp(jwn/2) [2 sin(w/2)]n
exp(jwk). (2.2.2)
Так как первые два множителя в выражении (2.2.2) равны 1, зависимость коэффициента передачи разностного оператора от частоты определяется вторым сомножителем (2 sin(w/2))n
и представлена на рисунке 2.2.1.
Рис. 2.2.1. Разностные фильтры.
Как следует из рисунка, разностные операторы подавляют постоянную составляющую сигнала и его гармоники в первой трети интервала Найквиста и увеличивают высокочастотные составляющие сигнала в остальной части интервала тем больше, чем больше порядок оператора. Как правило, эту часть главного интервала спектра сигналов занимают высокочастотные статистические шумы.
Шумы при анализе данных также могут представлять собой определенную информацию, например, по стабильности условий измерений и по влиянию на измерения внешних дестабилизирующих факторов. На рис. 2.2.2 приведен пример выделения интервалов интенсивных шумов в данных акустического каротажа, что может свидетельствовать о сильной трещиноватости пород на этих интервалах. Такая информация относится уже не шумовой, а к весьма полезной информации при поисках и разведке нефти, газа и воды.
Рис. 2.2.2.
Восстановление данных.
Разностные операторы имеют одну особенность: оператор n+1 порядка аннулирует полином степени n, т.е. свертка оператора n+1 порядка с полиномом n-ой степени дает нулевые значения: Dn+1
*
Pn
(k) = 0. Эту особенность можно использовать для создания очень простых и достаточно надежных операторов восстановления в массивах пропущенных и утраченных значений или для замены аннулированных при обработке величин (например, явных выбросов).
Пример.
P2
(k) = xk
= 1+2k-k2
, k = 0,1,2,... xk
= 1,2,1,-2,-7,-14,-23,-34,... yk
= xk*
D3
=0,0,0,0,...
Если считать, что отрезок данных, содержащий пропуск, является многочленом некоторой степени, то свертка данных с разностным оператором следующего порядка должна быть равна нулю. Так, при аппроксимации данных многочленом третьей степени для любой точки массива должно выполняться равенство:
D4
·(sk
) = sk-2
-4sk-1
+6sk
-4sk+1
+sk+2
= 0.
Интерполяционный фильтр восстановления утраченной центральной точки данных:
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. Эти значения показывают, что операция двойного дифференцирования может применяться только для данных, достаточно хорошо очищенных от шумов, с основной энергией сигнала в первой трети интервала Найквиста.
В принципе, вторую производную можно получать и последовательным двойным дифференцированием данных оператором первой производной. Однако для таких простых операторов эти две операции не тождественны. Оператор последовательного двойного дифференцирования можно получить сверткой оператора первой производной с самим собой:
2
h1 = h1*
h1 = {0.25, 0, -0.5, 0, 0.25},
и имеет коэффициент усиления дисперсии шумов всего 0.375. Частотная характеристика оператора:
2
H1(w) = 0.5[1-cos(2w)].
Графики 2
H1(w) и коэффициента соответствия 2
K1(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).
При численном интегрировании по формуле Симпсона уравнение фильтра имеет вид:
yk+1
= yk-1
+(sk+1
+4sk
+sk-1
)/3. (2.3.4)
Частотный анализ фильтра проведите самостоятельно. Контроль:
K(w) = (2+cos w)/[3 sin(w)/w].
|
Рис. 2.3.1. Коэффициенты соответствия. |
Графики функций К(w) приведены на рисунке 2.3.1. При интегрировании происходит накопление результатов по всему предыдущему циклу суммирования и в этих условиях значение коэффициента K(w) является более представительным и информационным, чем передаточная функция оператора для одной текущей точки.
Наиболее простые формулы цифрового интегрирования, трапеций и прямоугольников, ведут себя различным образом в главном частотном диапазоне. Формула прямоугольников завышает результаты на высоких частотах, а формула трапеций - занижает. Эти особенности легко объяснимы. Для одиночной гармоники площадь трапеции по двум последовательным отсчетам всегда меньше, чем площадь с выпуклой дугой гармоники между этими отсчетами, и разница тем больше, чем больше частота. В пределе, для гармоники с частотой Найквиста, отсчеты соответствуют знакочередующемуся ряду (типа 1, -1, 1, -1, ... или любые другие значения в зависимости от амплитуды и начального фазового угла) и при нулевых начальных условиях суммирование двух последовательных отсчетов в формуле (3.2.1) будет давать 0 и накопления результатов не происходит. Интегрирование по площади прямоугольников с отчетом высоты по центральной точке между двумя отсчетами всегда ведет к завышению площади прямоугольника относительно площади, ограниченной выпуклой дугой гармоники.
Формула Симпсона отличается от формул трапеций и прямоугольников более высокой степенью касания единичного значения, что обеспечивает более высокую точность интегрирования в первой половине главного диапазона. Однако на высоких частотах погрешность начинает резко нарастать вплоть до выхода на бесконечность на конце диапазона (полюс в знаменателе передаточной функции рекурсивного фильтра на частоте Найквиста).
Эти особенности интегрирования следует учитывать при обработке данных сложного спектрального состава.
Курсовая работа 4 - Разработка методики расчета полосовых фильтров интегрирования.
2.4. Расчет фильтра по частотной характеристике.
В качестве примера проведем расчет простого симметричного сглаживающего НЦФ исходя непосредственно из требуемой формы частотной характеристики. Расчет выполним для фильтра с окном в пять точек:
yk
= ask-2
+bsk-1
+csk
+bsk+1
+ask+2
. (2.4.1)
Полагаем sk
= exp(jwk), при этом yk
= H(w)exp(jwk). Подставляем значения входного и выходного сигнала в уравнение фильтра, сокращаем левую и правую части на общий член exp(jwk) и, объединяя комплексно сопряженные члены в правой части, получаем уравнение передаточной функции:
H(w) = 2a cos(2w)+2b cos(w)+ c.
Сокращаем количество параметров функции заданием граничных условий по частоте. Как правило, имеет смысл принять: H(0) = 1, H(p) = 0. Отсюда:
H(0) = 2a+2b+c = 1,
H(p) = 2a-2b+c = 0.
B = 1/4, c = 1/2-2a.
При этом функция H(w) превращается в однопараметровую:
H(w) = 2a(cos(2w)-1)+(cos(w)+1)/2.
По полученному выражению рекомендуется построить семейство кривых в параметрической зависимости от значений 'а' и выбрать фильтр, удовлетворяющий заданию. Пример семейства частотных характеристик приведен на рисунке 2.4.1.
Рис. 2.4.1. Частотные характеристики НЦФ.
Можно наложить еще одно дополнительное условие и определить все коэффициенты фильтра непосредственно. Так, например, если к двум граничным условиям задать третье условие сбалансированности: H(w=p/2) = 0.5, то из трех полученных уравнений сразу же получим все три коэффициента фильтра: a = 0, b = 1/4, c = 1/2 (фильтр сокращается до трех точек).
В принципе, таким методом можно задать любую произвольную форму частотной характеристики симметричного НЦФ с произвольным количеством N точек дискретизации, что определит полное уравнение (2.4.1) с окном 2N+1 точка и соответствующую передаточную функцию фильтра, по которой можно составить и решить N+1 уравнение для определения коэффициентов фильтра.
литература
24. Хемминг Р.В. Цифровые фильтры. – М.: Недра, 1987. – 221 с.
Главный сайт автора
¨
Лекции по ЦОС
¨
Практикум
Сайт автора-
http://prodav.narod.ru/
О замеченных опечатках, ошибках и предложениях по дополнению
:
davpro
@
yandex
.
ru
.
Copyright ©2005 Davydov А.V.