Использование фильтров Винера. Результаты восстановления смазанных изображений автомобильных номеров
Другая постановка задачи при расчете по критерию минимума дисперсии ошибки заключается в том, что ставится вопрос о нахождении оптимальной структуры и значений параметров системы автоматического управления, при которых обеспечивается получение теоретического минимума среднеквадратичной ошибки при заданных вероятностных характеристиках полезного сигнала (задающего воздействия) и помехи Эта задача будет решена, если найти, например, частотную передаточную функцию замкнутой системы для непрерывной системы или для дискретной системы (рис. 4.1). Передаточной функции может быть поставлена в соответствие передаточная функция а передаточной функции - функция Задача относится к категории вариационных задач в открытой области, т. е. без ограничений на фазовые координаты системы и управляющие воздействия.
Для решения этой задачи требуется знание статистических характеристик полезного входного сигнала (задающего воздействия) и помехи на входе системы (или помехи, пересчитанной на вход системы). При этом предполагается, что на систему управления заранее не налагается никаких ограничений в смысле обязательного использования реальных элементов (чувствительных элементов, усилителей, исполнительных элементов и др.) с заданными Характеристиками.
Запишем критерий оптимальности в задаче Н. Винера. При поступлении на вход системы аддитивной смеси полезного сигнала и помехи (рис. 4.1)
представляющих собой стационарные случайные функции с нулевыми математическими ожиданиями и известными
корреляционными функциями, требуется найти частотную передаточную функцию замкнутой системы или ей соответствующую физически реализуемую весовую функцию осуществляющую требуемое линейное преобразование входного сигнала
где - заданный линейный оператор, и обеспечивающую минимум дисперсии ошибки (4.4):
Если то это будет задача оптимального сглаживания, т. е. выделения сигнала из аддитивной смеси полезного сигнала и помехи. При равенстве помехи нулю решение задачи сглаживания имеет тривиальный вид:
Рис. 4.1. Оптимальный фильтр Винера
В задачах оптимального статистического упреждения где - время упреждения. Решение получается нетривиальным даже в случае отсутствия помехи. В задачах дифференцирования сигнала при наличии помех заданный линейный оператор имеет вид Но где - порядок отыскиваемой производной.
После нахождения оптимальной передаточной функции или конструктор должен попытаться реализовать ее посредством использования тех элементов, которыми он располагает и из которых должна быть построена система управления. Так как в большинстве практических случаев точное воспроизведение оптимальной передаточной функции оказывается невозможным, то приходится использовать квазиоптимальную, или субоптимальную, систему, более или менее близко совпадающую по своим параметрам с оптимальной.
Концепции оптимального линейного оценивания являются фундаментальными при любом рассмотрении адаптивных фильтров. Процесс адаптивной фильтрации включает два этапа проведения оценивания: 1) оценивание искомого выхода фильтра и 2) оценивание весов фильтра, необходимых для достижения вышеупомянутой цели. Второй из этих двух этапов необходим вследствие того, что в случае адаптивной фильтрации характеристики входного сигнала априорно не известны.
Наиболее широко распространенным типом структуры адаптивного фильтра является структура, в которой используется архитектура с конечной импульсной характеристикой (КИХ). Эти фильтры должны сходится к решению с помощью оптимального нерекурсивного устройства оценки, причем решение задается уравнением Винера – Хопфа.
Синтез КИХ и БИХ устройств оценки существенно зависит от определения стоимостной функции, в соответствии с которым качество оценивания характеризуется разностью между выходным сигналом устройства оценки и истинным параметром, подлежащим оцениванию:
Здесь e(n) – ошибка оценивания; x(n) – случайная величина, которую необходимо оценить и которая может быть детерминированной, а – оценка , выполненная с помощью нашей системы оценивания, причем
т.е. x(n) – линейная функция последовательности входных сигналов y(n) и набора весов фильтра h(n) . Наблюдаемую последовательность сигналов y(n) в общем виде можно представить как исходную последовательность x(n) , искаженную адаптивным белым шумом v(n) с дисперсией σ v 2:
. (5.26)
Наиболее употребительным при проведении оптимального оценивания является метод наименьших квадратов (МНК). Среднеквадратическая ошибка определяется как
Она минимизируется относительно весовых коэффициентов устройства оценки для получения оптимального оценивания по критерию МНК. Следует отметить, что можно применять не только описанную функцию стоимости. Альтернативными будут такие функции, как абсолютная величина ошибки и нелинейная пороговая функция. Такая функция ошибки используется в том случае, когда имеется приемлемый интервал ошибок (т.е. существует заданная допустимая ошибка). При использовании критерия наименьшего среднеквадратичного малые ошибки вносят меньший вклад, чем большие ошибки (в противоположность критерию абсолютной величины ошибки, который дает одинаковый вес для всех ошибок).
Рис. 5.9. Обобщенный нерекурсивный фильтр или устройство оценки.
В нерекурсивном устройстве оценки оценка x(n) определяется в виде конечного линейного полинома y(n) :
, (5.28)
где h k – отдельные веса в структуре нерекурсивного фильтра КИХ-типа, показанного на рис. 5.9. Выражение (5.28) можно переписать в матрично-векторной системе обозначений:
И 
,
а верхний индекс Т обозначает транспонирование матрицы. Тогда функция среднеквадратичной ошибки принимает вид
Это выражение описывает стандартную поверхность квадратичной ошибки с одним единственным минимумом. Дифференцирование (5.30) по дает
. (5.31)
а допуская, что (5.31) равно нулю, имеем
(5.32)
Полагая, что весовой вектор и вектор сигнала Y(n) не коррелированы, получаем
Члены математического ожидания, входящие в (5.33), можно определить следующим образом:
P=E{x(n)Y(n)} – взаимная корреляция между входным сигналом и оцениваемым параметром;
R=E{Y(n)Y T (n)} – автокорреляционая матрица входной сигнальной последовательности.
Тогда (5.33) можно переписать в виде
P T =H T opt R. (5.34)
Уравнение (5.34) является общеизвестным уравнением Винера – Хопфа, которое дает оптимальное (по методу наименьших квадратов) винеровское решение для H.
Здесь рассмотрим частный случай задачи, разобранный в предыдущем параграфе, и назовем его фильтром Винера. Для простоты ограничимся исследованием только непрерывного случая. Все выкладки, представленные в этом параграфе, являются частным случаем общей теории фильтра Калмана-Винера. Метод вычислений, основанный на алгоритме фильтра Калмама, вообще говоря, ближе к практическому осуществлению, чем метод, рассматриваемый здесь. С другой стороны, многие практически важные задачи оценивания можно отнести, по крайней мере с достаточным приближением, к стационарным, и методы, излагаемые в этом параграфе и разработанные раньше общей теории Калмана, уже нашли успешное применение в многочисленных практических задачах.
В 1949 г. была опубликована работа Винера «Экстраполяция, интерполяция и сглаживание стационарных временных последовательностей». Ее публикация явилась важной вехой не только потому, что результаты были новыми и вызвали к себе повышенный интерес, но и (что более важно) они возводили частную задачу в ранг теории, которая в то время нашла широкое применение, в частности, теории частотных фильтров. К сожалению, из-за того, что основные результаты были сформулированы на «частотном языке», они непосредственно не могли быть обобщены на нестационарные задачи. Хотя нестационарная задача и была сформулирована в общем виде, т. е. записано уравнение Винера-Хопфа, но было получено очень мало практических результатов, за исключением только работ Бутона , Заде и Рагазини , . До тех пор, пока не был разработан алгоритм фильтра Калмана, вычислительные трудности не были преодолены в общем нестационарном случае.
Стационарный фильтр Калмана. В стационарном варианте общей задачи оценивания состояния должны выполняться следующие три условия:
1. Модели сообщения и наблюдения не изменяются во времени, т. е. они описываются уравнениями с постоянными коэффициентами:
, (7.152)
где - матрицы постоянных коэффициентов.
2. Входной шум и шум измерений стационарны, по крайней мере, в широком смысле, т.е.
где и - матрицы постоянных коэффициентов. Кроме того, предполагается, что и - некоррелированные белые шумы с нулевым средним.
3. Интервал наблюдений начинается при . Очевидно, что это условие никогда не может выполняться на практике. Однако, поскольку момент начала наблюдений расположен достаточно далеко в прошлом, так что все переходные процессы успевают закончиться, то это допущение можно считать справедливым.
Если эти три допущения выполняются, то, очевидно,
задача оценивания уже не зависит от выбора начала отсчета времени в том смысле,
что допустимо любое конечное перемещение шкалы времени без изменения условий
задачи. Поэтому коэффициент усиления фильтра Калмана должен быть постоянным для
конечных значений времени , поскольку не существует причины его
изменения между двумя конечными моментами времени . Кроме того,
случайный процесс и,
следовательно, процесс стационарны,
так что
, а
уравнение дисперсии (7.105) записывается в виде
Постоянный коэффициент усиления Калмана в этом случае определяется как
. (7.156)
и, наконец, уравнение фильтрации имеет вид
Заметим, что в стационарном случае уравнение дисперсии превращается в вырожденное матричное уравнение Риккати.
Один из часто используемых способов решения ур-ния (7.155) (обычно с помощью ЦВМ) заключается в решении нестационарного уравнения дисперсии (7.105) с соответствующими постоянными значениями коэффициентов, из которых составлены матрицы и , и произвольной неотрицательно определенной матрицей начальных условий для в текущем времени до тех пор, пока полученное решение не достигнет постоянного установившегося значения. Это окончательное значение принимается за искомое решение ур-ния (7.155). Здесь алгебраическое уравнение преобразуется в дифференциальное, так как алгоритмы решения дифференциальных уравнений на цифровых (или аналоговых) вычислительных машинах, как правило, эффективнее алгоритмов решения нелинейных алгебраических уравнений. Другой подход, который может быть использован для отыскания решения ур-ния (7.155), связан с реализацией на ЦВМ поисковой процедуры, например, градиентного метода .
Пример 7.7. Определим стационарный фильтр, обеспечивающий минимальною дисперсию ошибки, для системы
Вырожденное уравнение Риккати (7.155) для этого примера записывается в виде
Если перемножить все указанные матрицы, то получим следующие три уравнения
Решение последнего уравнения относительно имеет вид . Если выбрать в качестве решения положительный корень этого уравнения, то получим и соответственно . Для того чтобы была положительно определенной, необходимо выбрать положительные значения для и - так, что в результате имеем
.
Таким образом, мы нашли, как и требовалось, действительные положительно определенные корни вырожденного уравнения Риккати.
Чтобы найти постоянный коэффициент усиления фильтра Калмана , достаточно подставить найденную матрицу в ур-ние (7.156). В результате получаем
.
На рис. 7.9, а показана «каноническая» реализация фильтра Винера. Если нас интересуют, как это часто бывает, только оценки состояния или , то могут быть использованы фильтры, реализованные по структурным схемам рис. 7.9б, 7.9в. Эти фильтры могут оказаться проще фильтра, изображенного на рис. 7.9а. Однако фильтр, структурная схема которого представлена на рис. 7.9а, одновременно формирует оценки и , которые в общем случае не связаны соотношением . Теперь рассмотрим классическую форму решения уравнения фильтра Винера, которое лежит в основе реализации фильтров, изображенных на рис. 7.9, б, 7.9, в.
Рис. 7.9. Структурные схемы фильтров, рассмотренных в примере 7.7
Фильтр Винера. Выше результаты решения стационарной задачи оценивания были получены путем введения дополнительных предположений, связанных со стационарностью задачи и позволивших упростить обобщенный алгоритм фильтрации Калмана. В частности, было установлено, что фильтр Калмана становится стационарным. Теперь сформулируем стационарную задачу оценивания в другой форме, достаточно близкой к первоначальной работе Винера. Сформулируем задачу на «частотном языке» с использованием таких понятий, как передаточные функции, спектральные плотности. На первый взгляд может показаться, что существует лишь незначительная связь между задачами Калмана и Винера. Однако ниже будет показано, что эти две задачи эквивалентны, хотя решение, полученное в форме фильтра Калмана, с точки зрения вычислений часто оказывается предпочтительнее. Задачу можно представить в виде структурной схемы, изображенной на рис 7.10. Сигнал искажается аддитивным шумом , причем и взаимно
Рис. 7.10. Представление многомерной задачи фильтрации Винера.
некоррелированные стационарные случайные процессы с
нулевым средним и со спектральными плотностями , . Наблюдение 
пропускается
через линейный фильтр с постоянными параметрами и передаточной функцией . Сигнал на выходе
фильтра обозначен как . Задача
состоит в выборе такого фильтра , на выходе которого формировалась бы
наилучшая, в смысле минимума дисперсии, оценка исходного сигнала , который получается при действии идеального оператора на сигнал . Часто
под идеальным оператором понимают тождественный (единичный) оператор, поэтому представляет собой
неискаженный сигнал . Короче говоря, необходимо выбрать передаточную функцию фильтра , которая обеспечивала
бы минимальное среднеквадратическое значение ошибки (СКО)
где 
.
Согласно теореме Парсеваля среднеквадратическая ошибка может быть выражена через матрицу спектральных плотностей ошибки :
. (7.159)
Это выражение, которое определяет среднеквадратическую ошибку как интеграл спектральной плотности ошибки, заданной в области комплексных частот, позволяет выполнить вывод уравнения фильтра в частотной области, оперируя лишь со спектральными плотностями. Этот частотный подход позволяет значительно упростить выкладки, но, очевидно, его применение ограничено лишь стационарными задачами.
Спектральная плотность ошибки, которая может быть найдена методами, рассмотренными в § 3.5, равна
Подставляя это выражение в ур-ние (7.159), получаем
Задача состоит в том, чтобы выбрать матричную передаточную функцию , минимизирующую среднеквадратическую ошибку. Для решения этой задачи представим как
. (7.162)
где - оптимальная передаточная функция, - произвольная матричная передаточная функция; - скалярная величина. Передаточная функция оптимального фильтра получается как решение уравнения
при произвольной .
Воспользуемся
этим методом для решения сформулированной задачи оценивания. Если считать, что
, то
среднеквадратическая ошибка выражается следующим образом
Здесь мы ввели два аргумента и , чтобы подчеркнуть, что СКО зависит как от , так и . Уравнение (7.163) теперь записывается в виде
(7.165)
Если воспользоваться свойством
симметрии матриц спектральной плотности и тождеством 
, то ур-ние (7.165) можно
записать в виде
Уравнение (7.166) будет удовлетворяться при произвольной , если
Это решение соответствует физически нереализуемому фильтру Винера, поскольку , вообще говоря, имеет полюса в правой полуплоскости комплексной переменной . Напоминаем, что наличие полюсов в правой полуплоскости свидетельствует не о неустойчивости системы, а скорее о физической нереализуемости системы, так как в такой системе отклик опережает воздействие.
Для того чтобы представляла собой
допустимое решение, , и должны быть физически реализуемыми или,
другими словами, должны иметь все полюса в левой полуплоскости комплексной
переменной .
Используя это ограничение на , можно выбрать , которая удовлетворяла бы
ур-нию (7.166) и была физически реализуемой.
Допустим, что матрица
спектральной плотности 
представляет собой спектр, который
допускает факторизацию в виде
где - матрица, для которой имеет все нули и полюса в левой полуплоскости комплексной переменной . Выполнение этого условия гарантирует, что и обратная матрица будут аналитическими функциями в правой полуплоскости комплексной переменной . Процедура нахождения представляет собой, вообще говоря, трудную вычислительную задачу, которую можно обычно довести до конца только численным методом при помощи достаточно сложных алгоритмов . Уравнение (7.166) теперь можно записать в виде
Представим 
в виде двух слагаемых
где объединяет все члены, имеющие полюса в левой полуплоскости, а - все члены, имеющие полюса в правой полуплоскости.
Матрицу , которую можно также
рассматривать как преобразование Лапласа части отклика фильтра, существующей
при положительном времени, назовем физически реализуемой частью фильтра и
обозначим
.
При этом полный импульсный отклик фильтра определяется как преобразование
Лапласа для величины, стоящей в правой части ур-ния (7.170).
Если подставить ур-ние (7.170) в (7.169), то необходимое условие оптимальности запишется в виде
Однако второй интеграл здесь
равен нулю, так как все полюса 
лежат в правой полуплоскости. Если
контур интегрирования заканчивается слева и ни один полюс не расположен внутри
контура, то значение интеграла равно нулю, так что ур-ние (7.171) принимает вид
. (7.172)
Поэтому оптимальный физически реализуемый фильтр имеет передаточную функцию
, (7.173)
которую можно также выразить через исходные величины
Итак, мы получили окончательное решение многомерной стационарной задачи оценивания в форме матричного фильтра Винера. Матричный фильтр Винера как решение многомерной стационарной задачи оценивания был получен Дарлингтоном , Янгом и Томасом , а также Девисом . Однако этот результат не нашел широкого применения в инженерной практике из-за достаточно трудных проблем вычислительного характера, связанных с необходимостью факторизации спектра, заданного в виде матрицы. Хотя в работах , , приведены вычислительные процедуры для факторизации спектра, которые основаны на решении матричных уравнений Риккати, широкое использование алгоритмов фильтрации Калмана вытеснило многих сторонников использования матричного фильтра Винера.
В том случае, когда сигнал и шум - некоррелированы, оптимальный фильтр имеет передаточную функцию
. (7.177)
В одномерном случае
и, наконец, когда сигнал и аддитивный шум некоррелированы, передаточная функция оптимального линейного фильтра имеет вид
Пример 7.8.
Рассмотрим
простую одномерную задачу. Спектральная плотность сигнала равна 
. Шум белый со
спектральной плотностью , причем сигнал и шум - некоррелированы.
Необходимо оценить сигнал , причем . В рассматриваемом случае
Факторизацию спектра легко выполнить, и в результате имеем:
; 
.
Отметим, что два полюса, расположенные в начале координат, были разделены так, что один из них был отнесен к правой полуплоскости, а другой - к левой полуплоскости. Используя ур-ние (7.179), получаем:
.
Рассмотрим числитель этого выражения. Разложение на элементарные дроби имеет вид
.
Функция не имеет существенного значения, так как в нее входит в первой степени, а нас интересует только та часть разложения на элементарные дроби, которая имеет полюса только в левой полуплоскости. Так как стоит в числителе , то эта функция должна иметь корни только в правой полуплоскости комплексной переменной . Следовательно, получаем
.
и передаточная функция оптимального фильтра
.
Минимальная величина дисперсии ошибки вычисляется путем подстановки выражения для в ф-лу (7.159) и последующего интегрирования по контуру. Эта часть работы значительно облегчается тем, что интегралы вида
. (7.180)
. (7.181)
табулированы для всех значений . При значения интегралов соответственно равны:
(7.182)
Однако сведение выражений к табличным интегралам часто требует выполнения громоздких алгебраических преобразований. Если сигнал и шум некоррелированны, то
Анализируя отдельно каждый член, входящий в это выражение, можно заметить, что разложение на множители, которое необходимо для сведения интегралов к табличным, легко выполняется просто разложением на множители спектральных плотностей и , которые очень часто уже заданы в факторизованном виде.
Полезным свойством этого метода является то, что он сразу дает возможность разделить среднеквадратическую ошибку на составляющую сигнала и составляющую шума. Если соответственно обозначить эти составляющие как и , то получим
. (7.186)
Пример 7.9. Воспользуемся представленным выше методом для определения минимальной величины среднеквадратической ошибки оптимального фильтра, синтезированного в примере 7.8.. Так как сигнал и шум некоррелированны, то воспользуемся упрощенными уравнениями (7.184) - (7.186). Среднеквадратическая ошибка определяется выражением
Это выражение записано в таком виде, который позволяет непосредственно использовать для вычислений табличный интеграл (7.182) с параметрами и . В результате получаем
.
Согласно выражению (7.186) шумовая составляющая среднеквадратической ошибки
.
И в этом случае воспользуемся табличным интегралом с параметрами: . В результате получаем
Наконец, находим суммарную среднеквадратическую ошибку:
Более полное исследование фильтра Винера, в том числе и ряд обобщений основной теории, интересующийся читатель может найти в литературе (см. например, , , ). Дискретный вариант фильтра Винера подробно изучен в .
Соотношение между стационарными фильтрами Калмана и Винера. В предыдущих пунктах этого параграфа были исследованы два различных метода решения стационарной задачи оценивания. Уравнение стационарного фильтра Калмана, или вырожденной формы обобщенного фильтра Калмана, было получено во временной области и выражено через переменные состояния. Уравнение фильтра Винера, напротив, было получено в частотной области и выражено через частотную характеристику. В обоих случаях вывод уравнения базировался на непосредственном использовании методов вариационного исчисления. При неглубоком анализе может показаться, что эти два подхода имеют мало общего. Однако это не так, и существует тесная связь между этими двумя подходами.
Основное отличие между задачами, сформулированными для фильтров Калмана и Винера, состоит в способе задания модели сообщения. При рассмотрении фильтра Калмана модель сообщения задается векторным дифференциальным уравнением первого порядка (7.151), а связанная с ней модель наблюдений - ур-нием (7.152). При рассмотрении фильтра Винера модель сообщения задается через спектральную плотность сигнала . Очевидно, что эти два подхода эквивалентны, так как можно найти спектральную плотность процесса, который связан с моделью сообщения, используемой при рассмотрении фильтра Калмана, как
Точно также, если задана
спектральная плотность сообщения, то можно всегда определить связанные с ней
векторные дифференциальные уравнения первого порядка, формирующие процесс с
заданной спектральной плотностью. В частности, для скалярного наблюдения можно
разложить на
два сомножителя: 
,
где имеет
все полюса и нули в левой полуплоскости, a имеет все полюса и нули в правой
полуплоскости. Если записать в виде
(7.188)
то можно построить модель сообщения с переменной фазой в виде [см. (7.151) и (7.152)]
; ; (7.189)
.
Чтобы установить эквивалентность стационарных фильтров Калмана и Винера, выберем модели сообщения и наблюдений в постановке задачи Калмана и найдем их эквивалентное спектральное представление, которое необходимо при использовании подхода Винера. Затем решим эти две задачи оценивания и сравним полученные результаты.
Предположим, что модели сообщения и наблюдений при решении задачи методом Калмана заданы уравнениями:
. (7.191)
где и - белые шумы с нулевым средним значением и ковариациями ; .
Предположим, что система, описываемая ур-нием (7.190), асимптотически устойчива и управляема. Эквивалентная спектральная плотность сигнала в постановке задачи Винера определяется как
где - резольвентная матрица;
, (7.193)
а спектральная плотность шума равна
Матрица предполагается положительно определенной, так что она допускает представление рассматривать матрицу то ур-ние (7.200) принимает вид. Используя специальное матричное тождество , получим. фильтра можно записать в виде
Вообще говоря, алгоритм Калмана в вычислительном отношении обладает преимуществом перед алгоритмом Винера главным образом благодаря тому, что он лучше приспособлен для вычислений на ЦВМ, особенно при решении многомерных стационарных задач или уравнений высокого порядка, в которых наблюдение является вектором, а также нестационарных задач. С другой стороны, имеется ряд задач, в которых факторизация спектра можжет быть выполнена в общем виде, и это позволяет глубже вникнуть в сущность полученного решения. Методом Винера могут быть также исследованы случаи небелого шума наблюдений, идеальные операции предсказания и задержки, а также учтены ограничения, связанные с насыщением и ограниченной полосой пропускания системы , . Заметим, однако, что после того, как найден оптимальный фильтр , для получения требуемой частотной характеристики необходимо еще определить его физически реализуемую часть, а это далеко не всегда простая задача.
Условие оптимальности фильтра. Фильтр Колмогорова-Винера является оптимальным фильтром формирования из входного сигнала x(t) выходного сигнала z(t) при известной форме полезного сигнала s(t), который содержится во входном сигнале в сумме с шумами. В качестве критерия его оптимизации используется среднее квадратическое отклонение сигнала y(t) на выходе фильтра от заданной формы сигнала z(t). Подставим уравнение свертки (12.2.1) в раскрытой форме весового суммирования в выражение (12.2.2") и получим отклонение e 2 выходного сигнала y(k) = h(n)③x(k-n) от заданной формы выходного сигнала z(k) по всем точкам массива данных:
В частном случае воспроизведения формы полезного сигнала в качестве функции z(k) принимается функция s(k). Минимум выражения (12.3.1) определяет значения коэффициентов h(n) оптимального фильтра. Для нахождения их значений продифференцируем выражение (12.3.1) по коэффициентам фильтра и приравняем полученные уравнения нулю. В итоге получаем:
где 
- корреляционная функция входного сигнала, 
- взаимная корреляционная функция сигналов z(k) и x(k). Отсюда:
h n R(m-n) = B(m), n = m = 0,1,2, ... , M. (12.3.2)
В краткой форме символической записи:
h(n) ③ R(m-n) = B(m). (12.3.3)
Другими словами, свертка функции отклика оптимального фильтра с функцией автокорреляции входного сигнала должна быть равна функции взаимной корреляции выходного и входного сигналов.
Система линейных уравнений фильтра. Выражение (12.3.2) может быть записано в виде системы линейных уравнений - однострочных равенств левой и правой части для фиксированных значений координаты m коэффициентов фильтра. При расчете симметричных фильтров, не сдвигающих фазы выходного сигнала, фильтр может вычисляться только одной половиной своих значений:
m=0: h o R(0)+ h 1 R(1)+ h 2 R(2)+ h 3 R(3)+ ...+ h M R(M) = B(0), (12.3.3")
m=1: h o R(1)+ h 1 R(0)+ h 2 R(1)+ h 3 R(2)+ ...+ h M R(M-1) = B(1),
m=2: h o R(2)+ h 1 R(1)+ h 2 R(0)+ h 3 R(1)+ ...+ h M R(M-2) = B(2),
..............................................................................................................
m=M: h o R(M)+ h 1 R(M-1)+ h 2 R(M-2)+ .... + h M R(0) = B(M).
Решение данной системы уравнений относительно значений h i дает искомый оператор фильтра.
При расчете каузальных (односторонних) фильтров выходной сигнал z(t) следует задавать со сдвигом вправо по оси координат таким образом, чтобы значимая часть функции взаимной корреляции B(m) полностью располагалась в правой части системы уравнений (12.3.3").
Отметим, что R(m) = R s (m)+R q (m), где R s - функция автокорреляции сигнала, R q - функция автокорреляции шума, а B(m) = B zs (m)+B zq (m), где B zs - функция взаимной корреляции сигналов z(k) и s(k), B zq - функция взаимной корреляции сигнала z(k) и помех q(k). Подставляя данные выражения в (12.3.3), получаем:
h(n) ③ = B zs (m)+B zq (m). (12.3.4)
Частотная характеристика фильтра находится преобразованием Фурье левой и правой части уравнения (12.3.4):
H(w) = W zs (w)+W zq (w),
H(w) = / , (12.3.5)
где W s (w) ó R s (m) и W q (w) ó R q (m) - энергетические спектры (плотности мощности) сигнала и помех, W zs (w) ó B zs (m) - взаимный энергетический спектр входного и выходного сигналов, W zq (w) ó B zq (m) - взаимный энергетический спектр выходного сигнала и помех.
В геофизической практике обычно имеет место статистическая независимость полезного сигнала, а, следовательно, и сигнала z(k), от шумов, при этом B zq = 0 и фильтр называют оптимальным по сглаживанию шумов при заданной форме выходного сигнала:
H(w) = W zs (w) / , (12.3.6)
Фильтр (12.3.6) оптимален в том смысле, что максимизирует отношение мощности сигнала к мощности шума по всему интервалу сигнала, но не в каждой индивидуальной точке.
Выражения (12.3.5-12.3.6), как правило, всегда имеют решение. Однако это не означает возможность формирования фильтром любой заданной формы выходного сигнала. Из чисто практических соображений можно сразу предполагать, что если спектр заданного сигнала z(t) больше значимой части спектра полезного сигнала s(t), то оператор фильтра попытается сформировать требуемые высокие частоты заданного сигнала из незначимых частот спектра полезного сигнала, что может потребовать огромных коэффициентов усиления на этих частотах, под действие которых попадут и частотные составляющие шумов. Результат такой операции непредсказуем. Эти практические соображения можно распространить и на все частотные соотношения входного и выходного сигналов линейных фильтров: значимые гармоники спектров выходных сигналов должны формироваться из значимых гармоник спектров входных сигналов.
Если заданная форма сигнала z(k) совпадает с формой полезного сигнала s(k), то B(m) = B ss = R s и фильтр называют фильтром воспроизведения полезного сигнала :
H(w) = W s (w) / , (12.3.7)
Выражения (12.3.6-12.3.7) достаточно наглядно демонстрируют физический смысл формирования передаточной функции фильтра. При воспроизведении сигнала частотная функция взаимной корреляции входного сигнала с выходным W zs (плотность взаимной мощности) повторяет частотную функцию автокорреляции W s (плотность мощности сигнала). Плотность мощности статистических шумов W q распределена по частотному диапазону равномерно, в отличие от плотности мощности сигнала W s , которая, в зависимости от формы сигнала, может занимать любые частотные интервалы спектрального диапазона. Частотная передаточная функция фильтра воспроизведения сигнала формируется отношением W s (w)/. На частотах, где сосредоточена основная энергия сигнала, имеет место W s (w)>>W q (w) и H(w) Þ 1 (как минимум, больше 0.5). Там, где значение W s (w) становится меньше W q , коэффициент передачи фильтра становится меньше 0.5, и в пределе H(w)=0 на всех частотах, где полностью отсутствуют частотные составляющие сигнала.
Аналогичный процесс имеет место и при формировании произвольного сигнала z(t) на выходе фильтра, только в этом случае из частот входного сигнала устанавливаются на выделение и усиление частоты взаимной мощности входного и выходного сигнала, необходимые для формирования сигнала z(t), причем коэффициент на этих частотах может быть много больше 1, а подавляться могут не только шумы, но и частоты основного сигнала, если их нет в сигнале z(t).
Таким образом, оптимальные фильтры учитывают особенности спектрального состава сигналов и способны формировать передаточные функции любой сложности на выделение полезных частот сигналов из любых диапазонов спектра с максимальных подавлением шумов на всех частотах спектрального диапазона, не содержащих полезных сигналов, при этом границы усиления-подавления устанавливаются автоматически по заданному уровню шумов.
Задание мощности шумов. Следует внимательно относиться к заданию функции шумов Wq(w). При полной неопределенности шума необходимо, как минимум, выполнять оценку его дисперсии s 2 и распространять на весь частотный диапазон с соответствующей нормировкой на его величину (2Wq(w) dw = s 2), т.е. считать его белым шумом. При известной функции полезного сигнала в зарегистрированной реализации значение дисперсии шумов в первом приближении может быть оценено по разности дисперсий реализации и функции полезного сигнала. Можно выполнить и выделение шумов из входного сигнала в отдельный шумовой массив, например, вейвлетным преобразованием. Однако использовать выделенный шум непосредственно для вычисления функции Wq(w) допустимо только по достаточно представительной реализации при условии стационарности и эргодичности шума. В противном случае функция Wq(w) будет отображать только распределение шумов в зарегистрированной реализации сигнала, а соответственно фильтр будет оптимален только к этой реализации, что не гарантирует его оптимальности к любой другой реализации. Но для обработки единичной зарегистрированной реализации сигнала такой метод не только вполне допустим, но и может существенно повысить точность формирования выходного сигнала.
Эффективность фильтра.
Из выражений (12.3.5-12.3.7) следует, что с позиции минимального искажения полезного сигнала при максимальном подавлении шумов фильтр Колмогорова-Винера эффективен в тем большей степени, чем больше отношение сигнал/шум на входе фильтра. В пределе, при W q (w)<
Пример. Расчет оптимального фильтра воспроизведения сигнала. Выполняется в среде Mathcad.
Форма входного сигнала считается известной и близка к гауссовой. На входной сигнал наложен статистический шум с равномерным распределением мощности по всему частотному диапазону (белый шум), некоррелированный с сигналом, и функцию Wzq принимаем равной нулю. Для наглядного просмотра связи параметров фильтра с параметрами сигнала задаем модели сигналов в двух вариантах:
K:= 1000 k:= 0 .. K A:= 50
s1 k:= A·exp[-0.0005·(k-500) 2 ] s2 k:= A·exp[-0.00003·(k-500) 2 ] Ü информационные сигналы
Q:= 30 q k:= rnd(Q) – Q/2 x1 k:= s1 k + q k x2 k:= s2 k + q k Ü входные сигналы
Рис. 12.3.1. Модельные сигналы.
В качестве выходных сигналов задаем те же самые функции s1 и s2. Быстрым преобразованием Фурье вычисляем спектры сигналов и формируем спектры плотности мощности:
S1:= CFFT(s1) S2:= CFFT(s2) Q:= CFFT(q) Ü спектры сигналов
Ü спектры мощности
Ds1:= var(s1) Ds2:= var(s2) Dx1:= var(x1) Dx2:= var(x2) Dq:= var(q) Ü дисперсии
Ds1 = 124.308 Ds2 = 310.264 Dx1 = 202.865 Dx2 = 386.78 Dq = 79.038 Ü информация
mean(Wq) = 0.079 Wq1:= (Dx1 – Ds1)/(K+1) Wq1 = 0.078 Ü информация
Wq2:= (Dx2 – Ds2)/(K+1) Wq2 = 0.076 Ü информация
Wq k:= Wq1 Ü замена на постоянное распределение
Для воспроизведения сигналов вычисления функций Wzs не требуется, т.к. Wzs = Ws. Вычисление Wq также имеет только информативный характер.
Передаточные функции оптимальных фильтров (приведены на рис. 12.3.2):
Рис. 12.3.2. Передаточные функции оптимальных фильтров
в сопоставлении с нормированными модулями спектров сигналов
Как следует из рисунка 12.3.2, для плавных монотонных функций, основная энергия которых сосредоточена в низкочастотной части спектра, передаточные функции оптимальных фильтров, по существу, представляют собой низкочастотные сглаживающие фильтры с автоматической подстройкой граничной частоты пропускания под основные частоты входного сигнала. Операторы фильтров можно получить обратным преобразованием Фурье:
h1:= ICFFT(H1)/(K+1) h2:= ICFFT(H2)/(K+!) Ü обратное преобразование Фурье
Рис. 12.3.3. Импульсные отклики фильтров.
Оператор фильтра, в принципе, бесконечен. В данном случае, при использовании БПФ максимальное число отсчетов равно К/2 = 500. Усечение размеров оператора может выполняться по типовым методам с применением весовых функций (в расчете операторы нормируются к 1, весовые функции не применяются).
N1:= 160 n1:= 0 .. N1 N2 ;= 500 n2:= 0 .. N2 Ü размеры и нумерация операторов
hm1:= h1 0 + 2·h1 n 1 hm1=0.988 h1:= h1/hm1 Ü нормировка
hm2:= h2 0 + 2·h2 n 2 hm2=1.001 h2:= h2/hm2 Ü нормировка
Ü свертка
Рис. 12.3.4. Проверка действия оптимальных фильтров.
Коэффициент усиления дисперсии помех Þ Kd:= (h 0) 2 + 2·h n Kd1=0.021 Kd2= 0.0066
Среднеквадратическое отклонение воспроизведения сигнала:
e1= 1.238 e2 = 0.701
Проверка действия оператора фильтра приведена на рис. 12.3.4.
Особую эффективность оптимальный фильтр имеет при очистке от шумов сигналов, имеющих достаточно сложный спектральный состав. Оптимальный фильтр учитывает конфигурацию спектра сигнала и обеспечивает максимальное подавление шумов, в том числе внутри интервала основного частотного диапазона сигнала. Это наглядно видно на рис. 12.3.5 для сигнала, близкого к прямоугольному, спектр которого имеет кроме основной низкочастотной части затухающие боковые осцилляции. Расчет фильтра выполнялся по методике, приведенной в примере 1.
Рис. 12.3.5. Оптимальная фильтрация сигнала со сложным спектральным составом.
Рис. 12.3.6. Оптимальная фильтрация радиоимпульса.
На рис. 12.3.6 приведен пример фильтрации оптимальным фильтром радиоимпульса. Основной пик спектра радиоимпульса находится в области несущей частоты, а боковые полосы определяются формой модулирующего сигнала, в данном случае – прямоугольного импульса. На графике модулей сигнала и передаточной функции фильтра можно видеть, что оптимальный фильтр превратился в полосовой фильтр, при этом учитывается форма боковых полос спектра сигнала.
Фильтры прогнозирования и запаздывания. Если в правой части уравнения (12.3.3) желаемым сигналом задать входной сигнал со сдвигом на величину kDt, то при этом B(m) = R(m+k) и уравнение принимает вид:
h(n) ③ R(m-n) = R(m+k). (12.3.8)
При k > 0 фильтр называется фильтром прогнозирования и вычисляет будущие значения сигнала по его предшествующим значениям. При k < 0 фильтр является фильтром запаздывания. Реализация фильтра заключается в решении соответствующих систем линейных уравнений для каждого заданного значения k. Фильтр может использоваться для интерполяции геофизических полей, в том числе в наперед заданные точки, а также для восстановления утраченных данных.