Адаптивные фильтры. Применение адаптивных фильтров в идентификации систем. Адаптивная фильтрация

В этой книге рассматриваются теория, расчет и применение адаптивных фильтров. Первый адаптивный, или самообучающийся, фильтр часто приписывают Лаки из-за разработанного им в 1966 г. обнуляющего корректирующего фильтра, компенсирующего искажения в системах передачи данных. Однако более ранняя работа по адаптивному распознаванию формы сигнала была выполнена в 1960 г. Яковацом и др. . В 1961 г. Глезер в США провел теоретическое исследование по адаптивным фильтрам, а Габор и др., в том же году, в Великобритании, воспользовались аналоговым лентопротяжным механизмом для подстройки весов нелинейного «обучающегося» фильтра. Мы можем считать, что название «обучающийся» относится к адаптивному процессору.

Яндекс.ДиректВсе объявленияОптимизация оплаты. Оклады Премии Семинар 26-27.02 в Москве, 11-12.03 в Киеве. Первым лицам скидки до 50%. lityagin.ru Решение задач по механике он-лайн Более 1000 преподавателей он-лайн! Решение задач по механике! liveexpert.ru

Большинство ранних работ по адаптивным фильтрам выполнено в ходе независимых исследований различными научно-исследовательскими организациями. Заслуживающие упоминания исследования проводились в Высшей технической школе г. Карлсруэ в ФРГ и в Станфордском университете, где в 1959 г. было начато создание адаптивных систем распознавания образов. В ходе совместной работы этих организаций в 1964 г. была произведена сравнительная оценка каждого метода , что впоследствии привело к разработке наиболее широко используемого алгоритма для подстройки весовых коэффициентов процессора. Дальнейшая работа в данном направлении одновременно проводилась в Институте автоматики и телемеханики в Москве. В середине 60-х годов прекрасный сводный обзор по адаптивным фильтрам и предварительные рекомендации по их применению для адаптивного или автоматического выравнивания был представлен в работе . Позднее были подготовлены несложные обзорные статьи по гашению отраженного сигнала в телефонии и адаптивному выравниванию .

Для получения оптимального решения существует много методов подстройки значений весовых коэффициентов фильтра. Применялись методы случайных возмущений , которые изменяли весовые коэффициенты фильтра; далее анализировался выходной сигнал для того, чтобы установить, приближает ли его случайное возмущение к искомому решению или отдаляет от него. В гл. 3 подробно рассматривается разработка адаптивного алгоритма метода наименьших квадратов (МНК), который использовался в работе Станфордского университета, по распознаванию образов и впервые был официально описан в 1967 г. Уидроу и др. для адаптивных антенных решеток, а в 1971 г. для адаптивных фильтров . В настоящее время этот алгоритм широко применяется для расчета весовых коэффициентов адаптивных фильтров, поскольку в нем используются градиентные методы, которые значительно эффективнее, чем другие, обеспечивают сходимость к оптимальному решению. Можно показать, что градиентный метод наименьших квадратов очень схож с методом максимизации отношения сигнал – шум, который параллельно разрабатывал Эпплбаум c целью применения в тех случаях, когда необходимо получить оптимальные весовые коэффициенты адаптивных антенных решеток. Было также показано, что обнуляющий корректирующий фильтр Лаки является упрощением более общего градиентного метода наименьших квадратов.

ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ

Digital signals processing

Тема 11. АДАПТИВНАЯ ФИЛЬТРАЦИЯ ЦИФРОВЫХ ДАННЫХ

Пусть они постараются подчинить себе обстоятельства, а не подчиняются им сами.

Гораций. Послания. Римский поэт, I в.д.н.э.

Если в этой теории Вы не увидите смысла, тем лучше. Можно пропустить объяснения и сразу приступить к ее практическому использованию.

Валентин Ровинский. Теория карточных игр.

Киевский геофизик Уральской школы, ХХ в.
Содержание

Введение.

1. Общие сведения об адаптивной . Основные области применения. Адаптивный шумоподавитель. Адаптивный фильтр Винера. Адаптивный алгоритм наименьших квадратов Уидроу-Хопфа. Рекурсивные схемы наименьших квадратов.

2. Основы статистической группировки информации. Предпосылки метода. Задача статистической группировки. Использование априорных данных. Эффективность метода.

Статистическая регуляризация данных. Проверка теоретических положений метода. Оценка сохранения разрешающей способности. Статистическая оценка регуляризации данных. Результаты моделирования. Частотное представление. Пример практического использования.

4. Статистическая группировка полезной информации. Сущность аппаратной реализации. Особенности аппаратной реализации. Реализация систем группировки информации. Пример исполнения системы группировки информации.

Введение

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

Адаптивные устройства обработки данных отличаются наличием определенной связи параметров передаточной функции с параметрами входных, выходных, ожидаемых, прогнозируемых и прочих дополнительных сигналов или с параметрами их статистических соотношений, что позволяет самонастраиваться на оптимальную обработку сигналов. В простейшем случае, адаптивное устройство содержит программируемый фильтр обработки данных и блок (алгоритм) адаптации, который на основании определенной программы анализа входных, выходных и прочих дополнительных данных вырабатывает сигнал управления параметрами программируемого фильтра. Импульсная характеристика адаптивных систем также может иметь как конечный, так и бесконечный характер.

Как правило, адаптивные устройства выполняются узкоцелевого функционального назначения под определенные типы сигналов. Внутренняя структура адаптивных систем и алгоритм адаптации практически полностью регламентируются функциональным назначением и определенным минимальным объемом исходной априорной информации о характере входных данных и их статистических и информационных параметрах. Это порождает многообразие подходов при разработке систем, существенно затрудняет их классификацию и разработку общих теоретических положений /л38/. Но можно отметить, что наибольшее применение при разработке систем для адаптивной обработки сигналов находят два подхода: на основе схемы наименьших квадратов (СНК) и рекурсивной схемы наименьших квадратов (РСНК).

^ 11.1. ОБЩИЕ СВЕДЕНИЯ ОБ АДАПТИВНОЙ ЦИФРОВОЙ ФИЛЬТРАЦИИ .

Основные области применения адаптивной фильтрации – очистка данных от нестабильных мешающих сигналов и шумов, перекрывающихся по спектру со спектром полезных сигналов, или когда полоса мешающих частот неизвестна, переменна и не может быть задана априорно для расчета параметрических фильтров. Так, например, в цифровой связи сильная активная помеха может интерферировать с полезным сигналом, а при передаче цифровой информации по каналам с плохими частотными характеристиками может наблюдаться межсимвольная интерференция цифровых кодов. Эффективное решение этих проблем возможно только адаптивными фильтрами.

Частотная характеристика адаптивных фильтров автоматически регулируется или модифицируется в соответствии с определенным критерием, позволяющем фильтру адаптироваться к изменениям характеристик входного сигнала. Они достаточно широко используются в радио- и гидролокации, в системах навигации, в выделении биомедицинских сигналов, и многих других отраслях техники. В качестве примера рассмотрим наиболее распространенные схемы адаптивной фильтрации сигналов.

Адаптивный шумоподавитель . Блок-схема фильтра приведена на рис. 11.1.1.

Рис. 11.1.1.
Фильтр состоит из блока цифрового фильтра с регулируемыми коэффициентами и адаптивного алгоритма для настройки и изменения коэффициентов фильтра. На фильтр одновременно подаются входные сигналы y(k) и x(k). Сигнал y(k) содержит полезный сигнал s(k) и некоррелированный с ним загрязняющий сигнал g(k). Сигнал x(k) какого-либо источника шума коррелирован с g(k) и используется для формирования оценки сигнала ğ(k). Полезный сигнал оценивается по разности:

š(k) = y(k) – ğ(k) = s(k) + g(k) – ğ(k). (11.1.1)

Возводим уравнение в квадрат и получаем:

š 2 (k) = s 2 (k) + (g(k) – ğ(k)) 2 + 2.s(k) (g(k) – ğ(k)). (11.1.2)

Вычислим математическое ожидание левой и правой части этого уравнения:

M[š 2 (k)] = M + M[(g(k) – ğ(k)) 2 ] + 2M. (11.1.3)

Последнее слагаемое в выражении равно нулю, поскольку сигнал s(k) не коррелирует с сигналами g(k) и ğ(k).

M[š 2 (k)] = M + M[(g(k) – ğ(k)) 2 ]. (11.1.4)

В этом выражении M = W(s(k)) – мощность сигнала s(k), M[š 2 (k)] = W(š(k)) – оценка мощности сигнала s(k) и общая выходная мощность, M[(g(k) – ğ(k)) 2 ] = W( g) - остаточная мощность шума, который может содержаться в выходном сигнале. При настройке адаптивного фильтра к оптимальному положению минимизируется мощность остаточного шума, а, следовательно, и мощность выходного сигнала:

Min W(š(k)) = W(s(k)) + min W( g). (11.1.5)

На мощность полезного сигнала настройка не влияет, поскольку сигнал не коррелирован с шумом. Эффект минимизации общей выходной мощности будет выражаться в максимизации выходного отношения сигнал/шум. Если настройка фильтра обеспечивает равенство ğ(k) = g(k), то при этом š(k) = s(k). Если сигнал не содержит шума, адаптивный алгоритм должен устанавливать нулевые значения всем коэффициентам цифрового фильтра.

Рис. 11.1.2.
Адаптивный фильтр Винера . Входной сигнал y(k) фильтра, приведенного на рис. 11.1.2, включает компоненту, коррелированную со вторым сигналом x(k), и полезную компоненту, некоррелированную с x(k). Фильтр формирует из x(t) сигнал ğ(k) - оптимальную оценку той части у(k), которая коррелированна с x(k), и вычитает ее из сигнала y(k). Выходной сигнал:

E(k) = y(k) - ğ(k) = y(k) - H T X k = y(k) -h(n) x(k-n),

Где H T и X k – векторы весовых коэффициентов фильтра и его входного сигнала.

Аналогично предыдущему методу, возводим в квадрат левую и правую части уравнения, находим математические ожидания обеих частей и получаем уравнение оптимизации  выходного сигнала:

   2P T H + H T RH , (11.1.6)

Где  2 = M – дисперсия y(k), P = M – вектор взаимной корреляции, R = M[X k X k T ] – автокорреляционная матрица.

Рис. 11.1.3.
В стационарной среде график зависимости  от коэффициентов H представляет собой чашеобразную поверхность адаптации (рис. 11.1.3). Градиент поверхности:

d / dH = -2P + 2RH .

Каждому набору коэффициентов h(n) на этой поверхности соответствует определенная точка. В точке минимума градиент равен нулю и вектор весовых коэффициентов фильтра является оптимальным:

H opt = R -1 P . (11.1.7)

Эта формула называется уравнением Винера-Хопфа. Задачей алгоритма автоматической настройки является подбор таких весовых коэффициентов фильтра, которые обеспечивают работу в оптимальной точке поверхности адаптации.

Однако практическое применение фильтра затрудняется использованием корреляционных матриц R и P, априори неизвестных, и которые могут изменяться со временем для нестационарных сигналов.

Адаптивный алгоритм наименьших квадратов Уидроу-Хопфа . По существу, это модификация фильтра Винера, в которой вместо вычисления коэффициентов (11.1.7) за один шаг используется алгоритм последовательного спуска в оптимальную точку при обработке каждой выборки:

H k +1 = H k - e k X k , (11.1.8)

E k = y k - H T X k . (11.1.9)

Условие сходимости к оптимуму:

0 <  > 1/ max , (11.1.10)

Где  - параметр скорости спуска,  m ax – максимальное собственное значение ковариационной матрицы данных. Блок-схема алгоритма приведена на рис. 11.1.4.

Рис. 11.1.4. Алгоритм адаптации методом наименьших квадратов.

На практике точка максимальной оптимальности флюктуирует около теоретически возможной. Если входной сигнал нестационарный, то изменение статистик сигнала должно происходить достаточно медленно, чтобы коэффициенты фильтра успевали следить за этими изменениями.

Рекурсивные схемы наименьших квадратов отличаются тем, что вычисление каждой последующей выборки коэффициентов h(n) производится не только по коэффициентам только одной предыдущей выборки, но и с определенной длиной постепенно затухающей памяти по предшествующим выборкам, что позволяет снижать флюктуации оценок при обработке стационарных сигналов.

^ 11.2. Основы статистической группировки информации.

При построении систем адаптивной фильтрации данных большое значение имеют статистические характеристики обрабатываемых сигналов и шумов, их стационарность, и наличие какой-либо дополнительной информации, коррелированной с основной. Возможность использования дополнительной информации при построении адаптивных систем рассмотрим на конкретном примере – системе адаптивной фильтрации данных непрерывных ядерногеофизических измерений.

