Тема 12. Оптимальные линейные цифровые фильтры ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ Тема 12. ОПТИМАЛЬНЫЕ ЛИНЕЙНЫЕ ЦИФРОВЫЕ ФИЛЬТРЫ. Как много дел считалось невозможными, пока они не были осуществлены. Гай Плиний Секунд (философ). Специалисты в науке подобны старателям. Стоит одному найти крупинку золота, как другие выроют в этом месте котлован. А тема оптимальности, это вообще золотое Эльдорадо, можно копать в любом направлении. Владимир Старцев. Уральский геофизик, ХХ в. Содержание Введение. 1. Случайные процессы и шумы. Белый шум. Модель белого шума. Фильтрация белого шума. 2. Критерии построения оптимальных фильтров. Среднее квадратическое отклонение. Амплитудное отношение сигнал/шум. Энергетическое отношение сигнал/шум. 3. Фильтр Колмогорова-Винера. Условие оптимальности фильтра. Система линейных уравнений фильтра. Частотная характеристика фильтра. Задание мощности шумов. Эффективность фильтра. Пример расчета оптимального фильтра воспроизведения сигнала. Фильтры прогнозирования и запаздывания. 4. Оптимальные фильтры сжатия сигналов. Условие оптимальности. Частотная характеристика. Примеры использования. 5. Фильтр обнаружения сигналов. Частотная характеристика. Система линейных уравнений. Эффективность фильтра. Согласованный фильтр. Обратный фильтр. 6. Энергетический фильтр. Критерий оптимальности. Расчет векторов операторов фильтров. ВВЕДЕНИЕ Результаты практических измерений, подлежащие обработке, содержат определенный полезный сигнал на фоне различного рода помех (шумов), при этом спектр помех в общем случае представлен по всему интервалу главного частотного диапазона и наложен на спектр полезного сигнала. В этих условиях ставится задача реализации оптимальных фильтров, которые позволяют достаточно надежно производить обнаружение сигнала, наилучшим образом выделять сигнал на фоне помех или подавлять помехи без существенного искажения сигнала. Главным критерием при проектировании оптимальных фильтров, как правило, является минимизация среднеквадратичной ошибки восстановления полезного сигнала. Линейные оптимальные фильтры, которые рассматриваются в настоящей теме, обычно базируются на оптимальном фильтре Колмогорова-Винера. 12.1. случайные процессы и шумы /12/. Случайные процессы и шумы описываются функциями автокорреляции и спектрами мощности. Модели случайных процессов и сигналов с заданными статистическими характеристиками обычно получают фильтрацией белого шума. Белый шум является стационарным случайным процессом q(t), у которого автокорреляционная функция описывается дельта - функцией Дирака, спектральная плотность мощности не зависит от частоты и имеет постоянное значение Wq (f) = s2 , равное дисперсии значений q(t). Все спектральные составляющие белого шума имеют одинаковую мощность. По существу, это идеализированный случайный процесс с бесконечной энергией. Но в случае постоянства спектральной плотности мощности случайного процесса в конечном диапазоне частот такая идеализация позволяет достаточно просто разрабатывать оптимальные методы фильтрации. Многие помехи в радиотехнике, в технике связи и в других отраслях, в том числе в информатике, рассматривают как белый шум, если эффективная ширина спектра сигналов Bs много меньше эффективной ширины спектра шумов Bq , а спектральная плотность мощности шумов слабо изменяется в интервале спектра сигнала. Понятие "белый шум" определяет только спектральную характеристику случайного процесса и под это понятие подпадают любые случайные процессы, имеющие равномерный энергетический спектр и различные законы распределения. Если частотный диапазон спектра, на котором рассматриваются сигналы и помехи, равен 0-В, то спектральная плотность шума задается в виде: Wq (f)=s2 , 0£ f £B; Wq (f)=0, f > B, (12.1.1) при этом корреляционная функция шума определяется выражением: Rq (t)= s2 B×sin(2pBt)/2pBt. (12.1.2) Эффективный интервал корреляции: Tk = 2 |Rq (t)|dt /Rq (0). (12.1.3) Рис. 12.1.1. Функции корреляции белого шума в частотном интервале 0-В. Реальный интервал корреляции целесообразно определять по ширине главного максимума функции Rq (t) (значения t при первых пересечениях нулевой линии), в котором сосредоточена основная часть энергии шумов, при этом Tk = 1/В и BTk = 1. Как следует из всех этих выражений и наглядно видно на рис. 12.1.1, при ограничении частотного диапазона в шумах появляется определенная корреляция между значениями. Чем меньше частотный диапазон шумов, тем больше их радиус корреляции. Ограничение шумов определенным частотным диапазоном эквивалентно фильтрации белого шума частотным фильтром с соответствующей шириной полосы пропускания, при этом корреляционная функция импульсного отклика фильтра свертывается с дельта – функцией белого шума. Модель белого шума q(t) можно формировать как случайную по времени (аргументу) последовательность дельта - импульсов d(ti ) со случайными амплитудными значениями ai : q(t) = Si ai d(t-ti ), (12.1.4) которая удовлетворяет условиям статистической однородности: постоянное среднее число импульсов в единицу времени и статистическая независимость появления каждого импульса от предыдущих. Такой поток импульсов называют пуассоновским, он является некоррелированным и имеет равномерный спектр плотности мощности: Wq (w) = c2 = Nsa 2 , где N - число импульсов на интервале Т реализации случайного процесса, sa 2 -дисперсия амплитуд импульсов. Фильтрация белого шума. Если на входе фильтра с импульсным откликом h(t) действует белый шум q(t), то сигнал на выходе фильтра: g(t) = h(t) ③ q(t) = h(t) ③ Si ai d(t-ti ) = Si ai h(t-ti ), (12.1.5) т.е. выходной сигнал будет представлять собой последовательность сигналов импульсной реакции фильтра h(t) с амплитудой ai , при этом автокорреляционная функция и спектр мощности выходного потока также становятся подобными ФАК и спектру мощности импульсной реакции фильтра, и в первом приближении определяются выражениями: Rg (t) N sa 2 Rh (t) = c2 Rh (t), (12.1.6) Wg (w) N sa 2 |H(w)|2 = c2 |H(w)|2 . (12.1.7) Этот результат известен как теорема Кэмпбелла. 12.2. Критерии построения оптимальных фильтров. В практике обработки данных используются три основных критерия построения оптимальных фильтров: минимум среднего квадратического отклонения профильтрованного сигнала от его действительного или заданного значения, максимум отношения сигнал/шум и максимум энергетического отношения сигнал/шум на выходе фильтра. Критерии исходят из вероятностно - статистической модели обрабатываемых данных. При анализе и синтезе фильтров используется аддитивная модель входного сигнала: x(k) = s(k)+q(k), где s(k) - полезная составляющая сигнала, q(k) - составляющая шумов и помех. Синтез оптимальных фильтров производится с максимальным использованием известной информации как о сигналах, которые необходимо выделить, так и о шумах и помехах. Как правило, используется информация о природе полезного сигнала и шума, об их спектральном составе, о корреляционных и взаимных корреляционных характеристиках. Наличие определенных особенностей (различий) в характеристиках сигнала и шума позволяет реализовать фильтр вообще и оптимальный фильтр в частности. В геофизической практике априорные данные о полезных сигналах являются достаточно определенными, особенно для активных методов геофизики (сейсмические методы, электроразведка на переменном токе, индукционные методы ядерной геофизики и пр.). Определение характеристик действующих помех представляет собой более сложную проблему, но даже при полной неопределенности можно допустить, что помеха является нормальным стационарным процессом с нулевым средним значением. Среднее квадратическое отклонение. При наличии помех точное выделение полезного сигнала методами линейной фильтрации, как правило, невозможно. Результат фильтрации y(k) = h(n) ③ x(k-n) (12.2.1) отличается от s(k) на величины e(k) = y(k)-s(k), которые являются абсолютными значениями погрешности воспроизведения полезного сигнала по координатам k. Качество фильтра оценивается средним значением квадрата величины e(k): . (12.2.2) Во многих задачах анализа данных не требуется восстановления исходной формы сигнала s(k), т.к. в процессе его дальнейшей обработки осуществляется преобразование сигнала s(k) в сигнал z(k), форма которого может быть более удобной для извлечения (измерения) каких-либо информационных параметров сигнала (амплитуды, частоты, длительности и т.п.). В этом случае фильтр может проектироваться непосредственно на получение выходного сигнала z(k). Качество таких формирующих фильтров оценивается средним значением квадрата величины e(k) получения сигнала заданной формы: . (12.2.2') Выражения (12.2.2) дают возможность определить значения h(k) фильтра по критерию минимума среднего квадратического отклонения выходного сигнала от его действительной или заданной формы. Амплитудное отношение сигнал/шум. При постановке задачи обнаружения (установления факта наличия) в экспериментальных данных сигнала известной формы для проектирования фильтра используется критерий максимума пикового отношения сигнал/шум на выходе фильтра: rа = yэкс /s, где yэкс - экстремальное значение амплитуды сигнала, s - средний квадратический уровень значений помех. Если в полезном сигнале отсутствует четко выраженный экстремум, а сам сигнал достаточно протяженный по аргументу, то в качестве критерия используется отношение средних квадратов амплитуд сигнала и шума: , (12.2.3) где y2 - средний квадрат амплитуды сигнала в пределах его формы. Энергетическое отношение сигнал/шум. При узко конкретной задаче обнаружения сигнала степень искажения самого сигнала может не ограничиваться. Если кроме обнаружения сигнала, как основной цели обработки данных, ставится и задача оценки его формы, то в этом случае для проектирования фильтра обычно используется критерий максимума энергетического отношения сигнал/шум: r = Еy /Eq , (12.2.4) где Еу и Eh - энергия соответственно сигнала и шума на выходе фильтра. 12.3. Фильтр Колмогорова-Винера. Условие оптимальности фильтра. Фильтр Колмогорова-Винера является оптимальным фильтром формирования из входного сигнала x(t) выходного сигнала z(t) при известной форме полезного сигнала s(t), который содержится во входном сигнале в сумме с шумами. В качестве критерия его оптимизации используется среднее квадратическое отклонение сигнала y(t) на выходе фильтра от заданной формы сигнала z(t). Подставим уравнение свертки (12.2.1) в раскрытой форме весового суммирования в выражение (12.2.2') и получим отклонение e2 выходного сигнала y(k) = h(n) ③ x(k-n) от заданной формы выходного сигнала z(k) по всем точкам массива данных: . (12.3.1) В частном случае воспроизведения формы полезного сигнала в качестве функции z(k) принимается функция s(k). Минимум выражения (12.3.1) определяет значения коэффициентов h(n) оптимального фильтра. Для нахождения их значений продифференцируем выражение (12.3.1) по коэффициентам фильтра и приравняем полученные уравнения нулю. В итоге получаем: , где - корреляционная функция входного сигнала, - взаимная корреляционная функция сигналов z(k) и x(k). Отсюда: hn R(m-n) = B(m), n = m = 0,1,2, ... , M. (12.3.2) В краткой форме символической записи: h(n) ③ R(m-n) = B(m). (12.3.3) Другими словами, свертка функции отклика оптимального фильтра с функцией автокорреляции входного сигнала должна быть равна функции взаимной корреляции выходного и входного сигналов. Система линейных уравнений фильтра. Выражение (12.3.2) может быть записано в виде системы линейных уравнений - однострочных равенств левой и правой части для фиксированных значений координаты m коэффициентов фильтра. При расчете симметричных фильтров, не сдвигающих фазы выходного сигнала, фильтр может вычисляться только одной половиной своих значений: m=0: ho R(0)+ h1 R(1)+ h2 R(2)+ h3 R(3)+ ...+ hM R(M) = B(0), (12.3.3') m=1: ho R(1)+ h1 R(0)+ h2 R(1)+ h3 R(2)+ ...+ hM R(M-1) = B(1), m=2: ho R(2)+ h1 R(1)+ h2 R(0)+ h3 R(1)+ ...+ hM R(M-2) = B(2), .................................................................................................... m=M: ho R(M)+ h1 R(M-1)+ h2 R(M-2)+ .... + hM R(0) = B(M). Решение данной системы уравнений относительно значений hi дает искомый оператор фильтра. При расчете каузальных (односторонних) фильтров выходной сигнал z(t) следует задавать со сдвигом вправо по оси координат таким образом, чтобы значимая часть функции взаимной корреляции B(m) полностью располагалась в правой части системы уравнений (12.3.3'). Отметим, что R(m) = Rs (m)+Rq (m), где Rs - функция автокорреляции сигнала, Rq - функция автокорреляции шума, а B(m) = Bzs (m)+Bzq (m), где Bzs - функция взаимной корреляции сигналов z(k) и s(k), Bzq - функция взаимной корреляции сигнала z(k) и помех q(k). Подставляя данные выражения в (12.3.3), получаем: h(n) ③ [Rs (m-n)+Rq (m-n)] = Bzs (m)+Bzq (m). (12.3.4) Частотная характеристика фильтра находится преобразованием Фурье левой и правой части уравнения (12.3.4): H(w)[Ws (w)+Wq (w)] = Wzs (w)+Wzq (w), H(w) = [Wzs (w)+Wzq (w)] / [Ws (w)+Wq (w)], (12.3.5) где Ws (w) - Rs (m) и Wq (w) - Rq (m) - энергетические спектры (плотности мощности) сигнала и помех, Wzs (w) - Bzs (m) - взаимный энергетический спектр входного и выходного сигналов, Wzq (w) - Bzq (m) - взаимный энергетический спектр выходного сигнала и помех. Обычно имеет место статистическая независимость полезного сигнала, а, следовательно, и сигнала z(k), от шумов, при этом Bzq = 0 и фильтр называют оптимальным по сглаживанию шумов при заданной форме выходного сигнала: H(w) = Wzs (w) / [Ws (w)+Wq (w)], (12.3.6) Фильтр (12.3.6) оптимален в том смысле, что максимизирует отношение мощности сигнала к мощности шума по всему интервалу сигнала, но не в каждой индивидуальной точке. Выражения (12.3.5-12.3.6), как правило, всегда имеют решение. Однако это не означает возможности формирования фильтром любой заданной формы выходного сигнала. Можно предполагать, что если спектр заданного сигнала z(t) шире значимой части спектра полезного сигнала s(t), то формирование требуемых высоких частот заданного сигнала из незначимых частот спектра полезного сигнала может потребовать огромных коэффициентов усиления на этих частотах, под действие которых попадут и частотные составляющие шумов. Результат такой операции непредсказуем. Эти практические соображения можно распространить на все частотные соотношения входного и выходного сигналов линейных фильтров: значимые гармоники спектров выходных сигналов должны формироваться из значимых гармоник спектров входных сигналов. Если заданная форма сигнала z(k) совпадает с формой полезного сигнала s(k), то B(m) = Bss = Rs и фильтр называют фильтром воспроизведения полезного сигнала : H(w) = Ws (w) / [Ws (w)+Wq (w)], (12.3.7) Выражения (12.3.6-12.3.7) достаточно наглядно демонстрируют физический смысл формирования передаточной функции фильтра. При воспроизведении сигнала частотная функция взаимной корреляции входного сигнала с выходным Wzs (плотность взаимной мощности) повторяет частотную функцию автокорреляции Ws (плотность мощности сигнала). Плотность мощности статистических шумов Wq распределена по частотному диапазону равномерно, в отличие от плотности мощности сигнала Ws , которая, в зависимости от формы сигнала, может занимать любые частотные интервалы спектрального диапазона. Частотная передаточная функция фильтра воспроизведения сигнала формируется отношением Ws (w)/[Ws (w)+Wq (w)]. На частотах, где сосредоточена основная энергия сигнала, имеет место Ws (w)>>Wq (w) и H(w) Þ 1 (как минимум, больше 0.5). Там, где значение Ws (w) меньше Wq , коэффициент передачи фильтра становится меньше 0.5. В пределе H(w)=0 на всех частотах, где полностью отсутствуют частотные составляющие сигнала. Аналогичный процесс имеет место и при формировании произвольного сигнала z(t) на выходе фильтра. В этом случае из частот входного сигнала устанавливаются на выделение и усиление частоты взаимной мощности входного и выходного сигнала, необходимые для формирования сигнала z(t), причем коэффициент усиления на этих частотах может быть много больше 1, а подавляться могут не только шумы, но и частоты основного сигнала, если их нет в сигнале z(t). Таким образом, оптимальные фильтры учитывают особенности спектрального состава сигналов и способны формировать передаточные функции любой сложности на выделение полезных частот сигналов из любых диапазонов спектра с максимальных подавлением шумов на всех частотах спектрального диапазона, не содержащих полезных сигналов, при этом границы усиления-подавления устанавливаются автоматически по заданному уровню шумов. Задание мощности шумов. Следует внимательно относиться к заданию функции шумов Wq(w). При полной неопределенности шума необходимо, как минимум, выполнять оценку его дисперсии s2 и распространять на весь частотный диапазон с соответствующей нормировкой (2 Wq(w) dw = s2 ), т.е. считать его белым шумом. При известной функции полезного сигнала в зарегистрированной реализации значение дисперсии шумов в первом приближении может быть оценено по разности дисперсий реализации и функции полезного сигнала. Можно выполнить и выделение шумов из входного сигнала в отдельный шумовой массив, например, вейвлетным преобразованием. Однако использовать выделенный шум непосредственно для вычисления функции Wq(w) допустимо только по достаточно представительной реализации при условии стационарности и эргодичности шума. В противном случае функция Wq(w) будет отображать только распределение шумов в зарегистрированной реализации сигнала, а соответственно фильтр будет оптимален только к этой реализации, что не гарантирует его оптимальности к любой другой реализации. Но для обработки единичной зарегистрированной реализации сигнала такой метод не только вполне допустим, но может существенно повысить точность формирования выходного сигнала. Эффективность фильтра. Из выражений (12.3.5-12.3.7) следует, что с позиции минимального искажения полезного сигнала при максимальном подавлении шумов фильтр Колмогорова-Винера эффективен тем в большей степени, чем больше отношение сигнал/шум на входе фильтра. В пределе, при Wq (w)<<Ws (w) имеем H(w)Þ 1 и фильтр воспроизводит входной (или заданный) сигнал без искажений (помех нет, исправлять нечего). Отметим также, что помеха, коррелированная с полезным сигналом, как это следует из (12.3.5), используется фильтром для повышения точности воспроизведения сигнала. С другой стороны, при Wq (w)>>Ws (w) имеем H(w) Þ 0 и сигнал будет сильно искажен, но никакой другой фильтр не обеспечит лучшего результата. Пример. Расчет оптимального фильтра воспроизведения сигнала. Выполняется в среде Mathcad . Форма входного сигнала считается известной и близка к гауссовой. На входной сигнал наложен статистический шум с равномерным распределением мощности по всему частотному диапазону (белый шум), некоррелированный с сигналом, и функцию Wzq принимаем равной нулю. Для наглядного просмотра связи параметров фильтра с параметрами сигнала задаем модели сигналов в двух вариантах: K := 1000 k := 0 .. K A := 50 s1k := A·exp[-0.0005·(k-500)2 ] s2k := A·exp[-0.00003·(k-500)2 ] Ü информационные сигналы Q := 30 qk := rnd(Q) – Q/2 x1k := s1k + qk x2k := s2k + qk Ü входные сигналы Рис. 12.3.1. Модельные сигналы. В качестве выходных сигналов задаем те же самые функции s1 и s2. Быстрым преобразованием Фурье вычисляем спектры сигналов и формируем спектры плотности мощности: S1 := CFFT(s1) S2 := CFFT(s2) Q := CFFT(q) Ü спектры сигналов Ü спектры мощности Ds1 := var(s1) Ds2 := var(s2) Dx1 := var(x1) Dx2 := var(x2) Dq := var(q) Ü дисперсии Ds1 = 124.308 Ds2 = 310.264 Dx1 = 202.865 Dx2 = 386.78 Dq = 79.038 Ü информация mean(Wq) = 0.079 Wq1 := (Dx1 – Ds1)/(K+1) Wq1 = 0.078 Ü информация Wq2 := (Dx2 – Ds2)/(K+1) Wq2 = 0.076 Ü информация Wqk := Wq1 Ü замена на постоянное распределение Для воспроизведения сигналов вычисления функций Wzs не требуется, т.к. Wzs = Ws. Вычисление Wq также имеет только информативный характер. Передаточные функции оптимальных фильтров (приведены на рис. 12.3.2): Рис. 12.3.2. Передаточные функции оптимальных фильтров в сопоставлении с нормированными модулями спектров сигналов Как следует из рисунка 12.3.2, для плавных монотонных функций, основная энергия которых сосредоточена в низкочастотной части спектра, передаточные функции оптимальных фильтров, по существу, представляют собой низкочастотные сглаживающие фильтры с автоматической подстройкой граничной частоты пропускания под основные частоты входного сигнала. Операторы фильтров можно получить обратным преобразованием Фурье: h1 := ICFFT(H1)/(K+1) h2 := ICFFT(H2)/(K+!) Ü обратное преобразование Фурье Рис. 12.3.3. Импульсные отклики фильтров. Оператор фильтра, в принципе, бесконечен. В данном случае, при использовании БПФ максимальное число отсчетов равно К/2 = 500. Усечение размеров оператора может выполняться по типовым методам с применением весовых функций (в расчете операторы нормируются к 1, весовые функции не применяются). N1 := 160 n1 := 0 .. N1 N2 ;= 500 n2 := 0 .. N2 Ü размеры и нумерация операторов hm1 := h10 + 2· h1n 1 hm1=0.988 h1 := h1/hm1 Ü нормировка hm2 := h20 + 2· h2n 2 hm2=1.001 h2 := h2/hm2 Ü нормировка Ü свертка Рис. 12.3.4. Проверка действия оптимальных фильтров. Коэффициент усиления дисперсии помех Þ Kd := (h0 )2 + 2· hn Kd1=0.021 Kd2= 0.0066 Среднеквадратическое отклонение воспроизведения сигнала: e1= 1.238 e2 = 0.701 Проверка действия оператора фильтра приведена на рис. 12.3.4. Особую эффективность оптимальный фильтр имеет при очистке от шумов сигналов, имеющих достаточно сложный спектральный состав. Оптимальный фильтр учитывает конфигурацию спектра сигнала и обеспечивает максимальное подавление шумов, в том числе внутри интервала основного частотного диапазона сигнала. Это наглядно видно на рис. 12.3.5 для сигнала, близкого к прямоугольному, спектр которого имеет кроме основной низкочастотной части затухающие боковые осцилляции. Расчет фильтра выполнялся по методике, приведенной в примере 1. Рис. 12.3.5. Оптимальная фильтрация сигнала со сложным спектральным составом. Рис. 12.3.6. Оптимальная фильтрация радиоимпульса. На рис. 12.3.6 приведен пример фильтрации оптимальным фильтром радиоимпульса. Основной пик спектра радиоимпульса находится в области несущей частоты, а боковые полосы определяются формой модулирующего сигнала, в данном случае – прямоугольного импульса. На графике модулей сигнала и передаточной функции фильтра можно видеть, что оптимальный фильтр превратился в полосовой фильтр, при этом учитывается форма боковых полос спектра сигнала. Фильтры прогнозирования и запаздывания. Если в правой части уравнения (12.3.3) желаемым сигналом задать входной сигнал со сдвигом на величину kDt, то при этом B(m) = R(m+k) и уравнение принимает вид: h(n) ③ R(m-n) = R(m+k). (12.3.8) При k > 0 фильтр называется фильтром прогнозирования и вычисляет будущие значения сигнала по его предшествующим значениям. При k < 0 фильтр является фильтром запаздывания. Реализация фильтра заключается в решении соответствующих систем линейных уравнений для каждого заданного значения k. Фильтр может использоваться для интерполяции геофизических полей, в том числе в наперед заданные точки, а также для восстановления утраченных данных. 12.4. Оптимальные фильтры сжатия сигналов. Условие оптимальности. Фильтр сжатия сигнала x(t), по существу, представляет собой фильтр формирования сигнала z(t) с эффективной шириной длительности, меньшей по сравнению с эффективной шириной длительности полезного сигнала s(t) во входном сигнале x(t). Расчет оптимального фильтра сжатия может выполняться непосредственно по выражениям (12.3.3). В предельном случае сжатия сигнала до импульса Кронекера положим, что z(k)=d(k) при статистической независимости сигнала и шума. Отсюда: Bsz (m) = d(m) ③ s(k+m) = s(-m). h(n) ③ (Rs(m-n)+Rq(m-n)) = s(-m). (12.4.1) H(w) = S*(w) / (|S(w)|2 +Wq (w)). (12.4.2) При некоррелированной помехе с дисперсией s2 система уравнений для определения значений коэффициентов h(n): ho (R(0)+s2 )+ h1 R(1)+ h2 R(2)+ h3 R(3)+ ...+ hM R(M) = s(0), (12.4.3) ho R(1) + h1 R(0)+ h2 R(1)+ h3 R(2)+ ...+ hM R(M-1) = 0, ho R(2) + h1 R(1)+ h2 R(0)+ h3 R(1)+ ...+ hM R(M-2) = 0, . . . . . . . . . . . . . ho R(M) + h1 R(M-1)+ h2 R(M-2)+ ....... + hM R(0) = 0. При расчете коэффициентов фильтра значение s(0) обычно принимается произвольным, чаще всего равным площади сигнала s(t). Тем самым делается попытка полного сжатия площади сигнала до весовой функции Кронекера, что возможно только для сигналов со спектром в главном диапазоне до частоты Найквиста. Рис. 12.4.1. Сжатие гладких сигналов с разным уровнем шумов. Для гладких и монотонных функций со спектром в низкочастотной части главного диапазона сжатие до импульса Кронекера невозможно, и в зависимости от уровня шумов фильтр поднимает насколько возможно высокие частоты сигнала, учитывая значение уровня шумов. При этом нарушаются условия нормировки площади оператора фильтра к 1, о чем можно судить по значению передаточной функции H(w) при w=0, которое становится меньше 1, и при обратном преобразовании H(w) Þ h(m) оператор h(m) требует нормировки к 1. Все эти факторы можно наглядно видеть на рис. 12.4.1. На рисунке приведены три сигнала с одной и той же базовой функцией, на которую наложены шумы разного уровня. При малом уровне шумов (сигнал х1) фильтр в максимальной степени использует высокие частоты сигнала (|H1| >>1 на этих частотах), сохраняя устойчивость работы фильтра при достаточно удовлетворительном (хотя и больше 1) коэффициенте усиления дисперсии помех при максимально возможном сжатии сигнала. При повышении уровня шумов (сигналы х2 и х3) подъем высоких частот сигнала уменьшается, и сжатие сигнала соответственно также уменьшается, предпочтение отдается максимальному подавлению шумов. Рис. 12.4.2. Сжатие сигнала с высокочастотным спектром На рис. 12.4.2. приведен пример сжатия сигнала, близкого к прямоугольному импульсу. Базовая функция сигнала s(k) имеет достаточно высокочастотный спектр мощности Ws(w), и при задании формы выходного сигнала сжатия в виде гауссовой функции z(k) передаточная функция фильтра H(w) обеспечивает уверенное сжатие сигнала (при уменьшении уровня шумов практически до заданной формы). В пределе, при Wq=0 фильтр сжатия превращается в фильтр деконволюции: H(w)= S*(w) / |S(w)|2 = 1/S(w), (12.4.4) На выходе такого фильтра имеем: Y(w) = H(w)X(w) → 1, при X(w) → S(w). Реализация фильтра возможна только при условии S(w) > 0 на всех частотах в главном частотном диапазоне. В противном случае, при S(wi ) → 0, H(wi ) → ∞ и фильтр становится неустойчивым. Для исключения возможности такого явления в фильтр (12.4.4) вводится стабилизатор a: H(w) = S*(w) / [|S(w)|2 +a], (12.4.5) где |S(w)|2 +a > 0 во всем частотном диапазоне. Фильтры деконволюции могут использоваться не только для повышения разрешающей способности данных, но и для интерпретации геофизических данных, если формирование полезного входного сигнала удовлетворяет принципу суперпозиции данных по зависимости от искомых параметров. 12.5. Фильтр обнаружения сигналов /42/. Фильтр используется при решении задач обнаружении сигналов известной формы на существенном уровне шумов, значение которых может превышать значения сигналов. В процессе фильтрации необходимо зафиксировать наличие сигнала в массиве данных, если он там присутствует (может не присутствовать), при этом сохранения формы сигнала не требуется. Сама форма сигнала полагается известной либо по теоретическим данным, либо по результатам измерений на моделях. Для уверенного обнаружения сигнала фильтр должен обеспечить максимально возможную амплитуду выходного сигнала над уровнем помех и выполняется на основе критерия максимума пикового отношения сигнал/шум. Частотная характеристика. Для расчета фильтра требуется задать известную форму полезного сигнала s(k) - S(w) и функцию автокорреляции или спектр мощности помех Rq (m) - Wq (w). Полный входной сигнал принимается по аддитивной модели: x(t) = s(t)+q(t). На выходе проектируемого фильтра h(n) - H(w) для составляющих выходного сигнала имеем: y(t) = H(w) S(w) exp(jwt) dw, (12.5.1) s2 = |H(w)|2 Wq (w) dw, (12.5.2) где s - средняя квадратическая амплитуда выходной помехи. Оптимальным в задаче обнаружения одиночного сигнала конечной длительности является фильтр, обеспечивающий на выходе максимальное отношение пиковой мощности сигнала к мощности шума в момент окончания импульса. Значения (12.5.1, 12.5.2) используются для задания критерия максимума отношения сигнал/шум (12.2.3) для произвольной точки ti : r = [y(ti )]2 /d2 . (12.5.3) Исследование функции (12.5.3) на максимум показывает, что он достигается при частотной характеристике фильтра: H(w) = exp(-jwti ) |S*(w)| / Wq (w), (12.5.4) Для физически реализуемых фильтров в качестве точки ti целесообразно использовать интервал длительности импульса t, при этом: H(w) = exp(-jwt) |S*(w)| / Wq (w) = exp(-jjs -jwt) |S(w)|/Wq (w). (12.5.4') Аргумент js в этом выражении компенсирует фазовые сдвиги составляющих спектра сигнала, а wt обеспечивает их задержку на время длительности сигнала. Таким образом, на концевой части сигнала фильтр выполняет синфазное суммирование всех полезных частотных составляющих входного сигнала с весами, пропорциональными отношению |S(w)|/Wq (w), что обеспечивает накопление амплитуды полезного сигнала на интервале всей длительности входного импульса и формирует максимум сигнала на момент его окончания. Вместе с тем фильтр ослабляет спектральные составляющие шума тем сильнее, чем меньше модуль |S(w)|, и полная мощность шума на выходе фильтра оказывается меньшей, чем на входе. Для получения линейных уравнений расчета коэффициентов фильтра без потери общности можно принять ti =0, при этом: H(w) = S*(w)/Wq (w) = |S(w)|exp(jjs (w)) / Wq (w). (12.5.5) При переходе во временную (координатную) область: H(w)Wq (w) = S*(w) - h(n) ③ Rq (n-m) = s(-m). (12.5.6) Система линейных уравнений для расчета фильтра: ho Rq (0)+ h1 Rq (1)+ h2 Rq (2)+ h3 Rq (3)+ ...+ hM Rq (M) = S(-M), ho Rq (1)+ h1 Rq (0)+ h2 Rq (1)+ h3 Rq (2)+ ...+ hM Rq (M-1)= S(-M+1), ho Rq (2)+ h1 Rq (1)+ h2 Rq (0)+ h3 Rq (1)+ ...+ hM Rq (M-2)= S(-M+2), . . . . . . . . . . . . . ho Rq (M)+ h1 Rq (M-1)+ h2 Rq (M-2)+ ..... + hM Rq (0) = S(0). При задании ti по центру симметричных входных сигналов можно получить симметричные двусторонние фильтры, не изменяющие фазы сигнала. На рис. 12.5.1 приведен пример фильтрации фильтром обнаружения сигнала радиоимпульса (информационный сигнал) в сумме с шумами (входной сигнал) при отношении сигнал/шум по средним амплитудным значениям на входе фильтра порядка 1. Аналогичное отношение сигнал/шум на выходе фильтра повышается до 7 по интервалу полезного сигнала в целом, и превышает 8 в центральной части интервала сигнала. Рис. 12.5.1. Фильтр обнаружения сигнала. Эффективность фильтра. Из выражения (12.5.5) можно видеть, что фильтр имеет максимальный коэффициент передачи на частотах доминирования сигнала и минимальный коэффициент передачи на частотах доминирования помех. Синфазность суммирования всех частотных составляющих выходного сигнала обеспечивает максимальную амплитуду выходного сигнала в заданный момент времени ti . Значение максимальной амплитуды можно оценить, приняв ti =0, при этом выходной сигнал: y(0)- S(w) H(w) = = . Коэффициент передачи фильтра прямо определяется спектром подлежащего обнаружению сигнала, его формой и длительностью. Для оценки эффективности фильтра зададим входной сигнал в виде прямоугольного импульса амплитудой u0 длительностью t на интервале 0-t. Спектральная плотность прямоугольного импульса при интегральном преобразовании Фурье: П(w) = (1-exp(-jwt))/jw, П*(w) = (exp(jwt)-1)/jw. При подстановке в (12.5.4'), принимая Wq (w) = const, коэффициент передачи фильтра: H(w) = a[(exp(jwt)-1)exp(-jwt]/jw = a(1-exp(-jwt))/jw, где a - коэффициент пропорциональности с размерностью, обратной спектральной плотности, для получения безразмерных значений коэффициента H(w). При a=1 (нормировка оператора фильтра производится, как правило, по коэффициенту усиления постоянной составляющей входного сигнала) сигнал на выходе фильтра: uвых (t) = (u0 /2p) П(w)H(w) dw = (u0 /2p) (1-exp(-jwt))2 exp(jwt) dw, uвых (t) = U0 {t|t >0 – 2(t-t)|t > t + (t-2t)|t >2 t }. Рис. 12.5.2. Как можно видеть на рис 12.5.2, выходной сигнал для входного прямоугольного импульса представляет собой треугольный импульс длительностью 2t по основанию с максимальным значением амплитуды на концевой части входного импульса. Это определяется тем, что при Wq (w)=1 оператор фильтра полностью повторяет форму входного сигнала (прямоугольного импульса), а выходной сигнал в отсутствие шумов представляет собой свертку двух одинаковых импульсов, максимальное значение которой достигается при полном входе сигнала в интервал оператора фильтра (t=t) и равно полной энергии входного импульса: U0 = п(t)·h(t) dt = п(t)2 dt = u0 2 ·t. Значение U0 определяется нормировкой оператора фильтра a. Что касается усиления дисперсии (мощности) шумов, то, как известно, дисперсия шума на выходе фильтра равна дисперсии входных шумов s2 , умноженной на интеграл квадрата импульсного отклика фильтра (для цифровых систем – сумма квадратов коэффициентов оператора фильтра): s2 вых = s2 h2 (t) dt = (s2 /2p) |H(w|2 dw. Для вычисления интеграла модуль передаточной функции фильтра для прямоугольного импульса может быть представлен в виде интегрального синуса: |H(w|2 dw = 2t u0 2 sinc2 (wt/2) d(wt/2) = 2p u0 2 t. Дисперсия шумов на выходе: s2 вых = s2 u0 2 t. С использованием этого выражения для отношения мощности сигнала к мощности шума для сигналов на входе и выходе фильтра имеем: rвх = u0 2 /s2 , rвых = u0 4 t2 /s2 u0 2 t = u0 2 t/s2 . Для отношения амплитудных значений сигнала к среднеквадратическим значениям шума: rвх = u0 /s, rвых = (u0 /s) . Отсюда следует, что эффективность фильтра тем выше, чем больше длительность взаимодействия сигнала с оператором фильтра. Фильтр жестко настраивается под форму сигнала, и любое изменение формы сигнала понижает его эффективность. Отметим также, что коэффициент передачи фильтра тем больше и эффективность его работы тем выше, чем больше различия в форме частотных спектров сигнала и шумов. При постоянной форме спектров сигнала и шума любой другой фильтр уступает данному фильтру, как по пиковому, так и по энергетическому отношению сигнал/шум на выходе фильтра. Согласованный фильтр. При помехах типа белого шума Wq (w) = s2 и из (12.5.5) следует H(w) = S*(w)/s2 . Постоянный множитель 1/s2 может быть опущен. Частотная характеристика фильтра с точностью до постоянного множителя определяется только спектром сигнала, вследствие чего он и получил название согласованного (со спектром сигнала). Рис. 12.5.3. В качестве примера на рис. 12.5.3 приведен радиоимпульсный сигнал s(t), форма которого хорошо известна априорно, и зашумленный входной сигнал y(t) = s(t) + q(t), при этом мощность шумового сигнала q(t) в 2 раза больше мощности радиоимпульса. Рис.12.5.4. На рис. 12.5.4 приведены модули спектров S(w) ↔ s(t), S* (w) и аргументы этих спектров F(w) и F* (w) в правой половине главного частотного диапазона. Согласованный фильтр имеет АЧХ, подобную фильтру сигнала и фазово-частотную характеристику, комплексно сопряженную с ФЧХ сигнала. При этом оператор фильтра определяется выражением: S* (w) ↔ s(-t’) = h(t’), (12.5.7) что представляет собой реверсирование известной формы сигнала в интервале его значимых значений (для данного примера на рис. 12.5.3 t’=0…100 из интервала t=150…250) Рис. 12.5.5 На рис. 12.5.5 приведены модули спектров исходного и сигнала и выходного сигнала согласованного фильтра: u(t) = y(t) ③ h(t’) = s(t) ③ h(t’) + q(t) ③ h(t’). (12.5.8) Первый член этой формулы s(t)③h(t’), с учетом реверса оператора h(t’) = s(-t’) при свертке относительно s(t’), это функция автокорреляции сигнала s(t’), сдвинутая относительно начала сигнала s(t) в 12.5.8 на длину оператора h(t’), т.е. с максимумом на концевой части сигнала s(t). Что касается второго члена формулы, то свертка q(t) ③ h(t’) представляет собой интегрирование по интервалу оператора случайных шумовых компонент, умножаемых на весовую функцию h(t’) c нулевым средним значением, что приводит к их уменьшению, хотя спектральная характеристика этих шумов перекрывается со спектром выделяемого сигнала. Рис. 12.5.6. Пример согласованной фильтрации приведен на рис. 12.5.6. Для исключения сдвига по фазе выходного сигнала относительно входного можно выполнять сдвиг оператора фильтра влево по координатам на половину интервала длины оператора (второй график свертки на рисунке). Фильтр мало эффективен при выделении коротких импульсных или длинных гармонических сигналов. Обратный фильтр. Допустим, что помеха имеет такой же частотный состав, что и полезный сигнал, т.е.: Wq = s2 |S(w)|2 . Выделение полезного сигнала в таких условиях сомнительно. Тем не менее, определим оптимальный фильтр: H(w) = S*(w) / [s2 |S(w)|2 ] = 1 / [s2 S(w)]. (12.5.9) Выражение (12.5.9) с точностью до постоянного множителя соответствует фильтру сжатия сигнала. Но если согласованный фильтр и фильтр сжатия рассматривать в качестве предельных случаев при полной неопределенности характеристики помех, то в качестве модели помех можно принять их суперпозицию: Wq = a2 |S(w)|2 +b2 . Подставляя это выражение в (12.5.5), с точностью до множителя получаем: H(w) = S*(w) / [|S(w)|2 +g2 ], (12.5.10) где g = b/a - отношение дисперсий шума и сигнала. Фильтр стремится к согласованному при больших g, и к обратному (фильтру сжатия) при малых. 12.6. Энергетический фильтр. Энергетический фильтр максимизирует отношение сигнал/помеха по всей длине фильтра (а не в отдельной точке), и если сигнал по своей протяженности укладывается в окно фильтра, то тем самым обеспечивается оценка формы сигнала. Фильтр занимает промежуточное положение между фильтром воспроизведения сигнала Колмогорова- Винера и согласованным фильтром и требует задания корреляционных функций сигнала и помех. Сигнал может быть представлен и в детерминированной форме с соответствующим расчетом его автокорреляционной функции. Критерий оптимальности. Энергия сигнала на выходе фильтра: Esh = Sk sk 2 = Sk (Sn hn sk- n )2 = Sk hk Sn hn Rs (k-n), (12.6.1) где Rs - функция автокорреляции сигнала. В векторной форме: Esh = . (12.6.2) Аналогично, выражение для энергии помех на выходе: Eqh = Sk hk Sn hn Rq (k-n) = , (12.6.3) где Rq - функция автокорреляции помех. При некоррелированной помехе Eqh = s2 . Подставим (12.6.2, 12.6.3) в выражение (12.2.4): r = / . (12.6.4) Расчет векторов операторов фильтров. Для определения значений вектора продифференцируем r по , и приравняем производную к нулю: . . (12.6.5) В системе уравнений (12.6.5) неизвестны собственные значения r матрицы и значения коэффициентов hn . Система имеет N+1 ненулевых решений относительно значений r и соответствующих этим значениям векторов . Для определения коэффициентов фильтра приравнивается к нулю и решается относительно r определитель матрицы , после чего максимальное значение rmax подставляется в (12.6.5) и система уравнений решается относительно коэффициентов hi вектора . При фильтрации сигнала вектор обеспечивает выделение первой по мощности главной компоненты сигнала, т.е. составляющей сигнала, которая имеет наибольшую энергию и отношение сигнал/шум. В сложных полях такая компонента, как правило, соответствует региональному фону. В принципе, расчет может быть продолжен и для других значений r<rmax , и определены значения коэффициентов векторов , и т.д., с использованием которых могут выделяться вторая и прочие компоненты сигнала. Наиболее эффективно такой метод используется для разделения сигналов (полей) при некоррелированных помехах. В этом случае корреляционная матрица помех является единичной (единицы по диагонали, остальное - нули) и уравнение (12.6.5) имеет вид: . (12.6.6) В развернутой форме: ho (Rs (0)-r)+ h1 Rs (1)+ h2 Rs (2)+ h3 Rs (3)+ ...+ hM Rs (M) = 0, ho Rs (1)+ h1(Rs(0)-r)+ h2 Rs (1)+ h3 Rs (2)+ ...+ hM Rs (M-1) = 0, ho Rs (2)+ h1 Rs (1)+ h2 (Rs(0)-r)+ h3 Rs (1)+ ...+ hM Rs (M-2) = 0, . . . . . . . . . . . . . ho Rs (M)+ h1 Rs (M-1)+ h2 Rs (M-2)+ ..... + hM (Rs(0)-r) = 0. Выражение (12.6.6) при малом уровне шумов позволяет вместо ФАК какого-либо определенного сигнала использовать ФАК непосредственно зарегистрированных данных. Если при этом в зарегистрированных данных кроме помех присутствуют два (и более) сигналов, например, региональный фон и локальная составляющая (аномалия), то расчет векторов hi приобретает конкретный практический смысл. После первой фильтрации оператором и выделения региональной составляющей, массив данных (исходный или с вычитанием из него региональной составляющей) может быть профильтрован повторно оператором , что позволит выделить и локальную аномалию (и т.д.). Разделение сигналов будет тем надежнее, чем сильнее они отличаются друг от друга по энергии и интервалу корреляции. В заключение отметим, что расчеты оптимальных фильтров могут производиться с использованием алгоритма Левинсона. литература 12. Канасевич Э.Р. Анализ временных последовательностей в геофизике. - М.: Недра, 1985.- 300 с. 42. Яневич Ю.М. Задачи приема сигналов и определения их параметров на фоне шумов: Курс лекций. / СПбУ. А.В.Давыдов. 07.02.10. Cайт автора Лекции Практикум О замеченных ошибках и предложениях по дополнению: davpro@yandex.ru.