Пример инверсии оператора через спектральное представление. Определение сигнала на входе линейной системы по значениям его выходного типа. Особенности выполнения деконволюции. Вычисление коэффициентов инверсного фильтра по значениям каузального оператора.
Пример инверсии оператора через спектральное представление приведен на рис. Функция H(z) в выражении (3) имеет особые точки - нули функции, которые становятся полюсами функции Так, например, фильтр, реализующий передаточную функцию (5), в самой общей форме может быть выполнен в виде включенных последовательно фильтров, каждый из которых имеет передаточную функцию следующего типа (для первого фильтра): H1-1(z) = 1/(a - z) = (1 z/a z2/a2 ...)/a. Коэффициенты при степенях 1/z являются, соответственно, коэффициентами инверсного фильтра с координатами (-n), т.е. фильтр оперирует с "будущими" отсчетами входного сигнала (см. рис. Вычисление коэффициентов инверсного фильтра по значениям каузального (одностороннего) оператора h[n] может быть проведено на основе выражения (2): h-1[k] h[n - k] = o(n), (6) для чего достаточно развернуть его в систему n-уравнений при = 0,1,2…, k ? n: n = 0: h-1[0] h[0] = 1, h-1[0] = 1/h[0]. n = 1: h-1[0] h[1] h-1[1]h[0] = 0, h-1[1] = h-1[0]h[1] / h[0]. n = 2: h-1[0]h[2] h-1[1]h[1] h-1[2]h[0] = 0, h-1[2] = {h-1[0]h[2] h-1[1]h[1]}/h[0], и т.д.
Введение
Деконволюция или обратная свертка. Применяется для сжатия сигналов с целью повышения временного или пространственного разрешения результатов измерений.
Понятие деконволюции
Определение. Прямая свертка сигнала x[k] c импульсным откликом h[n] линейной системы определяется как: y(k) = h(n) * x(k) o H(z)X(z) = Y(z).
Обратная задача, задача деконволюции - определение сигнала на входе линейной системы по значениям выходного сигнала: X(z) = Y(z)/H(z) = Y(z)H -1(z) o y[k] * h -1[n] = x[k]. (1) где * - оператор свертки.
Прямая деконволюция предполагает, что выполняются соотношения: H(z)H-1(z) = 1 o h[n] * h-1[n] = ?o(n), (2)
H-1(z) = 1/H(z) o h-1[n]. (3) где o(n) - импульс Кронекера (?o(n) = 1 при n = 0, ?o(n) = 0 при n ? 0).
Рис. 1.
Пример инверсии оператора через спектральное представление приведен на рис. 1 ,где h[n] - исходный оператор; H(?) - спектральная плотность; H-1(?) - инверсная спектральная плотность; h-1 [n] инверсный оператор на начальном интервале отсчетов.
Свойства деконволюции
Выражение (3) позволяет сделать некоторые выводы об особенностях выполнения деконволюции.
При ограниченной импульсной реакции h[n] инверсный оператор h-1[n] в общем случае не ограничен. Так, например, если импульсная реакция равна h[n] = {1,a} o (1 az) = H(z), то имеем: H-1(z) = 1/(1 az) = 1- az a2z2 - a3z3 .... h-1[n] = {1, -a, a2, -a3,....}.
Это действительно практически для любых операторов фильтров, энергия которых на каких-либо ограниченных участках главного частотного диапазона близка к нулевой. При инверсии спектральной функции таких операторов на этих участках возникают резкие энергетические пики, которые при обратном преобразовании Фурье дает медленно затухающие функции операторов. Пример такого явления приведен на рис. 2.
Рис. 2.
Отсюда следует, что для точного выполнения деконволюции необходимо располагать бесконечно длинным инверсным оператором фильтра.
Инверсия операторов обычно всегда связана с усилением высоких частот, что приводит к резкому повышению коэффициента усиления дисперсии помех инверсных фильтром.
Деконволюция выполняется, если инверсный оператор достаточно быстро затухает и может быть ограничен. Использование усеченных операторов приводит к появлению определенной погрешности деконволюции, величину которой следует контролировать.
Устойчивость фильтров деконволюции. Функция H(z) в выражении (3) имеет особые точки - нули функции, которые становятся полюсами функции
H-1(z) = 1/H(z) и определяют устойчивость инверсного фильтра. Для того чтобы фильтр деконволюции был устойчивым, ряд 1/H(z) должен сходиться, т.е. полюса функции должны находиться вне единичного круга на z-плоскости.
Многочлен H(z) порядка N может быть разложен на N простых сомножителей - двучленов (диполей): H(z) = (а - z)(b - z)(c - z)., (4) где а, b, с, - корни полинома. Обращение передаточной функции: H-1(z) = (5)
Функции первого типа. Если каждый из двучленов функции (4) является минимально-фазовой функцией, т.е. ее корни находится вне единичного круга на z-плоскости и модули нулевых членов всегда больше следующих за ними первых членов (в данном случае: |а| > 1, |b| > 1, |с |> 1), то и функция H(z) в целом также является минимально-фазовой функцией. При этом максимум энергии импульсного отклика сосредоточен в его начальной части и последовательность отсчетов представляет собой затухающий ряд. Соответственно, и функция 1/H(z) также будет представлять собой сходящийся ряд, и инверсный фильтр будет устойчив. Так, например, фильтр, реализующий передаточную функцию (5), в самой общей форме может быть выполнен в виде включенных последовательно фильтров, каждый из которых имеет передаточную функцию следующего типа (для первого фильтра): H1-1(z) = 1/(a - z) = (1 z/a z2/a2 ...)/a.
Модули всех корней больше 1, следовательно, полюсы инверсного полинома будут находиться за пределами единичной окружности на z-плоскости, и инверсный оператор устойчив. Форма исходного оператора и положение полюсов инверсного оператора на z-плоскости приведены на рис. 3.
Рис. 3.
Пример. Проверка устойчивости инверсного фильтра.
Оператор фильтра h[n] такой же как в предыдущей задаче, но сдвинут вправо на один отсчет и дополнен новым нулевым отсчетом h[0] = 0.048.
Форма исходного оператора и положение полюсов приведены на рис.4.
Рис. 4.
Два полюса находятся внутри единичного круга, следовательно, инверсный оператор не будет устойчивым (будет незатухающим расходящимся рядом).
Обращение функций второго типа. Если функция H(z) имеет корни составляющих его двучленов, которые находятся внутри и на единичном круге в z-плоскости, то устойчивое обращение H(z) имеет отрицательные степени z и для его использования необходимо располагать "будущими" значениями входного сигнала.
Пример.
Передаточная функция фильтра имеет вид: H(z) = 1 - 2 z.
Инверсная функция представляется как H-1(z) = 1 / (1 - 2 z).
Частотные спектры функций приведены на рис. 5.
Рис. 5.
Полюс функции zp = 1/2 и находится внутри единичного круга на z-плоскости.
Перепишем выражение для инверсного фильтра в следующем виде: H-1(z) = - (1 / 2z) [1 1/2z 1/(2z)2 ...].
Это выражение является разложением в ряд по степеням 1/z и сходится к кругу радиусом 1/2 при z > ?. Коэффициенты при степенях 1/z являются, соответственно, коэффициентами инверсного фильтра с координатами (-n), т.е. фильтр оперирует с "будущими" отсчетами входного сигнала (см. рис. 5).
Если двучлены (4) представляют собой функции и первого и второго типа то обращение определяется полным рядом Лорана: H-1(z) = ... h[-2] z- 2 h[-1] z -1 h[0] h[1] z 1 h[2] z 2 ..., т.е. оператор инверсного фильтра является двусторонним и не обязательно симметричным.
Вычисление коэффициентов инверсного фильтра по значениям каузального (одностороннего) оператора h[n] может быть проведено на основе выражения (2): h-1[k] h[n - k] = o(n), (6) для чего достаточно развернуть его в систему n-уравнений при = 0,1,2…, k ? n: n = 0: h-1[0] h[0] = 1, h-1[0] = 1/h[0]. n = 1: h-1[0] h[1] h-1[1]h[0] = 0, h-1[1] = h-1[0]h[1] / h[0]. n = 2: h-1[0]h[2] h-1[1]h[1] h-1[2]h[0] = 0, h-1[2] = {h-1[0]h[2] h-1[1]h[1]}/h[0], и т.д.
Продолжая последовательно, можно вычислить любое количество значений коэффициентов инверсного фильтра. Рекуррентная формула вычислений: h-1[n] = - [ h-1[k] h[n - k]] / h[0]. (7)
Если фильтр деконволюции устойчив и ряд h-1[n] сходится, то появляется возможность ограничения количества членов ряда с определенной ошибкой восстановления исходного сигнала. Метрика приближения Е (квадратичная норма разности) определяется выражением: Е2 = [?o(n) - (h[n] * h-1[n])]2. (8)
Ошибка восстановления исходного сигнала проявляется со сдвигом на длину прямого оператора фильтра.
Модули корней больше 1, инверсный фильтр должен быть устойчивым и быстро затухающим.
3. Двенадцать первых значений инверсного оператора при вычислении по (7): h-1[n]= {4.56, -6.033, 2.632, 0.417, -0.698, -0.062, 0.267, -0.024, -0.11, 0.051, 0.018, -0.019, 0.004}.
Значения прямого и инверсного оператора фильтра приведены на рис. 6
4. Значения свертки прямого оператора с инверсным при разных длинах N инверсного фильтра и метрики приближения имеют вид: N = 4, n= {1, 0, 0, 0, 0, 0.014, -0.04, -0.056, -0.028, -0.01, -0.002, 0, 0, …}, E = 0.077;
На рис. 7 приведены абсолютные значения ошибки деконволюции при разной длине N.
Рис. 7
Оптимальные фильтры деконволюции
Оптимальные фильтры деконволюции имеют метрики приближения много меньше, чем усеченные фильтры. Для получения общего уравнения оптимальной деконволюции будем считать, что число коэффициентов оператора h[n] равно M 1, a число коэффициентов инверсного оператора h[n]-1 равно N 1.
Выходная функция приближения при использовании уравнения свертки (2) с ограничением числа членов оператора фильтра равна: F = Е2 = [o(k) - x[k] ]2. x[k] = h[n] -1 h[k-n] (9)
Чтобы определить минимум функции, приравняем нулю частные производные от Е по неизвестным коэффициентам фильтра: DF/d h[j] -1 = -2 h[k - 1] [o(k) - h[n]-1 h[k-n] ] = 0. (10) h[k-j] h[n]-1 h[k-n] = h[k-j] o(k) = h[-j] - (11) hn-1 h[k-n] h[k-j] = h[n]-1 aj-n = h[-j], j = 0,1,2, ..., N, (12) где aj-n - функция автоковариации импульсной реакции h[n].
Учитывая также, что h[n] = 0 при n < 0 и aj = a-j (функция автоковариации является четной функцией), окончательное решение определяется следующей системой линейных уравнений: a0 h[0]-1 a1 h[1]-1 a2 h[2]-1 a3 h[3]-1 ... AN h[N]-1= h[0] (13) a1 h [0] -1 a0 h [1] -1 a1 h [2] -1 a2 h [3] -1 ... AN-1 h [N] -1= 0 a2 h [0] -1 a1 h [1] -1 a0 h [2] -1 a1 h [3] -1 ... AN-2 h [N] -1= 0
Метрика приближения оптимального оператора в 1.6 раза меньше усеченного.
Рис. 8
Как видно на рис. 8, оптимизация инверсного оператора заключается в центрировании ошибок приближения и тем самым уменьшении максимально возможных ошибок.
Уравнение оптимальной инверсии. Оптимальный инверсный фильтр может быть получен непосредственно с использованием z-образов импульсной реакции и автоковариационной функции прямого фильтра.
Для прямого фильтра z-образ автоковариационной функции фильтра (как z-отображение спектральной плотности мощности) представляет собой произведение: A(z) = H(z)H*(z), (14) где H*(z) -функция, комплексно сопряженная с передаточной функцией прямого фильтра H(z).
Заменяя H(z) для функции первого типа выражением H(z) = 1/H-1(z), получаем: А(z)H-1(z) = H*(z). (15)
Запишем последнее равенство в развернутом виде: (а-N z -N ... a-1 z -1 a0 a1 z 1 ... AN z N)?
?(h[0] -1 h[1] -1 z 1 h [2] -1 z 2 ... h [N] -1 z N ) =
= h[0]* h [1]* z -1 h [2] * z -2 ... h [N] * z N. (16)
В выражении (16) сумма коэффициентов при одинаковых степенях z в левой части равенства должна быть равна коэффициенту при соответствующей степени z в правой части равенства, что позволяет составить следующую систему из N уравнений для коэффициентов при степенях z0, z1, z2, ... , ZN: a 0 h [0] -1 a -1 h [1] -1 a-2 h [2] -1 a-3 h [3] -1 ... a -N h [N] -1 = h [0] *
(17) a 1 h [0] -1 a0 h [1] -1 a -1 h [2] -1 a -2 h [3] -1 ... a -N -1 h [N] -1 = 0 a 2 h [0 ]-1 a 1 h [1] -1 a 0 h [2] -1 a 5 h [3]-1 ... a -N -2 h [N] -1 = 0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a N h [0] -1 a N -1 h [1] -1 a N - 2 h [2]-1 a N - 3 h [3] -1 ... a 0 h [-N ] -1= 0
В случае вещественных фильтров, когда а i = a -i и h [0]* = h [0], уравнение (17) идентично уравнению (13).
Уравнение Левинсона
Перепишем уравнение (17) в матричной форме: (18)
Так как коэффициенты инверсного фильтра достаточно определить с точностью до произвольного масштабного множителя, приведем h [0]-1 к 1, a функцию автоковариации переведем в функцию коэффициентов корреляции делением обеих частей уравнения на ао. Обозначая
Аі = аі/ao, Wi = h [i]-1 / h [0] -1 и V = ho*/(h o -1 ao) = h o h o * / ao , получаем: (19) где для значений W и V введен индекс j номеров предстоящих итераций по циклу вычисления коэффициентов фильтра.
При нулевой итерации (N = 0, j = 0) имеем только одно уравнение: . (20)
Благодаря проведенной нормировке решения уравнения (20) не требуется: А0= 1, V0= 1, W00= 1.
Перепишем уравнение (20) в прямой форме: А0 W00 = V0. (22)
Запишем вспомогательную систему, для чего к уравнению (22) добавим вторую строку с постоянной Ej: A0 W00 A1·0 = V0, A1 W00 A0 ·0 = E1.
В матричной форме: (23)
Преобразуем уравнение (23): (24)
Вычтем (24) из (23) с неопределенным множителем Rj: (25)
Из верхней строки уравнения (24) можно получить значение Е1: Е1= A1W00. (26)
Уравнение (21) можно сделать равнозначным уравнению (25), если правую часть нижней строки уравнения (25) приравнять правой части нижней строки уравнения (21): E1 - R1V0 = 0, R1 = E1/V0. (27)
Из правых частей нижней и верхней строк уравнений (30, 33): R2 = E2 / V1, V2 = V1 - R2E2.
Новые коэффициенты из левых частей уравнений (30, 33): W02 = W01 - R2 ? 0= 1, W12 = W11 - R2W11, W22 = 0 - R2W01.
Анализ расчетов позволяет вывести следующие рекуррентные формулы: Ej = Ai Wj-i,j , j = 1,2,..., M. (34)
Rj = Ej / V j-1, Vj = Vj-1 - Rj Ej, Wi,j = Wi,j-1 - RJWJ-1,j-1, i = 0,1,.., j.
Рекурсивная деконволюция
Запишем уравнение (3) для инверсного фильтра в развернутой форме: H-1(z) = 1/(h [0] h [1] z h [2] z2 ...). (35)
Так как для минимально-фазового оператора всегда выполняется условие h[0]? 0, приведем (35) к виду: H-1(z) = (1/ h [0])/(1 h [1] z/ h [0] h [2] z2/ h [0] ...) = G / (1 q1z q2z2 ...), (36) где: G = 1/ h [0], q1 = h [1] / h [0], q2 = h [2] / h [0] и т.д.
Но уравнение (36) есть не что иное, как уравнение передаточной функции рекурсивного фильтра, где цепь обратной связи образована коэффициентами нормированного оператора h[n]. Алгоритм вычислений: yk = G·xk - qn·yk-n.
Выражение (36) уникально по своим возможностям. В принципе, оно может реализовать оператор инверсной фильтрации с бесконечным импульсным откликом. На практике оно может использоваться вместо медленно затухающих инверсных операторов, модуль одного из полюсов которого очень близок к 1 (менее 1.1) при высоких требованиях к метрике приближения.
Пример рекурсивной деконволюции.
Рис.9
Оператор фильтра h [n] = {0.41, 0.791, 0.401, -0.193, -0.367, -0.166, 0.032, 0.068, 0.027, -0.001}, N = 9.
1. Модуль одного из корней фильтра равен 1.032, что приводит к очень слабому затуханию инверсного оператора. Метрика приближения даже при N = =100 для усеченного оператора составляет 0.3. Форма операторов приведена на рис. 9.
Рис.10
2. При использовании оптимального инверсного оператора с N=100 значение погрешности приближения уменьшается более чем в 20 раз, что позволяет уменьшить длину оператора до N=35 при погрешности приближения порядка 0.1 (рис 7.4.2(А)), при этом абсолютные значения погрешностей приближения не превышают 0.03 (рис. 10 (В)).
3. Расчет коэффициентов фильтра рекурсивной деконволюции: G = 1/ h [0] = 2.441 gn = h [n] G. gn = {1.932, 0.978, -0.472, -0.896, -0.405, 0.077, 0.165, 0.065, -0.003}, n =1, 2, 3, …, 9.
Рис. 11
На рис. 11 приведены результаты рекурсивной деконволюции оператора h[n]. Как и следовало ожидать, деконволюция абсолютно точно, с нулевой метрикой, восстанавливает импульс Кронекера, хотя собственный импульсный отклик рекурсивного оператора повторяет оператор h[n]-1 при его вычислении по формуле (7) и длительность его значимой части близка к 200. Коэффициент усиления дисперсии шумов при данной операции вычисляется по значениям импульсного отклика и весьма существенен.
Фильтры неполной деконволюции
Рассмотренные выше методы относятся к расчету так называемых идеальных инверсных фильтров, т.е. фильтров полной прямой инверсии системных операторов. Однако использование идеальных инверсных фильтров на практике не всегда возможно, т.к. регистрируемые данные обычно осложнены влиянием помех (шумов), а инверсные фильтры обычно имеют коэффициент усиления шумов значительно больше 1. В этом случае задача точной деконволюции (восстановления истинной первоначальной формы сигнала), как правило, не ставится, и инверсные фильтры считаются оптимальными с точки зрения максимального приближения к форме полезного сигнала с определенным допустимым коэффициентом усиления дисперсии помех. Такие фильтры называются фильтрами неполной (частичной, ограниченной) деконволюции. При проектировании фильтров неполной деконволюции учитываются статистические характеристики помех во входном сигнале и их соотношение со статистическими характеристиками самого входного сигнала. Передаточная функция фильтра неполной деконволюции с учетом помех во входном сигнале определяется выражением: H-1(z) = H*(z)/[|H(z)|2 g2], (36) где g2 = k·h2 - дисперсия шумов в единицах дисперсии оператора h[n] h2 - дисперсия значений оператора h[n] (при условии суммы значений оператора, равной 1), k - отношение дисперсии шумов к дисперсии оператора h[n].
Коэффициент g2 играет роль регуляризирующего фактора при выполнении операции деконволюции информации.
Рис. 12
На рис. 12 пример формы оператора h[n] и спектральных функций (36) при разных значениях параметра g. При g = 0 выражение (36) обращается в идеальный инверсный фильтр 1/H(z). Во втором крайнем случае, при g2 >> |H(z)|2, фильтр (36) переходит в фильтр, согласованный с сигналом по частотному спектру: H-1(z) = H*(z)/g2, который только максимизирует отношение сигнал/помеха.
На рис. 13 приведена форма инверсных операторов, соответствующая их частотным характеристикам на рис. 12(В), и результаты свертки инверсных операторов с прямым (для лучшего просмотра графики прямой оператор при свертке сдвинут вправо на 2 значения ?t). При g = 0 коэффициент усиления дисперсии шумов равен 11, при g = 0.4 ?h2 равен 4.6. Однако снижение усиления дисперсии шумов сопровождается увеличением погрешности приближения, что можно видеть на рис. 13(В), при этом уменьшается амплитуда восстановления импульса Кронекера и появляются осцилляции после импульса. Но при наличии шумов и правильном выборе параметра g общее отношение амплитудных значений сигнал/ шум для оператора по (36) больше, чем для прямой инверсии по (3), что объясняется более существенным уменьшением коэффициента усиления дисперсии шумов при увеличении параметра g, чем увеличением погрешности приближения.
Рис. 13
Операторы оптимальных фильтров также могут вычисляться с учетом помех. Если сигнал s(k) и помеха статистически независимы, то функция автоковариации сигнала на входе фильтра: аі = asi bi, (37) где asi и bi - функции автоковариации сигнала и помех.
При помехе типа белого шума функция автоковариации помех представляет собой весовую дельта-функцию в точке 0: bi = c2i, (38) где с2- дисперсия помех.
С учетом этого фактора расчет оптимальных инверсных фильтров может проводиться по вышеприведенным формулам (13, 17) с изменением значения коэффициента ао:
ao= ao c2. (39)
Рис. 14
На рис. 14(А) приведены примеры операторов оптимальных инверсных фильтров, вычисленные по прямому оператору, приведенному на рис. 12(А). Значения коэффициента с2 заданы в долях дисперсии прямого оператора. деконволюция сигнал инверсия фильтр