Предпосылки метода. Физической величиной, регистрируемой в процессе ядерно-физических измерений в геофизике, обычно является частота импульсных сигналов на выходе детекторов ионизирующего излучения в интегральном или дифференциальном режиме амплитудной селекции. Значения измеряемой величины, как статистически распределенной по своей природе, могут быть определены только путем усреднения числа актов регистрации ионизирующих частиц по интервалам времени. Зарегистрированное количество импульсов определяет статистическую погрешность единичного измерения, а временной интервал усреднения, обеспечивающий нормативную погрешность – их производительность. Для методов с непрерывной регистрацией информации во времени (или в пространстве) временное окно измерений определяет также временную (или пространственную, с учетом скорости перемещения детектора) разрешающую способность интерпретации результатов измерений, при этом эффективность регистрации информации обычно ограничена условиями измерений и/или техническими средствами их исполнения. Типичный пример - каротаж скважин, где возможности увеличения интенсивности потоков информации ограничены параметрами эффективности регистрации и чувствительности детекторов излучения, которые зависят от их типа и размеров. Размеры детекторов, естественно, существенно зависят от размеров скважинных приборов, которые, в свою очередь, ограничены диаметрами скважин.

Ниже рассматривается возможность повышения точности и производительности непрерывных ядерно-физических измерений, для наглядности, применительно к условиям измерений в варианте скважинного гамма-опробования, хотя в такой же мере она может быть использована в авто- и аэрогаммасъемке, при радиометрическом обогащении руд, в рентгенорадиометрии и других методах ядерной геофизики. Предполагается, что регистрация данных производится в цифровой форме с накоплением отсчета по постоянным интервалам дискретизации данных (по времени и по пространству, при условии постоянной скорости перемещения детектора).

В общем случае полезная (целевая) информация может присутствовать в нескольких энергетических интервалах спектра излучения. Рабочими интервалами измерений обычно счи-таются участки спектра, где полезная информация присутствует в "чистом" виде либо в смеси с помехами (фоном), значение которых может быть учтено при обработке результатов изме-рений. Так, например, при гамма-опробовании пород на содержание естественных радионук-лидов (ЕРН) регистрируется излучение с энергией более 250-300 кэВ, представленное в основном первичными и однократно рассеянными квантами, плотность потока которых про-порциональна массовой доле ЕРН в породах. Плотность потока излучения в низкоэнергети-ческом интервале спектра (20-250 кэВ, в основном многократно рассеянное излучение) также зависит от массовой доли ЕРН, но эта зависимость является параметрически связанной с эф-фективным атомным номером излучающе-поглощающей среды в области детектора, вариации которого по стволу скважины могут приводить к большой погрешности интерпретации ре-зультатов измерений. Между тем плотность потока информации (относительно массовой доли ЕРН) в интервале 20-250 кэВ много выше, чем в интервале более 250 кэВ, осо-бенно при регистрации излучения сцинтилляционными детекторами малых объемов, которые имеют повышенную чувствительность именно к низкоэнергетической части спектра излуче-ния.

Задача статистической группировки информации в потоках сигналов в общей и наибо-лее простой форме может быть сформулирована следующим образом. Полезная информация присутствует в двух статистически независимых потоках сигналов (в двух неперекрывающихся интервалах спектра излучения). В первом потоке сигналов, условно- основном, полезная информация присутствует в "чистом" виде: плотность потока сигналов пропорциональна опреде-ляемой физической величине. Во втором потоке, условно-дополнительном, на полезную ин-формацию наложено влияние дестабилизирующих факторов, значение которых неизвестно. При отсутствии дестабилизирующих факторов коэффициент корреляции средних значений плотностей потоков в этих двух потоках сигналов постоянен и близок к 1. Для снижения ста-тистической погрешности измерений требуется осуществить извлечение полезной информа-ции из дополнительного потока сигналов и ее суммирование с основным потоком.

Обозначим потоки, а равно и частоты основного и дополнительного потоков сигналов индексами n и m (импульсов в секунду), связь потоков по частотам индексом х = m/n. Определению подлежит частота потока n. Значение х может изменяться за счет влияния дестабилизирующих факторов на поток m и в общем случае представляет собой случайную величину, распределенную по определенному закону с плотностью вероятностей Р(х), мате-матическим ожиданием , и дисперсией D x .

На основе теоремы Байеса, плотность вероятностей распределения частоты n по измеренному за единичный интервал t числу отсчетов сигнала N определяется выражением:

P N (n) = P(n) P n (N) P(N), (11.2.1)

P n (N) = (nТ) N e -n  N! , (11.2.2)

P(N) =P n (N) P(n) dn, (11.2.3)

Где: P(n) - априорная плотность вероятностей частоты n, P n (N) - апостериорное распределе-ние вероятностей числовых отсчетов N (закон Пуассона). Принимая в дальнейшем в качестве искомой величины значения отсчетов z=n по интервалам  (экспозиция цифровых отсчетов или скользящее временное окно аналоговых данных) и подставляя (11.2.2, 11.2.3) в (11.2.1), получаем:

P N (z) = P(z) z N e -z P(z) z N e -z dz. (11.2.4)

При неизвестном распределении значений z априорная плотность распределения P(z) принимается равномерной от 0 до , при этом из выражения (11.2.4) следуют общеизвестные выражения:

Z = D z = N+1  N, (11.2.5)

 z 2 = D z z 2 = 1 (N+1)  1N. (11.2.6)

Значениями единиц в выражениях пренебрегаем, что не только корректно в условиях "хорошей" статистики, но и необходимо в режиме последовательных непрерывных измерений для исключения смещения средних значений.

Как следует из теории гамма-каротажа (ГК) и достаточно хорошо подтверждено практикой гамма-опробования, пространственная разрешающая способность гамма-каротажных измерений при интерпретации результатов ГК на содержание естественных радиоактивных элементов в породах по стволу скважин в среднем составляет 10 см, а в скважинах малого диаметра может даже повышаться до 5-7 см. Однако реализация такой разрешающей способности возможна только в условиях достаточно "хорошей" статистики. Коэффициент усиления дисперсии помех цифровых фильтров деконволюции, которые используются при интерпретации ГК, в среднем порядка 12 и изменяется от 4 до 25 в зависимости от плотности пород, диаметра скважин, диаметра скважинных приборов и пр. Отсюда следует, что для достижения разрешающей способности в 10 см при нормативной погрешности дифференциальной интерпретации не более 10-20 % статистическая погрешность измерений не должна превышать 3-7 %. А это, в свою очередь, определяет объем отсчета за единичную экспозицию не менее 200-1000 импульсов. При гамма-каротаже последнее возможно только для пород с относительно высоким содержанием ЕРН (более 0.001 % эквивалентного урана), при использовании детекторов больших размеров (с эффективностью регистрации более 10 имп/сек на 1 мкР/час) и при низкой скорости каротажа (не более 100-300 м/час). В той или иной мере эта проблема характерна для всех методов ядерной геофизики, и особенно остро стоить в спектрометрических модификациях измерений.

Вместе с тем следует отметить, что процесс непрерывных измерений имеет определенную физическую базу как для применения методов регуляризации результатов интерпретации данных, так и для регуляризации непосредственно самих статистических данных (массивов отсчетов N) при их обработке.

Простейшим способом подготовки цифровых данных для интерпретации является их низкочастотная фильтрация методом наименьших квадратов (МНК) или весовыми функциями (Лапласа-Гаусса, Кайзера-Бесселя и др.). Однако любые методы низкочастотной фильтрации данных снижают пространственную разрешающую способность интерпретации, так как кроме снижения статистических флюктуаций приводят к определенной деформации частотных составляющих полезной части сигнала, спектр которого по условиям деконволюции должен иметь вещественные значения вплоть до частоты Найквиста. В определенной мере ликвидировать этот негативный фактор позволяет метод адаптивной регуляризации данных (АРД).

Выражения (11.2.5-6) получены в предположении полной неизвестности априорного распределения P(z)для отсчетов в каждой текущей экспозиции . Между тем, при обработке данных непрерывных измерений, и тем более каротажных данных, которые обычно являются многопараметровыми, для каждого текущего отсчета при обработке данных может проводиться определенная оценка распределения P(z). Как минимум, можно выделить два способа оценки распределения P(z).

Способ 1. По массивам данных параллельных измерений каких-либо других информационных параметров, значения которых достаточно четко коррелированны с обрабатываемым массивом данных либо в целом по пространству измерений, либо в определенном скользящем интервале сравнения данных. К таким массивам относятся, например, предварительные каротажные измерения в процессе бурения скважин, измерения другим прибором, с другой скоростью каротажа, в другом спектральном интервале излучения, и даже другим методом каротажа. При гамма-опробовании оценка распределения P(z) может производиться по параллельным измерениям интенсивности потока m в низкочастотном интервале спектра горных пород.

Способ 2. При единичной диаграмме ГК оценка распределения P(z) в каждой текущей точке обработки данных может выполняться по ближайшим окрестностям данной точки, захватывающим более широкий пространственный интервал по сравнению с интервалом отсчетов.

Использование априорных данных. Допустим, что кроме основного массива данных N, подлежащего обработке (подготовке к интерпретации), мы располагаем дополнительным массивом данных M, значения которого в определенной степени коррелированы с массивом N. При отсутствии дополнительных массивов способ 2 позволяет получить массив М обработкой массива N цифровым фильтром МНК (или любым другим весовым фильтром) со скользящим временным окном T  3 (M(k) = m(k)сглаженного сигнала m(k) = n(k) ③ h, где h – оператор симметричного цифрового фильтра). Отметим также, что 2-ой способ всегда может использоваться для регуляризации данных независимо от наличия данных для 1-го метода.

Массив М позволяют дать оценку статистических характеристик распределения P(z). Так, если для тех же интервалов времени  в массиве М имеются отсчеты М = m k  (или приведенные к ним отсчеты какого-либо другого параметра), то можно записать:

P M (z) =
, (11.2.7)

Где Р(х) – априорная плотность распределения значений x k = m k /n k , которые в общем случае также могут быть случайными. При равномерном распределении Р(х) от 0 до  для отсчета М равновероятно любое значение z, т.е. эффект от измерений в потоке m отсутствует. Однако по исходным условиям задачи в потоке m обязательно присутствие полезной информации, а, следовательно, и суще-ствование, как минимум, определенных границ распределения Р(х) от х min > 0 до x max << , и среднего значения по пространству измерений. При этом из выражения (11.2.7) следует, что наиболее вероятное значение z a , "априорное" для отсчетов z=n в потоке n по измерениям в потоке m (отсчетам М), должно быть равно:

Z a = (M+1)  М. (11.2.8)

При статистической независимости величин х и М относительная средняя квадратическая погрешность определения значений z a по отсчетам в массиве М:

 za 2 =  M 2 +  x 2 . (11.2.9)

Отсюда дисперсия распределения значений z a:

D za = (D M +M 2  x 2) 2 = D(M)  2 , (11.2.10)

D(M) = D M +M 2  x 2 = D M +D xm , (11.2.11)

D M = М+1  М, D xm = M 2  x 2 ,

Где значение дисперсии D M определяется статистикой отсчетов в массиве М при х = const, значение D xm представляет собой дисперсию значений М за счет флюктуаций величины х, а сумма D(M) определяет полную дисперсию отсчетов М.

Влияние Р(х) на форму распределения Р М (z) сказывается в его "растягивании" по координате z относительно модального значения, при этом решение интеграла (11.2.7) в первом приближении может быть представлено в следующем виде:

P M (z)  b
e -bz . (11.2.12)

Для данного распределения:

= z a = ab, (11.2.13)

D za = ab 2 , (11.2.14)

С учетом выражений (11.2.8) и (11.2.10):

A = MD M (D za 2) = MD M D(M), (11.2.15)

B = D M (D za ) = D М D(M). (11.2.16)

Значение "а" в выражении (11.2.15) принимается целочисленным. Выражение (11.2.12) может быть принято для распределения (11.2.4) в качестве априорного распределения вероятностей Р(z), при этом:

P N (z) = (b+1)
e -z(b+1) . (11.2.17)

Отсюда, математическое ожидание и дисперсия z:

Z = (N+a)(b+1), (11.2.18)

D z = (N+a)(b+1) 2 . (11.2.19)

C использованием выражений (11.2.15-16):

Z = N+(1-)M, (11.2.20)

Где  и (1-) – весовые коэффициенты доверия отсчетам N и M:

 = D(M)(D N 2 +D(M)). (11.2.21)

Дисперсия и относительная средняя квадратическая погрешность отсчетов z:

D z = D(M)
, (11.2.22)

 z 2 =1(N+MD M D(M)). (11.2.23)

Эффективность метода. Сравнение выражений (11.2.20-23) и (11.2.5-6) позволяет дать оценку эффекта использования дополнительной информации из статистически независимого от N потока М (произвольная дополнительная информация).

1. При  const имеет место  х 2  0, D xm  0 и дисперсия отсчетов в массиве М определяется только статистикой потока:

D(M)  D M = M, z = (N+M) (+1),

 z 2  1(N+M) <  N 2 = 1N, (11.2.24)

 =  N 2  z 2 = N  1+MN,

Что соответствует определению z по двум независимым измерениям и эффект использования дополнительной информации максимален. Так, при M  N,   2 и погрешность измерений уменьшается в
1.4 раза.

2. В общем случае D xm  0, при этом D(M) > D М и положительный эффект снижается. В пределе:  x  , D xm  , D(M)  ,   1, z  N,  z   N и положительный эффект полностью вырождается. Во всех остальных случаях  > 1 и  z <  N . Отсюда следует, что при наличии коррелированной информации в массиве М положительный эффект, в той или иной мере, всегда имеет место.

3. Положительный эффект тем больше, чем больше значение x = m/n, меньше флюктуации х (величина  х), и меньше значения отсчетов N = n. Положительный эффект увеличивается именно в тех случаях, когда особенно остро ощущается недостаток информации: при малых значениях плотности потока излучения и/или экспозиции измерений.

Аналогичный эффект будет иметь место и при формировании отсчетов M по окрестностям текущих точек обработки данных путем определения их среднего значения (низкочастотное сглаживание массива n). Предварительное низкочастотное сглаживание может применяться и для статистически независимого дополнительного массива m, что будет повышать достоверность прогнозных отсчетов и увеличивать глубину регуляризации, если это сглаживание при регуляризации по формулам (11.2.20 и 21) не сказывается на изменении формы основного сигнала. Последнее определяется соотношением частотных спектров основного сигнала и оператора сглаживания.

Возможны два способа реализации уравнения (11.2.20): непосредственно в процессе измерений методом статистической группировки полезной информации (СГПИ) в реальном масштабе времени, или методом статистической регуляризации данных (СРД), зарегистрированных в виде временного (пространственного) распределения в параллельных массивах отсчетов.

^ 11.3. Статистическая регуляризация данных.

Как следует из выражения (11.2.21), для практического использования информации из дополнительных потоков данных необходимо установить значения и диспер-сию D(M), причем, исходя из задания последней по выражению (11.2.11), должно быть известно значение  x - относительной средней квадратической флюктуации вели-чины х.

Применительно к СРД определение значений и  x по зарегистрированным мас-сивам данных не представляет затруднений как в целом по пространству измерений, так и в виде распределений в скользящем окне усреднения данных. Последнее эквивалентно приведе-нию D xm => 0 для текущей точки обработки данных по информации ее ближайших окрестно-стей и позволяет производить максимальное извлечение полезной информации из дополни-тельных потоков сигналов, если частотный спектр распределения величины х по пространству измерений много меньше частотного спектра полезного сигнала. Отметим, что информация о распределении х также может иметь практическое значение (в частности, при гамма-опробовании с дополнительным потоком сигналов в низкоэнергетическом диапазоне спектра излучения - для оценки эффективного атомного номера горных пород).

Проверка теоретических положений метода CРД проводилась путем статистического моделирования соответствующих массивов данных и их обработки цифровыми фильтрами.

В таблице 1 приведены 4 группы результатов обработки по формулам (11.2.20-21) двух статистически независимых и постоянных по средним значениям массивов данных n и m (модели постоянных полей) при различных установках СРД по скользящему окну К с счета текущих значений = m i /n i и D i (М) по массиву m. Текущая точка обработки данных – по центру окна. Количество отсчетов в каждом массиве – 1000, распределение значений отсчетов соответствует закону Пуассона. Определение прогнозных отсчетов М i по массиву m для использования в уравнении (11.2.20) проводилось со сглаживанием отсчетов в скользящем окне K s низкочастотного цифрового фильтра (вариант без сглаживания при K s = 1). В качестве низкочастотного фильтра в алгоритме СРД используется (здесь и в дальнейшем) весовое окно Лапласа-Гаусса. Теоретическое значение D z.т. дисперсии результатов z определялось по выражению (11.2.22) с расчетом дисперсии D(M) по выражению D(M) =
. При сглаживании прогнозных отсчетов значение D M в выражении (11.2.22) принималось равным D M . = H s , где H s – коэффициент усиления сглаживающим фильтром дисперсии шумов (сумма квадратов коэффициентов цифрового фильтра). Дополнительно в таблице приводятся зарегистрированные средние значения коэффициента снижения статистических флюктуаций  =  n 2 / z 2 .

Таблица 1. Статистика результатов моделирования СРД.

(Основной массив = 9.9, D n = 9.7, дополнительный массив = 9.9, D m = 9.9, 1000 отсчетов.)

K c	K s	z	D z	Dz.т.		K c	K s	z	D z	Dz.т.	
3	1	9,7	5,7	6,19	1,7	11	3	9,6	3,6	3,80	2,8
5	1	9,7	5,4	5,78	1,8	11	5	9,6	3,3	3,55	3,0
11	1	9,6	5,1	5,36	1,9	11	11	9,6	3,1	3,22	3,2
21	1	9,6	5,0	5,18	2,0	11	21	9,6	3,0	3,11	3,3
51	1	9,6	5,0	5,05	2,0	11	51	9,6	3,0	2,99	3,3
3	3	9,7	4,1	4,71	2,4	3	11	9,8	4,5	4,26	2,2
5	5	9,7	3,6	4,01	2,8	5	11	9,7	3,5	3,78	2,8
11	11	9,6	3,1	3,22	3,2	11	11	9,6	3,1	3,22	3,2
21	21	9,6	2,9	2,91	3,4	21	11	9,6	3,1	3,12	3,2
51	51	9,6	2,7	2,66	3,7	51	11	9,6	3,1	2,99	3,2

Как видно из данных таблицы, практические результаты фильтрации достаточно хорошо совпадают с ожидаемыми по данным теоретических расчетов. Некоторое уменьшение среднего значения z по отношению к исходному среднему значению n определяется асимметричностью пуассоновского типа модели. При малых средних значениях модельных отсчетов в массиве m это приводит к определенной статистической асимметрии в работе СРД, т.к. при (+ m) 2 > (- m) 2 среднестатистическое доверие к дополнительной информации с отсчетами M i + меньше, чем с отсчетами M i -. Этим же фактором, по-видимому, вызвано и большее расхождение между теоретическими и фактическими значениями D z при малых значениях окна К с. Можно также заметить, что по значению коэффициента  фильтрация выходит на теоретические значения ( 1+MN) только при достаточно точном определении значений и D i (М), что требует увеличения окна К с счета этих параметров для полного использования дополнительной информации.

Таблица 2.

Эффект использования дополнительной информации, в полном соответствии с выражением (11.2.22), усиливается при предварительном сглаживании статистических вариаций отсчетов M i и при увеличении значений отсчетов дополнительного массива (материалы по последнему случаю не приводятся, т.к. не имеют какой-либо дополнительной информации). В спокойных по динамике полях еще большая глубина регуляризации может быть достигнута при счете значений и D m по сглаженному массиву М, что позволяет повысить вес прогнозных отсчетов M i . Результаты моделирования данного варианта в тех же условиях, что и для таблицы 1, приведены в таблице 2. Такой же эффект, в принципе, может достигаться и непосредственным введением дополнительного коэффициента веса в выражение (11.2.20) в качестве множителя для значения D(M), что позволяет осуществлять внешнее управление глубиной регуляризации.

Оценка сохранения разрешающей способности полезной информации была проведена на фильтрации детерминированных сигналов n и m предельной формы – в виде прямоугольных импульсов. Оценивались два фактора: сохранение формы полезного сигнала и подавление статистических шумов, наложенных на полезный сигнал.

При установке СРД без усреднения данных по массиву М (К s = 1, прогноз М i по текущим значениям массива М) при любых значениях окна К с выходной массив Z без всяких изменений повторяет массив N, т.е. не изменяет полезный сигнал и полностью сохраняет его частотные характеристики. Естественно, при условии, что массив М пропорционален массиву N.

При К s > 1 форма выходных кривых несколько изменяется и приведена на рис. 11.3.1. В индексах выходных кривых z приведена информация по установкам окон СРД: первая цифра - окно счета дисперсии D M и текущего значения (в количестве точек отсчетов), вторая цифра (через флеш) - окно сглаживания отсчетов М весовой функцией Лапласа-Гаусса и определения прогнозных отсчетов М i . Для сравнения с результатами типовой низкочастотной фильтрации на рисунке приведена кривая n25 отсчетов N, сглаженных весовой функцией Лапласа-Гаусса с окном 25 точек.

Рис. 11.3.1. СРД прямоугольного импульса. Счет D m по несглаженному массиву М.

На рис. 11.3.1а приведен результат СРД прямоугольного импульса с амплитудным значением 10 на фоне 5 при отношении m/n = 1 (равные значения отсчетов N и М). Дисперсия D N в выражении (11.2.21) принималась равной значению отсчетов N (статистика Пуассона). Как видно на рисунке, при сохранении фронтов сигнальной функции сглаживание прогнозных значений М i приводит к появлению искажения формы сигнала по обеим сторонам скачка, интервал которого тем больше, чем больше значение K s . Амплитудное значение искажений, как это и следует из выражения (11.2.21), в первую очередь зависит от соотношения текущих значений D N и D(M) и в меньшей степени от глубины сглаживания прогнозных отсчетов.

Максимальную величину искажения для точек скачка в первом приближении можно оценить из следующих соображений. Значения D(M) между точками скачка равны D(M) = А 2 /4, где А - амплитуда скачка, при этом значения коэффициента  для нижней и верхней точек скачка определяются выражениями   А 2 /(4D N +A 2), где D N = N точки скачка (для статистики Пуассона). Отсюда, при прогнозном значении М  N+А/2 для нижней точки скачка и M  N-A/2 для верхней точки относительная величина изменений N определится выражением   1/(2N/A+A), т.е. будет тем меньше, чем больше значения А и N и больше отношение N/A, что можно наглядно видеть на рис. 11.3.1в. Из этого выражения также следует, что максимальные искажения скачков, вносимые системой СРД, будут всегда в несколько раз меньше, чем статистические флюктуации непосредственных отсчетов  = 1/
на краях скачков.

При увеличении глубины регуляризации введением счета дисперсии D(M) по сглаженному массиву М картина искажений несколько изменяется и приведена на рис. 11.3.2. Реакция СРД на сглаживание дисперсии D(M) проявляется в своеобразной компенсации абсолютных отклонений отсчетов непосредственно по сторонам скачка отклонениями противоположного знака в более дальней зоне от скачка. Максимальные значения искажений остаются примерно на таком же уровне, как и для работы по несглаженной дисперсии D(M), с несколько меньшей зависимостью от увеличения значений N и А.

Рис. 11.3.2. СРД прямоугольного импульса. Счет D m по сглаженному массиву М.

В приведенных примерах значение окна счета К с принималось равным значению окна сглаживания К s дополнительного массива М. При К с > K s картина процесса практически не изменяется. При обратном соотношении размеров окон вступает в действие второй фактор - отклонение от фактических значений счета текущих значений x i = m/n в малом окне К с по массиву отсчетов, сглаженных с большим окном K s . На расстояниях от скачка функции, больших К с /2, СРД переходит в режим предпочтения сглаженных значений массива М, т.к. D(M)  0, что при К с < K s может приводить к появлению существенной погрешности – выбросов на расстояниях  К с /2 от скачков. Естественно, что при практических измерениях таких условий наблюдаться не будет и эффект резко уменьшится, но для полного его исключения вариант K c  K s можно считать предпочтительным.

Рис. 11.3.3. СРД сигнала N по массиву M. Рис. 11.3.4. Коэффициент .

(Счет D m по несглаженному массиву М). (Среднее статистическое по 50 циклам)

На рис. 11.3.3 приведен пример регистрации рандомизированного модельного сигнала в виде прямоугольного импульса амплитудой 40 на фоне 10, на котором виден принцип работы СРД. Как и следовало ожидать, СРД производит сглаживание статистических флюктуаций фона и сигнала за пределами зоны К с от скачка, отдавая предпочтение сглаженным прогнозным значениям М i , и не изменяет значения фона и сигнала в пределах этой зоны в связи с резким возрастанием текущих значений D(M) в выражении (11.3.21). Изменение коэффициента  в зоне скачка, управляющего формированием выходных отсчетов, приведено на рис. 11.3.4 (среднестатистическое по 50-ти циклам рандомизации для модельного импульса на рис. 11.3.3) и наглядно показывает принцип адаптации СРД к динамике изменения значений обрабатываемых сигналов.

Статистическая оценка регуляризации данных по прямоугольным импульсам проводилась по 50-ти циклам рандомизации исходных массивов N и M. В качестве примера на рисунках 11.3.5 и 6 приведены результаты обработки статистики массивов N и Z. Кроме статистики циклов рандомизации проводилась суммарная обработка всех циклов по общей статистике фона и вершины импульсов. Результаты обработки для тех же установок фильтров приведены в таблице 3.

Рис. 11.3.5. Статистика сигнала N Рис. 11.3.6. Статистика сигнала Z

(Измерения по 50-ти циклам). (50 циклов. Счет D m по несглаженному М)

Таблица 3.

Статистика значений фона и вершины импульсов (50 циклов).

Результаты моделирования подтверждают преимущество СРД перед простыми методами сглаживания. В числовой форме это наглядно проявляется в снижении дисперсии отсчетов выходного массива Z при практическом сохранении средних значений массива N и для фоновых отсчетов, и для амплитудных значений сигнала. При простом сглаживании "развал" фронтов сигнала (подавление высокочастотных составляющих спектра сигнала), как и должно быть при использовании низкочастотных фильтров, вызывает снижение по отношению к исходному массиву средних значений в максимумах и повышение фоновых значений сигнала, которое тем больше, чем больше окно весовой функции. Этот эффект особенно отчетливо проявляется в интервале окна фильтра по обе стороны от резких изменений сигнала.

При отсутствии дополнительных массивов М, коррелированных с регуляризируемым массивом N, формирование прогнозных значений М i может производиться по ближайшим окрестностям текущих значений N i в скользящем окне K s . При строго корректном подходе текущая точка N i не должна включаться в число счета прогнозных значений M i , но, как показало моделирование, это практически не влияет на результаты регуляризации. При прогнозировании M i по всем точкам окна K s массив М формируется любым методом сглаживания из массива N, и все особенности работы СРД по сглаженным массивам М, рассмотренные выше, остаются без изменений при условии счета значений D m в окне К с по массиву М. Для исключения выбросов по обе стороны от скачков полезного сигнала счет D m как дисперсии прогнозных значений M i необходимо выполнять непосредственно по массиву N.

Фундаментальной особенностью СРД является возможность последовательной многократной фильтрации данных, при которой может осуществляться преимущественное повышение степени регуляризации данных с минимальными искажениями формы полезного сигнала. Для выполнения последнего размер окна К с счета x i и D m устанавливается минимальным (3-5 точек), а глубина регуляризации данных (степень подавления шумов) устанавливается количеством последовательных операций фильтрации (до 3-5 проходов). Пример регуляризации модельного массива N в три прохода приведен на рис. 11.3.7.

Рис. 11.3.7. СРД одиночного массива N (3 прохода. Счет D m по массиву n)

Для сравнения пунктиром на рисунке приведено сглаживание массива 5-ти точечным фильтром Лапласа-Гаусса, который имеет коэффициент подавления шумов, эквивалентный 3-х проходному СРД (см. рис. 11.3.9).

На рисунках 11.3.8 и 11.3.9 приведены результаты статистической обработки 3-х проходной СРД для 25 циклов моделирования в сравнении с 1-м проходом и с 5-ти точечным фильтром Лапласа-Гаусса (кривая n5).

Рис. 11.3.8. Статистика средних значений Рис. 11.3.9. Статистика дисперсий

(25 циклов. Счет D m по массиву n) (25 циклов. Счет D m по массиву n)

Количество проходов может ограничиваться в автоматическом режиме, например, по среднеквадратическому значению корректирующих отсчетов z i = N i - z i в каждом проходе по сравнению с предыдущим проходом, которое сначала резко уменьшается за счет сглаживания флюктуаций, а затем, в зависимости от динамики сигнальной функции, стабилизируется или даже начинает увеличиваться за счет искажения самого сигнала.

Частотное представление работы СРД хорошо видно на рис. 11.3.10, где приведены модули спектров рандомизированного сигнала в виде меандра (средние значения в минимуме - 20, в максимуме - 100, 25 периодов по 40 отсчетов, всего 1000 отсчетов) и результатов его обработки СРД (окно К с = 3, окно К s = 3).

Рис. 11.3.10. Модули спектров модельных сигналов. Рис.11.3.11. Участок спектра.

(1– входной массив N, 2– выходной массив Z, один цикл CРД,

3– выходной массив Z, три цикла CРД), 4 – массив нерандомизированного меандра).

Модуль спектра основного полезного сигнала (в данном случае чистого меандра) представляет собой последовательность отдельных частотных гармоник по всему диапазону спектра. В спектре рандомизированного меандра эти частотные гармоники суммируются со спектром шума, статистически равномерно распределенным по всему частотному диапазону (спектр шума на рисунке для наглядности сглажен). СРД осуществляет подавление шумовых составляющих сигнала, практически не затрагивая частотных гармоник меандра и не изменяя их по амплитуде. Последнее можно видеть на рис. 11.3.11, где представлен отрезок спектра сигналов в высокочастотной части главного диапазона в области одной гармоники меандра (частотные составляющие шума не сглажены). При 3-х цикловом СРД высокочастотные составляющие шумов подавляются практически на порядок.

Пример практического использования СРД приведен на рис. 11.3.12 при опробовании участка скважины, пересекающей пласты каменной соли, на содержание сильвинита по гамма-излучению Калия-40. По данным геологического опробования пласты сильвинита в толще вмещающих пород (галита) имеют достаточно резкие границы и однородны по содержанию сильвинита в пределах пластов. Исходная диаграмма ГК (детектор CsJ(Tl) со свинцовым фильтром толщиной 2 мм) и результаты фильтрации исходного массива данных ГК с использованием СРД и низкочастотного фильтра с весовым окном Лапласа-Гаусса приведены на рис. 11.3.12.

Рис. 11.3.12. Диаграммы ГК.

Результаты интерпретации диаграмм ГК симметричным деконволюционным цифровым фильтром (окно 13 точек) приведены на рис. 11.3.13. Как видно на рисунке, деконволюция по несглаженной диаграмме ГК дает существенные вариации содержания сильвинита в пределах пластов. Применение низкочастотной фильтрации диаграммы ГК снимает флюктуации содержания в пределах пластов, но существенно сглаживает границы пластов. Использование СРД позволяет устранить этот недостаток.

Рис. 11.3.13. Результаты интерпретации диаграмм ГК.

В заключение отметим, что СРД может использоваться для регуляризации не только ядернофизических данных, но и любых других числовых массивов непрерывных измерений, если радиус их корреляции не менее 3-5 отсчетов. В качестве примера на рис. 11.3.14 приведена диаграмма акустического каротажа, зарегистрированная с шагом дискретизации данных 20 см, сглаживание которой проведено СРД без потери пространственного разрешения.

Рис. 11.3.14. Диаграмма акустического каротажа и результат ее обработки СРД

(5 циклов, K c = K s = 3, физическое окно 0.6 м).

Курсовая работа 17-07. Модернизация адаптивного фильтра сглаживания данных, статистических распределенных по закону Пуассона.

^ 11.3. Статистическая группировка полезной информации.

Что касается аппаратных способов реализации СГПИ, то он может быть выполнен в реальном масштабе времени, если информация представлена потоком импульсов и основным информативным параметром является скорость следования импульсов.

Сущность аппаратной реализации заключается в статистической (близкой к статистической) нормированной выборке импульсов из дополнительного потока m и их суммировании с основным потоком n с заданием условий выборки по отношению частоты следования импульсов в потоках. Пола-гая для непрерывного режима измерений M+1 = М, перепишем выражение (5.2.20) с подстановкой значения  в следующем виде:

Z = N + (M/-N)·M/(M+D(M)). (11.3.1)

Умножим левую и правую части выражения на нормировочный коэффициент размножения выходного потока K = l+R:

Z = K·z= N + RN+(M/-N)·KM/(M+D(M). (11.3.2)

Заменим отсчеты RN выборкой сигналов из потока m:

RN = Р в М, (11.3.3)

Где Р в - вероятность выборки сигналов из потока m. Если вероятность выборки сигналов под-держивать равной значению

P в = R/, (11.3.4)

То при этом будет иметь место

M/-N = Р в M/R-N  0, (11.3.5)

И соответственно для выражения (11.3.2) имеем:

(M/-N)·KM/(M+D(M)  0, (11.3.6)

Z = N+P в M  N+RN. (11.3.7)

При статистической независимости величины х от частоты потоков n и m приведен-ные выражения действительны при определении значения как в целом по пространству измерений, так и для скользящих окон текущих значений по определенным интервалам предшествующих измерений. Действительно и обратное заключение: если по определенному интервалу измере-ний выражение (11.3.5) обращается в нуль, то установленная вероятность выборки соответствует условию (11.3.4). На этом принципе может проводиться аппаратная реализация СГПИ с авто-матической адаптацией к условиям измерений: управление процессом выборки импульсов из потока m и направление их на суммирование с потоком n по сигналам обратной связи с уст-ройства, следящего за обращением в нуль выражения (11.3.5).

Особенности аппаратной реализации СГПИ с автоматической адаптацией под условия измерений заключаются в следующем.

Значение вероятности выборки Р в не может быть больше 1. Отсюда из (11.3.3) следует, что для любых интервалов измерений должно выполняться условие М ≥ RN, а соответственно по всему пространству измерений должно выполняться условие ≥ R, чем и обуславливается выбор коэффициента R. Значение коэффициента R принципиально ограничивает степень положительного эффекта СГПИ (k max  1+R), в отличие от СРД, где такого ограничения не имеется.

Относительная статистическая погрешность измерений выходного потока отсчетов Z соответ-ствует выражению (11.2.23) при условии постоянного значения величины Р в, т.е. при установке значения Р в по среднему значению величины в целом по пространству измерений. При автоматической адаптации под условия измерений значение вероятности Р в по текущему среднему значению отношения n/m определенного предшествующего интервала измерений также является статистически флюктуирующей величиной с дисперсией распределения (без учета изменений действительного значения х):

D p = R 2 (n+m)n/(m 3 T), (11.3.8)

Где Т- интервал усреднения информации при определении текущего значения . Соответственно, дисперсия и средняя квадратическая погрешность текущих отсчетов Z:

D z = D N + P в D M +M 2 D p = N+Р в М+М 2 D р, (11.3.9)

 z 2 = (N+Р в М+М 2 D р)/(N+Р в М) 2 . (11.3.10)

При постоянной экспозиции измерений  положительный эффект возрастает с увеличением значения Т:

K = K 2 /(K+R 2 (n+m)/mT). (11.3.11)

K max  1+R,  z 2  1/(N+Р в М) при Т  . (11.3.12)

В общем случае, с учетом средней квадратической ошибки прогнозирования  xi значе-ний x i для текущих точек измерений по значениям в предшествующих интервалах при Т > :

D z = N+Р в М+M 2 (D p +P в 2  xi 2). (11.3.13)

Формирование значения Р в на основе информации по средним значениям интер-валов измерений, предшествующих текущим, определяет СГПИ как динамическую систему с соответствующей постоянной времени реакции на изменение условий измерений. Учитывая, что, во-первых, для любой точки пространства измерений должно выполняться условие m > nR, и, во-вторых, увеличение интервала Т приводит к возрастанию времени реакции на из-менение условий измерений, значение Т целесообразно ограничивать величиной порядка (5-10) значений текущих экспозиций. Чем меньше пространственная частота рас-пределения х по отношению к распределению n, тем большее значение Т допустимо.

Реализация систем СГПИ значительно облегчается при чисто практическом ограниче-нии целевой задачи: получение максимального положительного эффекта в экстремально небла-гоприятных условиях производства измерений (при низких значениях регистрируемой плот-ности потока излучения, при высокой скорости измерений) с вырождением положительного эффекта по мере снижения статистической погрешности измерений в основном потоке. Так, например, если при скважинном гамма-опробовании статистическая погрешность измерений основного потока сигналов в зонах с повышенной интенсивностью излучения снижается до 2-3%, то ее дальнейшее уменьшение не имеет практического смысла, т.к. основная погрешность каротажной радиометрической аппаратуры обычно не превышает 5%.

Использование данного целевого ограничения позволяет применить формирование параметра Р в не в скользящем окне временного или пространственного усреднения информа-ции, а по определенному зарегистрированному объему предшествующей информации, т.е. с автоматической вариацией интервала усреднения информации и постоянной регулирования P в в зависимости от частоты потоков сигналов, при этом объем информации формирования P в может задаваться с учетом характера вариаций величины и допустимого значения ди-намической погрешности измерений.

Для реализации такой возможности преобразуем выражение (11.3.5) по интервалу усред-нения t к виду:

P в mt/R-nt+Q = q, (11.3.14)

P в = nR/m = q/, (11.3.15)

Q  Q при t  ,

Где Q- средний уровень смещения числового эквивалента сигнала обратной связи системы АРВ - автоматического регулирования вероятности выборки Р в, при котором обеспечивается вы-полнение равенства (11.3.15), - коэффициент пропорциональности преобразования цифрового сигнала АРВ в сигнал Р в. Дифференциальное уравнение для системы АРВ:

Dq/dt = n-mq/R. (11.3.16)

Решение дифференциального уравнения при начальных условиях t = 0 и q = О (переходная функция АРВ):

Q = R(n/m) . (11.3.17)

P в = R(n/m) = R(n/m) . (11.3.18)

Как видно из этих выражений, значение сигнала обратной связи АРВ пропорционально отношению (n/m) частот потоков, а постоянная времени АРВ R/m прямо пропорцио-нальна значению коэффициента преобразования  при обратной пропорциональности от зна-чения частоты дополнительного потока m, равно как и, с учетом (11.3.15), прямо пропорциональ-на текущему значению сигнала обратной связи q при обратной пропорциональности от зна-чения частоты основного потока n. Первое полностью эквивалентно второму при (n/m)  const и q = Rn/m  Q. В первом приближении, с использованием выражения (11.3.8) и экви-валентности значения статистических флюктуаций при Т≈2 для скользящего прямоугольного временного окна и окна интенсиметра с экспоненциальной переходной функцией, для отно-сительных флюктуации значения Р в получаем:

 р 2 = (n+m)/(2Rn)= (n+m)/(2qm). (11.3.19)

Выражение действительно для прямого измерения 2-интенсиметром отношения (n/m) и является максимальной оценкой. Для более точной оценки следует учитывать, что в данном случае интенсиметр является устройством с отрицательной обратной связью по цепи АРВ, что несколько уменьшает значение флюктуации. Точная оценка может быть проведена с использованием формулы Кэмпбелла для дисперсии случайной величины x(t), образованной сложением импульсов пуассоновского потока , раздельно для потока n при m = const и потока m при n = const, с последующим сложением квадратов относительного среднего квадратического значения флюктуации. Так, для схемы, приведенной ниже, получено значение  р 2 ≈ (R+1)m/(2nR 2).

При выбранном для пространства измерений значении коэффициента R ≤ (m/n) min с использованием выражения (11.3.19) параметры системы АРВ (коэффициент  и среднее значе-ние Q для средней по пространству величины отношения n/m) могут устанавливаться под заданное значение допустимых флюктуаций вероятности выборки импульсов Р в:

 ≤ (l+(m/n) max)/(2R p 2). (11.3.20)

В процессе измерений АРВ осуществляет непрерывную адаптацию под текущие усло-вия измерений (nq, m  mR, P в  q/) с регулированием текущего значения P в по объему информации q = (n/m)R = n предшествующего интервала измерений путем соответ-ствующего изменения постоянной времени интегрирования этой информации в зависимости от изменения частот потоков сигналов. При n/m  const последнее имеет абсолютный ха-рактер:  р  const,   (l/n + l/m)/(2 p 2).

Следует отметить, что во многих методах геофизики существуют достаточно благоприятные условия использования как СГПИ, так и СРД. Так, например, применительно к скважинному гамма-опробованию с извлечением дополнительной информации из низкоэнергетической час-ти спектра излучения условия достаточно точной реакции на изменения параметра по стволу сква-жины являются весьма хорошими, т.к. основной фактор вариации значений x - эф-фективный атомный номер среды, изменяется в небольшом диапазоне с низкой простран-ственной частотой вариаций, причем в зонах расположения активных пород, где требуется наиболее высокая точность интерпретации результатов измерений и возможны значительные изменения атомного номера пород, за счет увеличения плотностей потоков излучения посто-янная времени АРВ будет существенно уменьшаться, а пространственная разрешающая спо-собность измерений соответственно увеличиваться. Аналогичные условия характерны, как пра-вило, и для других методов ядерной геофизики.

Пример исполнения системы СГПИ для двух импульсных потоков сигналов приведен на рис. 11.3.1. Функциональная схема СГПИ содержит реверсивный счетчик импульсов 1, на вход суммирования которого подаются импульсы основ-ного потока n, а на вход вычитания - импульсы дополнительного потока m, предварительно проходящие через схему выборки импульсов 3 и счетчик-делитель частоты следования импуль-сов 4 с коэффициентом пересчета R.

Рис. 11.3.1. Базовая функциональная схема СГПИ.

1- реверсивный счетчик импульсов, 2- блок формирования сигнала выборки импульсов, 3- схема выборки импульсов, 4- счетчик-делитель частоты на R, 5- блок суммирования потоков импульсов.
Информация о состоянии счетчика 1 (сигнал q) с выходов счетчика подается на блок формирования сигнала выборки импульсов 3. В простейшем случае этот блок может представлять собой пороговое устройство (по коду числа Q), открывающее схему 3, однако выборка в этом случае имеет характер, близкий к статистическому, только при достаточно малых различиях частоты потоков n и m/R (порядка n

Импульсы основного потока n и импульсы выборки из потока m, частота которых равна Р в m = R·n, поступают на вход блока 5 суммирования потоков сигналов. Интенсивность потока импульсов на выходе блока 5 равна z = n+Р в m = (1+R)n. Блок 5 может содержать пе-ресчетную схему с коэффициентом K=(1+R), при этом выходной поток будет приводиться к масштабу основного потока n и появляется возможность синхронного переключения коэффи-циентов пересчета схем 4 и 5 под различные условия измерений, при этом установка опти-мального значения коэффициента R может быть переведена в режим автоматической с управ-лением по текущему значению (в определенном интервале) информационного кода схемы 1. Альтернативное решение - подача на вход суммирования схемы 5 потока импульсов с выхода схемы 4, при этом частота потока z будет всегда в 2 раза больше потока n.

Попутно отметим, что при выводе информации q = R(n/m) в цифровом коде со счет-чика 1 данная схема может выполнять функции универсального цифрового интенсиметра: средней частоты импульсов (n-var, m-const от генератора тактовой частоты), среднего вре-менного интервала между импульсами (m-var, n-const) и отношения частот n/m двух статис-тически распределенных потоков импульсов.

литература

38. Адаптивные фильтры. /Под ред. К.Ф.Н.Коуэна и П.М.Гранта. – М.: Мир, 1988, 392 с.

43. Айфичер Э., Джервис Б. Цифровая обработка сигналов. Практический подход. / М., "Вильямс", 2004, 992 с.

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

(ФГБОУ ВПО)

«ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ»

ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ, ИНФОРМАТИКИ И МЕХАНИКИ

КАФЕДРА ТЕХНИЧЕСКОЙ КИБЕРНЕТИКИ И АВТОМАТИЧЕСКОГО РЕГУЛИРОВАНИЯ

Применение адаптивных фильтров для идентификации систем

Курсовая работа

Выполнил:

студент 5-го курса в/о

Литикова А.С.

Научный руководитель:

д.ф.-м.н. Аверина Л.И.

Воронеж 2014

Введение

2.1 Оптимальный фильтр Винера

2.2 Адаптивный алгоритм LMS

2.4 Адаптивный алгоритм RLS

4. Реализация моделей адаптивных фильтров, сравнение результатов

Приложения

Библиографический список

Введение

Цифровая обработка сигналов в сложных радиоэлектронных системах характеризуются рядом особенностей, осложняющих решение радиотехнических задач. В высококачественных системах цифровой обработка сигналов часто ведется в условиях неопределенности системных характеристик. Априорная и текущая информация о параметрах системы является неполной. Кроме того, при эксплуатации системы параметры её объектов и среды функционирования могут изменяться непредвиденным заранее образом и тогда нестационарность выступает как один из видов неопределенности. Стохастичность поведения также является важной чертой, характеризующей обработку сигналов в сложных системах. Это обусловлено наличием источников случайных помех и всякого рода второстепенных (с точки зрения решения задачи) процессов с непредсказуемым поведением. Таким образом, стохастичность и нестационарность объектов и условий их функционирования определяют фактор сложности системы.

Одним из наиболее перспективных путей преодоления трудностей, порождаемых этим фактором сложности является применение методов адаптации. При адаптивной цифровой обработки информация об объекте и внешних воздействиях собирается в ходе обработки и используется для изменения параметров блоков обработки. Адаптивные устройства обработки сигнала действуют по принципу замкнутого контура. Входной сигнал фильтруется или взвешивается в программируемом фильтре для получения выходного сигнала, который затем сравнивается с полезным, стандартным или обучающим сигналом для нахождения сигнала ошибки. Затем этот сигнал ошибки используется для корректировки весовых параметров процессора обработки с целью постепенной минимизации ошибки. При этом адаптивной системе требуется минимальный объект исходной информации о поступающем сигнале.

Процессоры радиотехнической обработки можно разделить на два больших класса: адаптивные фильтры и адаптивные антенны. Адаптивные антенны осуществляют пространственную обработку сигнала с помощью антенной решетки, формируя диаграмму направленности таким образом, чтобы в ней создавался главный максимум в направлении прихода сигнала и происходила генерации нулей в направлении источников помех. В адаптивных цифровых фильтрах используется программируемый фильтр, частотная характеристика или передаточная функция которого изменяется, или адаптируется, таким образом, чтобы пропустить без искажения полезные составляющие сигнала и ослабить нежелательные сигналы или помехи, т.е. уменьшить искажения входного полезного сигнала. Адаптивный фильтр действует по принципу оценивания статистических параметров поступающего сигнала и подстройки собственной переходной характеристики таким образом. чтобы минимизировать некоторую функцию стоимости. Такие адаптивные фильтры часто используются для восстановления на выходах каналов сигналов с меняющимися во времени характеристиками.

Цель настоящей курсовой работы - исследовать и изучить научную литературу по рассматриваемой теме, исследовать работу основных алгоритмов адаптации (LMS-алгоритма, RLS-алгоритма и RLS-алгоритма с экспоненциальным забыванием) на примере идентификации фильтра 31-го порядка с импульсной характеристикой в виде экспоненциально затухающей синусоиды. Кроме того, будет рассмотрена реакция алгоритмов на скачкообразные изменения параметров системы.

1. Основная идея адаптивной обработки сигнала

Общая структура адаптивного фильтра показана на рис. 1. Входной дискретный сигнал х(k) обрабатывается дискретным фильтром, в результате чего получается выходной сигнал у(k). Этот выходной сигнал сравнивается с образцовым сигналом d(k), разность между ними образует сигнал ошибки е(k), Задача адаптивного фильтра - минимизировать ошибку воспроизведения образцового сигнала. С этой целью блок адаптации после обработки каждого отсчета анализирует сигнал ошибки и дополнительные данные, поступающие из фильтра, используя результаты этого анализа для подстройки параметров (коэффициентов) фильтра.

Размещено на http://www.allbest.ru/

Рис. 1. Общая структура адаптивного фильтра.

Возможен и иной вариант адаптации, при котором образцовый сигнал не используется или требуемые параметры сигнала напрямую не могут быть измерены. Такой режим работы называется адаптацией по косвенным данным (blind adaptation) или обучением без учителя (unsupervised learning). Исследование последнего варианта адаптации не входит в рамки данной курсовой работы.

В качестве фильтра в структуре, показанной на рис. 1, чаще всего используется нерекурсивный цифровой фильтр. Одним из главных достоинств этого варианта является то, что нерекурсивный фильтр является устойчивым при любых значениях коэффициентов. Однако следует помнить, что алгоритм адаптации в любом случае вносит в систему обратную связь, вследствие чего адаптивная система в целом может стать неустойчивой.

Далее будут рассмотрены два адаптивных алгоритма с использованием образцового сигнала, часто применяемых на практике в различных системах обработки информации. Для упрощения математических выкладок предположим, что сигналы и фильтры являются вещественными. Однако результирующие формулы легко обобщаются на случай комплексных сигналов и фильтров.

2. Алгоритмы адаптивной фильтрации

2.1 Оптимальный фильтр Винера

Прежде чем рассматривать собственно алгоритмы адаптации, необходимо определить те оптимальные параметры фильтра, к которым эти алгоритмы должны стремиться. Подход к задаче оптимальной фильтрации может быть как статистическим, так и детерминированным. Сначала рассмотрим статистический вариант.

Пусть входной дискретный случайный сигнал {} обрабатывается нерекурсивным дискретным фильтром порядка N, с коэффициентами {wn},

n=0,1,…,n. Выходной сигнал фильтра равен

Кроме того, имеется образцовый (также случайный) сигнал d(k). Ошибка воспроизведения образцового сигнала равна

Необходимо найти такие коэффициенты фильтра {wn}, которые обеспечивают максимальную близость выходного сигнала фильтра к образцовому, то есть минимизируют ошибку. Поскольку также является случайным процессом, в качестве меры ее величины разумно принять средний квадрат. Таким образом, оптимизируемый функционал выглядит так

Для решения поставленной задачи перепишем (2) в матричном виде. Для этого обозначим вектор-столбец коэффициентов фильтра как w, а вектор-столбец содержимого линии задержки на k-м шаге как x(k):

Тогда (2) примет вид:

Квадрат ошибки равен

Статистически усредняя это выражение, получаем следующее:

адаптивный сигнал фильтр алгоритм

Входящие в полученную формулу усредненные величины имеют следующий смысл:

с - средний квадрат образцового сигнала, не зависит от коэффициентов фильтра, потому может быть отброшено;

с - вектор-столбец взаимных корреляций между k-м отсчетом образцового сигнала и содержимым линии задержки фильтра на k-м шаге. Будем считать случайные процессы x(k) и совместно стационарными, тогда вектор взаимных корреляций не зависит от номера шага k. Обозначим этот вектор как p:

с - корреляционная матрица сигнала, имеющая размер. Для стационарного случайного процесса корреляционная матрица имеет вид матрицы Теплица, то есть вдоль ее диагоналей стоят значения корреляционной функции:

где - корреляционная функция (КФ) входного сигнала.

С учетом введенных обозначений (4) принимает следующий вид:

Данное выражение представляет собой квадратичную форму относительно w и потому при невырожденной матрице R имеет единственный минимум, для нахождения которого необходимо приравнять нулю вектор градиента:

Отсюда получаем равенство:

Умножив слева обе части на обратную корреляционную матрицу, получим искомое решение для оптимальных коэффициентов фильтра:

Такой фильтр называется фильтром Винера. Подстановка (8) в (5) дает минимально достижимую дисперсию сигнала ошибки:

Несложно также показать, что и, то есть что сигнал ошибки для фильтра Винера не коррелирован со входным и выходным сигналами фильтра.

2.2 Адаптивный алгоритм LMS

Один из наиболее распространенных адаптивных алгоритмов основан на поиске минимума целевой функции методом наискорейшего спуска. При использовании данного способа оптимизации вектор коэффициентов фильтра w зависит от номера итерации k: , и должен на каждой итерации смещаться на величину, пропорциональную градиенту целевой функции в данной точке, следующим образом:

где - положительный коэффициент, называемый размером шага. Подробный анализ сходимости данного процесса приведен, например, в [фильтр Винера]. Показано, что алгоритм сходится, если

где - максимальное собственное число корреляционной матрицы R. Скорость сходимости при этом зависит от разброса собственных чисел корреляционной матрицы R- чем меньше отношение тем быстрее сходится итерационный процесс.

Однако для реализации алгоритма требуется вычислять значения градиента, а для этого необходимо знать значения матрицы R и вектора р. На практике могут быть доступны лишь оценки этих значений, получаемые по входным данным. Простейшими такими оценками являются мгновенные значения корреляционной матрицы и вектора взаимных корреляций, получаемые без какого-либо усреднения:

При использовании данных оценок формула принимает следующий вид:

Выражение, стоящее в скобках, согласно (3), представляет собой разность между образцовым сигналом и выходным сигналом фильтра на k-м шаге, то есть ошибку фильтрации е(k). С учетом этого, выражение для рекурсивного обновления коэффициентов фильтра оказывается очень простым:

Алгоритм адаптивной фильтрации, основанный на данной формуле, получил название LMS (Least Mean Square, метод наименьших квадратов). Можно получить ту же формулу и несколько иным образом: использовав вместо градиента статистически усредненного квадрата ошибки градиент его мгновенного значения.

Анализ сходимости алгоритма LMS оказывается чрезвычайно сложной задачей, для которой не существует точного аналитического решения. Достаточно подробный анализ с использованием приближений приведен в книге С. Хайкина. Он показывает, что верхняя граница для размера шага в данном случае является меньшей, чем при использовании истинных значений градиента. Эта граница примерно равна:

где - собственные числа корреляционной матрицы R, а - средний квадрат входного сигнала фильтра.

На вышеприведенной формуле основан нормированный (normalized) LMS-алгоритм, в котором коэффициент на каждом шаге рассчитывается, исходя из энергии сигнала, содержащегося в линии задержки:

где µ0 - нормированное значение µ, лежащее в диапазоне от 0 до 2, а е -- малая положительная константа, назначение которой -- ограничить рост µ при нулевом сигнале на входе фильтра.

Даже если LMS-алгоритм сходится, дисперсии коэффициентов фильтра при k>? не стремится к нулю -- коэффициенты флуктуируют вокруг оптимальных значений. Из-за этого ошибка фильтрации в установившемся режиме оказывается больше, чем ошибка для винеровского фильтра:

где Eex - средний квадрат избыточной ошибки (excess error) алгоритма LMS.

В той же книге С.Хайкина приводится следующая приближенная формула для так называемого коэффициента расстройства (misalignment), равного отношению средних квадратов избыточной и винеровской ошибок:

Значение коэффициента µ влияет на два главных параметра LMS-фильтра: скорость сходимости и коэффициент расстройки. Чем больше µ, тем быстрее сходится алгоритм, но тем больше становится коэффициент расстройки и наоборот.

Основным достоинством алгоритма LMS является предельная вычислительная простота - для подстройки коэффициентов фильтра на каждом шаге нужно выполнить N + 1 пар операций «умножение-сложение». Платой за простоту является медленная сходимость и повышенная (по сравнению с минимально достижимым значением) дисперсия ошибки в установившемся режиме.

2.3 Детерминированная задача оптимальной фильтрации

При рассмотрении статистической задачи оптимизации входной сигнал считался случайным процессом и минимизировалась дисперсия ошибки воспроизведения образцового сигнала. Однако возможен и иной подход, не использующий статистические методы.

Пусть, как и раньше, обработке подвергается последовательность отсчетов {x(k)}, коэффициенты нерекурсивного фильтра порядка N образуют набор {wn}, а отсчеты образцового сигнала равны {d(k)}. Выходной сигнал фильтра определяется формулой (1), а ошибка воспроизведения образцового сигнала - формулой (2) или, в векторном виде, (3).

Теперь оптимизационная задача формулируется так: нужно отыскать такие коэффициенты фильтра {wn}, чтобы суммарная квадратичная ошибка воспроизведения образцового сигнала была минимальной:

Для решения задачи необходимо перейти в формуле (3) к матричной записи вдоль координаты k, получив формулы для векторов-столбцов выходного сигнала у и для ошибки воспроизведения входного сигнала е:

Здесь d - вектор-столбец отсчетов образцового сигнала, а X - матрица, столбцы которой представляют собой содержимое линии задержки фильтра на разных тактах:

Выражение (18) для для суммарной квадратичной ошибки можно переписать в матричном виде следующим образом:

Подставив (19) в (20), имеем:

Для нахождения минимума необходимо вычислить градиент данного функционала и приравнять его нулю:

Отсюда легко получается искомое оптимальное решение:

В формуле прослеживается близкое родство с формулой (8), описывающей оптимальный в статистическом смысле фильтр Винера. Действительно, если учесть, что дает оценку корреляционной матрицы сигнала, полученную по одной реализации эргодического случайного процесса путем временного усреднения, а является аналогичной оценкой взаимных корреляций между образцовым сигналом и содержимым линии задержки фильтра, то формулы (8) и (21) совпадут.

2.4 Адаптивный алгоритм RLS

В принципе, в процессе приема сигнала можно на каждом очередном шаге пересчитывать коэффициенты фильтра непосредственно по формуле (21), однако это связано с неоправданно большими вычислительными затратами. Действительно, размер матрицы X постоянно увеличивается и, кроме того, необходимо каждый раз заново вычислять обратную матрицу.

Сократить вычислительные затраты можно, если заметить, что на каждом шаге к матрице X добавляется лишь один новый столбец, а к вектору d - один новый элемент. Это дает возможность организовать вычисления рекурсивно. Соответствующий алгоритм называется рекурсивным методом наименьших квадратов (Recursive Least Square, RLS).

На k-м шаге оптимальный вектор коэффициентов фильтра, согласно (21), равен

К матрице X на каждом шаге добавляется новый столбец x(k+1), а к вектору d - новый элемент d(k+1) :

При использовании алгоритма RLS производится рекурсивное обновление оценки обратной корреляционной матрицы.

При переходе к следующему шагу имеем:

Теперь воспользуемся матричным тождеством

где А и С - произвольные квадратные невырожденные матрицы, а В и D - произвольные матрицы совместимых размеров.

Установим между (24) и (25) соответствие:

с (квадратная матрица);

с (вектор-столбец);

с (скаляр);

с (вектор-строка).

В результате (24) можно записать следующим образом:

Заметим, что фрагмент выражения, возводимый в минус первую степень, представляет собой скаляр. Используем (23) и (26) в выражении (22) для коэффициентов фильтра:

Раскроем скобки:

Первое слагаемое в полученной формуле, согласно (22), представляет собой коэффициенты оптимального фильтра для k-го шага -- w(k). Этот же вектор можно выделить в качестве множителя во втором слагаемом. В третьем и четвертом слагаемых можно выделить общий множитель P(k)x(k+1)d(k+1):

Выполнив вычисление в скобках и вынеся общий множитель, получаем:

Заметим, что произведение есть результат обработки поступившего на k-м шаге входного сигнала фильтром со старыми коэффициентами w(k). То есть, это произведение представляет собой выходной сигнал адаптивного фильтра y(k+1). Соответственно, разность в скобках -- это ошибка фильтрации e(k+1):

Теперь введем обозначение для вынесенного за скобки векторного множителя:

С учетом этого формула для коэффициентов фильтра примет вид:

Вектор K(k+1) называют вектором коэффициентов усиления.

Итак, при использовании адаптивного алгоритма RLS необходимо на на каждом такте выполнять следующие шаги:

1. При поступлении новых входных данных x(k) производится фильтрация сигнала с использованием текущих коэффициентов фильтра w(k-1) и вычисление величины ошибки воспроизведения образцового сигнала:

2. Рассчитывается вектор-столбец коэффициентов усиления (следует отметить, что вектор K всякий раз рассчитывается заново, т. е. Вычисления не рекурсивны, а также что знаменатель дроби является скаляром):

3. Производится обновление оценки обратной корреляционной матрицы сигнала:

Наконец, производится обновление коэффициентов фильтра:

Что касается рекурсивно обновляемых матрицы P и вектора w, обычно вектор w принимается нулевым, а анализ матрицы Р показывает, что после заполнения линии задержки фильтра отсчетами сигнала результат не будет зависеть от начальных условий, если

На практике диагональ заполняется большими положительными числами, например, 100/у2x .

По сравнению с LMS-алгоритмом, алгоритм RLS требует значительно большего числа вычислительных операций (при оптимальной организации вычислений 2,5N2+4N пар операций «умножение-сложение»). Под оптимальной организацией здесь понимается учет симметрии матрицы P. Таким образом, число операций в алгоритме RLS квадратично возрастает с увеличением порядка фильтра.

Зато алгоритм RLS сходится значительно быстрее алгоритма LMS. Строго говоря, он даже не является алгоритмом последовательного приближения, поскольку на каждом шаге дает оптимальные коэффициенты фильтра, соответствующие формуле (21).

2.5 Алгоритм RLS с экспоненциальным забыванием

В формулах (18) и (20) значениям ошибки на всех временных тактах придается одинаковый вес. В результате, если статистические свойства входного сигнала со временем изменяются, это приведет к ухудшению качества фильтрации. Чтобы дать фильтру возможность отслеживать нестационарный входной сигнал, можно применить в (18) экспоненциальное забывание (exponential forgetting), при котором вес прошлых значений сигнала ошибки экспоненциально уменьшается:

где л -- коэффициент забывания,

При использовании экспоненциального забывания формулы (28) и (29) принимают следующий вид:

3. Адаптивные фильтры в идентификации систем

Адаптивные фильтры в настоящее время нашли применение во многих радиотехнических и телекоммуникационных системах обработки сигналов. Однако все способы использования адаптивных фильтров так или иначе сводятся к решению задачи идентификации, то есть определения характеристик некоторой системы.

Возможны два варианта идентификации - прямая и обратная. В первом случае адаптивный фильтр включается параллельно с исследуемой системой (рис. 2, а).

Размещено на http://www.allbest.ru/

Рис.2 Идентификация систем с помощью адаптивного фильтра:

а- прямая, б- обратная.

Входной сигнал является общим для исследуемой системы и адаптивного фильтра, а выходной сигнал системы служит для адаптивного фильтра образцовым сигналом. В процессе адаптации временные и частотные характеристики фильтра будут стремиться к соответствующим характеристикам исследуемой системы.

При обратной идентификации адаптивный фильтр включается последовательно с исследуемой системой (рис. 2, б). Выходной сигнал системы поступает на вход адаптивного фильтра, а входной сигнал системы является образцом для адаптивного фильтра. Таким образом, фильтр стремится компенсировать влияние системы и восстановить исходный сигнал, устранив внесенные системой искажения.

4. Реализация моделей адаптивных фильтров

В данной работе мы реализуем решение задачи прямой идентификации системы согласно рис.2.а с помощью трех алгоритмов - LMS, RLS и RLS с экспоненциальным забыванием.

Входным сигналом будет служить дискретный белый гауссов шум, а сама система будет представлять собой нерекурсивный фильтр 31-го порядка с импульсной характеристикой в виде экспоненциально затухающей синусоиды, на период колебаний которой приходится четыре отсчета. Кроме того, после обработки половины входного сигнала параметры системы скачкообразно изменятся -- у ее импульсной характеристики поменяется знак. Это позволит нам посмотреть на то, как реагируют алгоритмы на резкое изменение статистических свойств обрабатываемого сигнала.

Текст программы, реализующей алгоритмы, см. в приложении 1.

Сравнение моделей

Итак, на рис.3 можно видеть, что на начальном этапе RLS-алгоритмы показывают гораздо более быстрый переход и меньшую ошибку по сравнению с LMS-алгоритмом. Теперь рассмотрим подробнее ошибку от времени на начальном этапе (Рис.4.).

На рис.4 видно, что LMS и RLS алгоритмы дают примерно одинаковую ошибку в установившемся режиме, тогда как алгоритм RLS с экспоненциальным забыванием ошибки практически не дает.

Вернемся к рис.3 и обратимся к скачкообразному изменению характеристик системы. Здесь мы видим, что LMS-алгоритм после некоторого переходного процесса свел ошибку фильтрации к тому же уровню, который она имела раньше, а вот RLS так и не смог приспособиться к изменившимся условиям.

Рис.3. Результаты идентификации системы с помощью алгоритмов LMS (верхняя строка), RLS (средняя строка) и RLS с экспоненциальным забыванием (нижняя строка): слева -- зависимость ошибки от времени, справа -- итоговые импульсные характеристики фильтров.

Рис.4. Увеличенный график ошибки от времени до изменения характеристик системы.

Однако RLS-алгоритм с экспоненциальным забыванием (который имеет возможность отслеживать изменения статистических свойств обрабатываемого сигнала) отлично справился с изменением характеристик системы, и даже с меньшим переходным процессом, чем LMS. Это подтверждают итоговые импульсные характеристики фильтров (рис.3, справа) -- характеристики LMS-фильтра и RLS с экспоненциальным забыванием хорошо совпадают с экспоненциально затухающими колебаниями импульсной характеристики анализируемой системы (рис.6), а характеристика RLS-фильтра выглядит случайным набором значений.

Если увеличить графики ошибки фильтрации после изменения характеристик системы (рис.6), то можно видеть, что алгоритм RLS с экспоненциальным забыванием, как и прежде, дает меньшую ошибку, чем LMS-алгоритм. Алгоритм RLS на данном этапе мы не рассматриваем, поскольку его ошибку видно на рис.3, и она не сравнима с ошибками других алгоритмов.

Рис.5. Увеличенный график ошибки фильтрации после скачкообразного изменения характеристик системы

Рис.6. Импульсная характеристика системы

Проделанная работа показала, что алгоритм LMS неплохо справляется как со стационарным сигналом, так и с сигналом с изменяющимися характеристиками. Алгоритм RLS в «чистом» виде дает отличные результаты по стационарному сигналу, однако с изменением характеристик системы справиться не способен. Алгоритм RLS с экспоненциальным забыванием отлично справляется с обоими вышеперечисленными сигналами и показывает гораздо меньший переходный процесс и меньшую ошибку по сравнению с остальными рассмотренными алгоритмами.

Однако не стоит забывать, что RLS алгоритмы требуют гораздо больше вычислительных мощностей, чем LMS (2,5N2+4N пар операций «умножение-сложение» по сравнению с N + 1 пар операций на каждом шаге для LMS).

Соответственно, можно сделать вывод, что для случаев, не требующих высокой точности, удобнее всего использовать алгоритм LMS. Если же требуется высокая точность и достаточно вычислительных мощностей, лучше подойдут RLS-алгоритмы. Кроме того, нужно помнить, что RLS-алгоритм в «чистом» виде можно применять только в системах со стационарным сигналом. Если же характер сигнала может меняться, необходимо использовать RLS-алгоритм с экспоненциальным забыванием.

Приложения

Приложение 1

Реализация адаптивных фильтров в MATLAB

Для начала формируем входной и выходной сигналы идентифицируемой системы:

x=randn(2000,1); %дискретный белый гауссов шум

b=exp(-t/5).*cos(t*pi/2); % импульсная характеристика системы

% генерируем первую половину выходного сигнала

Filter(b, 1, x(1:1000));

% генерируем вторую половину выходного сигнала

y(1001:2000) = filter(-b, 1, x(1001:2000), state);

Далее создадим объекты адаптивных фильтров с помощью соответствующих конструкторов. Для LMS-алгоритма предварительно рассчитаем значение коэффициента м, выбрав его в два раза меньшим предельного значения, определяемого формулой (15), а для RLS-алгоритма с экспоненциальным забыванием возьмем л=0,6. Для всех остальных параметров используем значение по умолчанию:

%создаем объекты LMS- и RLS-адаптивных фильтров

N = 32; % длина фильтров

mu = 1/N/var(y); %размер шага для LMS

ha_lms = adaptfilt.lms(N,mu); % создание объекта для LMS

ha_rls = adaptfilt.rls(N); % создание объекта для RLS

ha_erls = adaptfilt.rls(N,0.6); % создание объекта для RLS с экспоненциальным забыванием

Теперь реализуем фильтрацию с помощью функции filter. В соответствии с рис.2, а входной сигнал адаптивного фильтра совпадает с входным сигналом исследуемой системы (x), а образцовый сигнал -- это выходной сигнал системы (y):

%реализуем фильтрацию:

Filter(ha_lms, x, y);

Filter(ha_rls, x, y);

Filter(ha_erls, x, y);

Теперь построим импульсную характеристику системы (рис.3) и графики зависимости сигнала ошибки от времени, а также выведем на экран импульсные характеристики фильтров, полученные на момент завершения обработки сигнала (рис.4):

% строим графики:

title("LMS error")

axis()

impz(ha_lms.coefficients)

title("LMS impulse response")

title("RLS error")

axis()

impz(ha_rls.coefficients)

title("RLS impulse response")

axis()

title("eRLS error")

axis()

impz(ha_erls.coefficients)

title("eRLS impulse response")

title("System impulse response")

Библиографический список

1. Сергиенко А.Б. Цифровая обработка сигналов.-- 3-е изд. -- СПб.: БХВ-Петербург, 2011. -- 768 с.

2. Хайкин С. Нейронные сети. Полный курс. -- М.: Вильямс, 2006. -- 1104 с.

3. Адаптивные фильтры: Пер. С англ./[Под ред. К.Ф.Н.Коуэнам и П.М.Гранта]. -- М.: Мир, 1988. -- 392 с.

4. Шахтарин Б.И. Случайные процессы в радиотехнике: Цикл лекций. -- М.: Радио и связь, 2000. -- 584 с.

5. Основы цифровой обработки сигналов: Курс лекций / А.И.Солонина, Д.А.Улахович, С.М.Арбузов, Е.Б.Соловьева. -- 2-е изд., перераб. и испр. -- Спб.: БХВ-Петербург, 2005. -- 768 с.

6. Солонина А.И. Цифровая обработка сигналов. Моделирование в MATLAB / А.И.Солонина, С.М. Арбузов. -- Спб.: БХВ-Петербург, 2008. -- 816 с.

7. Джиган В. Адаптивные фильтры. Современные средства моделирования и примеры реализации. / В.Джиган // Электроника: наука технология бизнес. -- 2012. -- №7(00121). -- С. 106-125.

8. Сеогиенко А.Б. Алгоритмы адаптивной фильтрации: особенности реализации в MATLAB. / А.Б.Сергиенко // Exponenta Pro. -- 2003. -- №1(1). -- С. 18-28.

Размещено на Allbest.ru

Подобные документы

Классификация адаптивных систем. Достоинства и недостатки типов и классов адаптивных, самонастраивающихся систем. Разработка оригинальной схемы адаптивной системы. Системы со стабилизацией основного контура, идентификатором или уточняемой моделью объекта.

статья , добавлен 24.07.2013


Жесткий и гибкий пороги фильтрации речевого сигнала. Графики вейвлет-разложения речевого сигнала. Блок схема алгоритма фильтрации с гибким порогом. Статистический метод фильтрации речевого сигнала. Оценка качества восстановленного речевого сигнала.

реферат , добавлен 01.12.2008


Исследование цифровой обработки сигналов и её применения в различных сферах деятельности. Изучение достоинств и недостатков медианной фильтрации. Анализ принципов работы медианных фильтров. Реализация медианной фильтрации при помощи MatLab712 R2011a.

курсовая работа , добавлен 04.07.2013


Исходные данные для расчета пассивных RC-фильтров. Расчет параметров элемента фильтра. Частотные фильтры электрических сигналов предназначены для повышения помехоустойчивости различных электронных устройств и систем. Параметры реальных фильтров.

контрольная работа , добавлен 04.10.2008


Характеристика основных требований к методам и алгоритмам фильтрации. Предпосылки возникновения помех и искажений. Особенности фильтров на основе ортогональных и дискретного косинусного преобразований. Применение фильтра со сменным размером окна.

курсовая работа , добавлен 08.12.2011


Общие амплитудно-частотные характеристики (АЧХ) различных типов фильтров. Построение схемы фильтра верхних и нижних частот: активные и пассивные фильтры первого и второго порядка. Принципы действия, функции и применение полосовых и режекторных фильтров.

реферат , добавлен 18.12.2011


Классификация систем радиочастотной идентификации (РЧИ) и области их применения. Состав системы РЧИ, физические принципы работы. Преимущества и недостатки радиочастотной идентификации. Характеристики систем РЧИ и её элементов, международные стандарты.

реферат , добавлен 15.12.2010


Классификация акустических локационных систем по назначению и типу первичного преобразователя, по характеру частотного спектра сигнала, по типу модулирующего воздействия, по избирательности. Область применения датчиков локации. Алгоритм идентификации.

курсовая работа , добавлен 11.08.2010


Конструкция электрических фильтров, технология их изготовления, принцип действия. Меры передачи и параметры фильтров. Использование их в системах многоканальной связи, радиоустройствах, устройствах автоматики, телемеханики. Фильтры нижних частот.

контрольная работа , добавлен 07.04.2016


Сущность и принцип работы сглаживающих фильтров. Классификация и виды. Величины, которые характеризуют качество фильтра. Расчет коэффициента сглаживания. Проектирование активных и пассивных сглаживающих фильтров: достоинства, недостатки, применение.

Основным свойством адаптивной системы является изменяющееся во времени функционирование с саморегулированием. Необходимость такого функционирования очевидна из следующих рассуждений. Если разработчик проектирует «неизменяемую» систему, которую он считает оптимальной, то это означает, что разработчик предвидит все возможные условия на ее входе, по меньшей мере в статистическом смысле, и рассчитывает, что система будет работать при каждом из этих условий. Далее разработчик выбирает критерий, по которому должно оцениваться функционирование, например среднее число ошибок между выходным сигналом реальной системы и выходным сигналом некоторой выбранной модели или «идеальной» системы. Наконец, разработчик выбирает систему, которая оказывается лучшей в соответствии с установленным критерием функционирования, обычно из некоторого априорно ограниченного класса (например, из класса линейных систем).

Однако во многих случаях весь диапазон входных условий может быть не известен точно даже в статистическом смысле или условия могут время от времени изменяться. Тогда адаптивная система, которая, используя регулярный процесс поиска, постоянно ищет оптимум в пределах допустимого класса возможностей, имеет преимущества по сравнению с неизменяемой системой.

Адаптивные системы по своей природе должны быть изменяющимися во времени и нелинейными. Их свойства зависят, помимо всего прочего, от входных сигналов. Если на вход подается сигнал х 1 , то адаптивная система настроится на него и сформирует выходной сигнал – назовем его у 1 . Если на вход подается другой сигнал х 2 , то система настроится на этот сигнал и сформирует выходной сигнал – назовем его у 2 . В общем случае структура и процессы коррекции адаптивной системы будут различными для двух различных входных сигналов.

Для получения оптимального решения существует много методов подстройки значений весовых коэффициентов фильтра. Применялись методы случайных возмущений, которые изменяли весовые коэффициенты фильтра; далее анализировался входной сигнал для того, чтобы установить, приближается ли его случайное возмущение к искомому решению или отдаляет от него. В настоящее время для расчета весовых коэффициентов адаптивных фильтров широко применяется адаптивный алгоритм, основанный на методе наименьших квадратов (МНК), поскольку в нем используются градиентные методы, которые значительно эффективнее, чем другие, обеспечивают сходимость к оптимальному решению. Можно показать, что градиентный метод наименьших квадратов очень схож с методом максимизации отношения сигнал – шум, который разрабатывался с целью применения в тех случаях, когда необходимо получить оптимальные весовые коэффициенты адаптивных антенных решеток. Было также показано, что обнуляющий корректирующий фильтр Лаки является упрощением более общего градиентного метода наименьших квадратов.

Таким образом, адаптивный фильтр - это фильтр, передаточная функция (или частотная характеристика) которого адаптируется, т.е. изменяется таким образом, чтобы пропустить без искажений полезные составляющие сигнала и ослабить нежелательные сигналы или помехи . Схема адаптивного фильтра представлена на рис.5.5.

Рис.5.5. Адаптивный фильтр

Подобный фильтр действует по принципу оценки статистических параметров сигнала и подстройки собственной передаточной функции таким образом, чтобы минимизировать некоторую целевую функцию. Эту функцию обычно формируют с помощью “эталонного” сигнала на задающем входе. Этот эталонный сигнал можно рассматривать как желаемый сигнал на выходе фильтра. Задача блока адаптации состоит в подстройке коэффициентов цифрового фильтра таким образом, чтобы свести к минимуму разность n = n - n , определяеющую ошибку в работе фильтра.

Важнейшей функцией, выполняемой адаптивным фильтром, является моделирование системы. Это иллюстрируется на рис. 5.6, где первичный сигнал с равномерной спектральной плотностью подается непосредственно либо на вход s , либо на вход y адаптивного фильтра. Первичный сигнал поступает на вход системы с импульсной характеристикой H(n) , выход системы соединен со вторым входом адаптивного фильтра. Для получения оптимальных весовых векторов H opt адаптивного фильтра можно применить два разных подхода, которые приведут к совершенно различным результатам. Это имеет место в следующих случаях:

1. Неизвестная система H(n) подключена ко входу y адаптивного фильтра (рис. 5.6, а ). В этом случае оптимальная импульсная характеристика адаптивного фильтра является точной моделью соответствующей характеристики системы H(n) .

2. Неизвестная система H(n) подключена ко входу s адаптивного фильтра (рис. 5.6, б). В этом случае оптимальная импульсная характеристика адаптивного фильтра является функцией, обратной соответствующей характеристике неизвестной системы.

Рис. 5.6. Применение адаптивного фильтра для прямого моделирования системы: H opt =H(n) (а ) и обратного моделирования системы: H opt =H -1 (n) (б ).

Практическим примером, иллюстрирующим работу адаптивного фильтра первого типа (т.е. прямое моделирование системы), является подавление отраженного сигнала в гибридной телефонной линии.

Примером, которым можно воспользоваться для иллюстрации принципа действия адаптивного фильтра, моделирующего обратную характеристику системы, является коррекция искажений при передаче данных по телефонным линиям. В этом случае вход телефонной линии возбуждается известным сигналом, а искаженный сигнал с выхода линии поступает на вход s(n) адаптивного фильтра. Затем фильтр перестраивается с помощью подачи на вход y(n) последовательной серии известных (неискаженных) первичных сигналов. Адаптивный фильтр моделирует импульсную характеристику, обратную характеристике линии, для получения на выходе отфильтрованных (свободных от искажений) данных.

Следующая область применения адаптивных фильтров – подавление шумов. В этой схеме первичный сигнал, содержащий искомую информацию наряду с мешающим сигналом, приложен к входу y(n) . Затем от другого источника, не содержащего никаких составляющих исходного сигнала, поступает независимый коррелированный сигнал – образец мешающего сигнала. Если этот коррелированный сигнал поступает непосредственно на вход s(n) адаптивного фильтра, фильтр формирует импульсную характеристику, обеспечивающую получение выходного сигнала y(n) , который когерентно вычитает из y(n) нежелательную составляющую, оставляя на выходе e(n) лишь искомый сигнал.

Одним из примеров использования этого метода является регистрация сердцебиения плода. Первичный сигнал поступает от преобразователь, расположенного на поверхности живота матери. Этот преобразователь вырабатывает сигнал, содержащий импульсы сердечных сокращений плода, которые, однако, существенно маскируются сердцебиением матери. Затем от второго преобразователя, расположенного на груди матери, получают вторичный сигнал, регистрирующий только сердцебиение матери. Далее адаптивный фильтр моделирует тракт искажений от преобразователя, расположенного на груди, до преобразователя расположенного на животе, для получения сигнала, который когерентно вычитается из сигнала с поверхности живота. Адаптивные фильтры применяются и в других случаях, например для устранения шума двигателя в микрофоне пилота в кабине самолета или для подавления акустических шумов окружающей среды, например на крупных электростанциях.

Еще одно применение адаптивных фильтров – это реализация самонастраивающегося фильтра, используемого для выделения синусоиды, маскируемой широкополосным шумом. Такое применение в адаптивном линейном усилителе (АЛУ) осуществляется путем подачи сигнала непосредственно на вход фильтра y(n) и подачи модификации сигнала с временной задержкой на вход фильтра s(n) . В том случае, если задержка превышает величину, обратную ширине полосы пропускания фильтра, шумовые составляющие на двух входах не будут коррелированы. Адаптивный фильтр дает на выходе синусоиду с увеличенным отношением сигнал – шум, тогда как на выходе сигнала ошибки синусоидальные составляющие уменьшаются.

Адаптивные фильтры БИХ-типа в основном применялись для решения таких проблем, как ослабление влияния многолучевого распространения сигнала в системах радиолокации и радиосвязи. В этом случае принятый сигнал содержит исходный передаваемый сигнал, свернутый с импульсной характеристикой канала, которая при многолучевом распространении содержит только нули. Тогда для исключения интерференционных помех адаптивный приемник моделирует характеристику, обратную характеристике канала (рис. 5.6, б ). Это наиболее эффективно осуществляется путем использования модели адаптивного фильтра с характеристикой, содержащей только полюсы, причем положения полюсов подбираются таким образом, чтобы они совпадали с положением нулей в характеристике канала.

При конструировании адаптивного фильтра КИХ-типа можно также учесть эту модель, но более экономично воспользоваться рекурсивной структурой, поскольку она реализует инверсную структуру фильтра при его более низком порядке и с меньшими весами. Отсюда с полным основанием можно сказать, что такая структура будет обеспечивать более быструю сходимость, чем ее трансверсальный аналог. Однако для обеспечения устойчивости адаптивного рекурсивного фильтра необходима высокая степень точности при расчете цифровой схемы. Метод адаптивной обработки сигналов на основе фильтров БИХ-типа применяется в электронных радиолокационных измерительных приемниках для выделения импульсов. Адаптивные же фильтры Калмана представляют интерес для идентификации типов радиолокационных колебаний, генерируемых определенными типами излучателей. Они также находят применение при фильтрации и ослаблении влияния многолучевого распространения в высокочастотных (от 3 до 30 МГц) каналах цифровой связи, где первоочередное значение имеет присущая этим фильтрам высокая скорость сходимости.

Большинство КИХ-фильтров строится с учетом довольно простых общепринятых допущений. Эти допущения приводят к хорошо известным несложным алгоритмам адаптации (например, МНК), реализация которых подробно разработана в отношении скорости сходимости, остаточной ошибки и т.д. Таким подходом шире всего пользуются при применении адаптивных фильтров в системах дальней связи, например для выравнивания и гашения отраженного сигнала.

В 1971 г. Чанг внес значительный вклад в классификацию типов фильтров: он предпринял попытку объединить все подходы и создать одну обобщенную структуру выравнивателя, или корректирующего фильтра (рис. 5.7.). Эта структура содержит набор произвольных фильтров, подключенных к линейной взвешивающей и комбинирующей цепи. Фильтр КИХ-типа можно получить из этой обобщенной структуры путем замены произвольного фильтра линией задержки с отводами, дающей на выходах серию выборок сигнала с временной задержкой. Фильтр БИХ-типа благодаря наличию элементов рекурсивной обратной связи осуществляет дальнейшую обработку сигнала до получения выборок сигнала с временной задержкой, которые последовательно поступают на взвешивающую и комбинирующую цепь.

Вечная проблема любых измерений их низкая точность. Основных способов повышения точности два, первый состоит в повышении чувствительности к измеряемой величине, однако при этом как правило растёт чувствительность и к неинформативным параметрам, что требует принятия дополнительных мер по их компенсации. Второй способ состоит в статистической обработке многократных измерений, при этом среднеквадратичное отклонение обратно пропорциональна корню квадратному из числа измерений.

Статистические методы повышения точности разнообразны и многочисленны, но и они делятся на пассивные для статических измерений и активные для динамических измерений, когда измеримая величина изменяется во времени. При этом сама измеряемая величина так же, как и помеха являются случайными величинами с изменяющимися дисперсиями.

Адаптивность методов повышения точности динамических измерений следует понимать, как использование прогнозирования значений дисперсий и погрешности для следующего цикла измерений. Такое прогнозирование осуществляется в каждом цикле измерений. Для этой цели применяются фильтры Винера, работающие в частотной области. В отличии от фильтра Винера, фильтр Калмана работает во временной, а не в частотной области. Фильтр Калмана был разработан для многомерных задач, формулировка которых осуществляется в матричной форме. Матричная форма достаточно подробно описана для реализации на Python в статье , . Описание работы фильтра Калмана, приведенная в указанных статьях, рассчитана на специалистов в области цифровой фильтрации. Поэтому возникла необходимость рассмотреть работу фильтра Калмана в более простой скалярной форме.

Немного теории

Рассмотрим схему фильтра Калмана для его дискретной формы.

Здесь G(t) блок работа которого описывается линейными соотношениями. На выходе блока вырабатывается неслучайный сигнал y(t). Этот сигнал суммируется с шумом w(t), который возникает внутри контролируемого объекта. В результате такого сложения, получаем новый сигнал x(t). Этот сигнал представляет сумму неслучайного сигнала и шума и является случайный сигналом. Далее сигнал x(t) преобразуется линейным блоком H(t), суммируясь с шумом v(t), распределённым по-другому чем w(t) закону. На выходе линейного блока H(t) получаем случайный сигнал z(t), по которому и определяется неслучайный сигнал y(t). Следует отметить, что линейные функции блоков G(t) и Н(t) могут тоже зависеть от времени.

Будем считать, что случайные шумы w(t) и v(t) - это случайные процессы с дисперсиями Q, R и нулевыми математическими ожиданиями. Сигнал x(t) после линейного преобразования в блоке G(t) распределён во времени по нормальному закону. С учётом изложенного, соотношение для измеряемого сигнала примет вид:

Постановка задачи

После фильтра нужно получить максимально возможное приближение y"" к неслучайному сигналу y(t).

При непрерывном динамическом измерении каждое следующее состояние объекта, а, следовательно, и значение контролируемой величины отличается от предыдущего по экспоненциальному закону с постоянного времени T в текущем интервале времени ,

Ниже приведена программа на Python, в которой решается уравнение для неизвестного не зашумлённого сигнала y(t). Процесс измерения рассматривается для суммы двух псевдослучайных величин, каждая из которых образуется как функция нормального распределения от равномерного распределения.

Программа для демонстрации работы дискретного адаптивного фильтра Калмана

#!/usr/bin/env python #coding=utf8 import matplotlib.pyplot as plt import numpy as np from numpy import exp,sqrt from scipy.stats import norm Q=0.8;R=0.2;y=0;x=0#начальные дисперсии шумов(выбраны произвольно) и нулевые значения переменных. P=Q*R/(Q+R)# первая оценка дисперсий шумов. T=5.0#постоянная времени. n=;X=;Y=;Z=#списки для переменных. for i in np.arange(0,100,0.2): n.append(i)#переменная времени. x=1-exp(-1/T)+x*exp(-1/T)#модельная функция для x. y=1-exp(-1/T)+y*exp(-1/T)# модельная функция для y. Y.append(y)#накопление списка значений y. X.append(x)# накопление списка значений x. norm1 = norm(y, sqrt(Q))# нормальное распределение с #математическим ожиданием – y. norm2 = norm(0, sqrt(R))#))# нормальное распределение с #математическим ожиданием – 0. ravn1=np.random.uniform(0,2*sqrt(Q))#равномерное распределение #для шума с дисперсией Q. ravn2=np.random.uniform(0,2*sqrt(R))# равномерное распределение #для шума с дисперсией R. z=norm1.pdf(ravn1)+norm2.pdf(ravn2)#измеряемая переменная z. Z.append(z)# накопление списка значений z. P=P-(P**2)/(P+Q+R) #переход в новое состояние для x. x=(P*z+x*R)/(P+R)# новое состояние x. P=(P*R)/(P+R)# прогноз для нового состояния x. plt.plot(n, Y, color="g",linewidth=4, label="Y") plt.plot(n, X, color="r",linewidth=4, label="X") plt.plot(n, Z, color="b", linewidth=1, label="Z") plt.legend(loc="best") plt.grid(True) plt.show()

В чём отличие предложенного алгоритма от известного

Мною был усовершенствован алгоритм работы фильтра Калмана, приведенный в методических указаниях для Mathcad:

В результате преждевременной смены состояния для сравниваемой переменной x(t) возрастала погрешность на участке резких изменений:

Тогда как в моем алгоритме используется начальная прогнозная оценка влияния шумов. Это дало возможность уменьшить ошибку измерений v(t).

В приведённом алгоритме используются заданные – модельные экспоненциальные функции, поэтому для наглядности приведём их отдельно на общем графике работы фильтра Калмана.

Код программы для графического анализ работы фильтра

#!/usr/bin/env python #coding=utf8 import matplotlib.pyplot as plt import numpy as np from numpy import exp,sqrt from scipy.stats import norm Q=0.8;R=0.2;y=0;x=0#начальные дисперсии шумов(выбраны произвольно) и нулевые значения переменных. P=Q*R/(Q+R)# первая оценка дисперсий шумов. T=5.0#постоянная времени. n=;X=;Y=;Z=#списки для переменных. for i in np.arange(0,100,0.2): n.append(i)#переменная времени. x=1-exp(-1/T)+x*exp(-1/T)#модельная функция для x. y=1-exp(-1/T)+y*exp(-1/T)# модельная функция для y. Y.append(y)#накопление списка значений y. X.append(x)# накопление списка значений x. norm1 = norm(y, sqrt(Q))# нормальное распределение с #математическим ожиданием – y. norm2 = norm(0, sqrt(R))#))# нормальное распределение с #математическим ожиданием – 0. ravn1=np.random.uniform(0,2*sqrt(Q))#равномерное распределение #для шума с дисперсией Q. ravn2=np.random.uniform(0,2*sqrt(R))# равномерное распределение #для шума с дисперсией R. z=norm1.pdf(ravn1)+norm2.pdf(ravn2)#измеряемая переменная z. Z.append(z)# накопление списка значений z. P=P-(P**2)/(P+Q+R) #переход в новое состояние для x. x=(P*z+x*R)/(P+R)# новое состояние x. P=(P*R)/(P+R)# прогноз для нового состояния x. plt.subplot(221) plt.plot(n, Y, color="g",linewidth=2, label="Модельная функция \n не зашумленой \n переменной") plt.legend(loc="best") plt.grid(True) plt.subplot(222) plt.plot(n, X, color="r",linewidth=2, label="Модельная функция \n сравниваемой \n переменной") plt.legend(loc="best") plt.grid(True) plt.subplot(223) plt.plot(n, Z, color="b", linewidth=1, label="Измеряемая функция \n псевдослучайных переменных") plt.legend(loc="best") plt.grid(True) plt.subplot(224) plt.plot(n, Y, color="g",linewidth=2, label="Y") plt.plot(n, X, color="r",linewidth=2, label="X") plt.plot(n, Z, color="b", linewidth=1, label="Z") plt.legend(loc="best") plt.grid(True) plt.show()

Результат работы программы

Выводы

В статье описана модель простой скалярной реализации фильтра Калмана средствами условно бесплатного языка программирования общего назначения Python, что расширит область её применения в целях обучения.