Пример расчета оптимального фильтра деконволюции.

Повторим инверсию оператора, приведенного в последнем примере (рис. 7.2.1 и 7.2.2), при N=6.

1. Значения оптимального инверсного оператора в сопоставлении с усеченным:

h^-1(n) = {4.56, -6.033, 2.632, 0.417, -0.698, -0.062, 0.267} – прямой расчет по (7.2.2).

h^-1(n) = {4.557, -6.026, 2.633, 0.397, -0.693, -0.009, 0.145} – расчет по (7.3.5).

2. Значения свертки инверсных операторов с прямыми и метрики приближения:

Рис. 7.3.1.

Оператор по (7.2.2) – рис. 7.3.1(А): s_n= {1, 0, 0, 0, 0, 0, 0, 0.005, 0.031, 0.027, 0.013, 0.004, 0.001, 0, 0,…}. E=0.044.

Оператор по (7.3.5) – рис. 7.3.1(В): s_n= {0.999, <0.001, 0.002, -0.003, -0.003, 0.013, -0.008, -0.012, 0.011, 0.013, 0.007, 0.002, <0.001, 0, 0, …}. E=0.027.

Метрика приближения оптимального оператора в 1.6 раза меньше усеченного.

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

Уравнение оптимальной инверсии. Оптимальный инверсный фильтр может быть получен непосредственно с использованием z-образов импульсной реакции и автоковариационной функции прямого фильтра. Если для прямого фильтра мы имеем передаточную функцию H(z), то z-образ автоковариационной функции фильтра (как z-отображение спектральной плотности мощности) представляет собой произведение:

A(z) = H(z)H*(z), (7.3.6)

где H*(z)- функция, комплексно сопряженная с H(z). Заменяя H(z) для функции диракоидного типа выражением H(z) = 1/H^-1(z), получаем:

А(z)H^-1(z) = H*(z). (7.3.7)

Запишем последнее равенство в развернутом виде:

(а_-Nz^-N+ ... +a_-1z^-1+a₀+a₁z¹+ ... +a_Nz^N)(h₀^-1+h₁^-1z¹+h₂^-1z²+ ... +h_N^-1z^N) =

= h₀*+h₁*z^-1+h₂*z^-2+ ... +h_N*z^N. (7.3.8)

В выражении (7.3.8) сумма коэффициентов при одинаковых степенях z в левой части равенства должна быть равна коэффициенту при соответствующей степени z в правой части равенства, что позволяет составить следующую систему из N уравнений для коэффициентов при степенях z⁰, z¹, z², ... , z^N:

(7.3.9)

a₀ h₀^-1 + a_-1 h₁^-1 + a_-2 h₂^-1 + a_-3 h₃^-1 + ... + a_-N h_N^-1 = h₀*

a₁ h₀^-1 + a₀ h₁^-1 + a_-1 h₂^-1 + a_-2 h₃^-1 + ... + a_-N-1 h_N^-1 = 0

a₂ h₀^-1 + a₁ h₁^-1 + a₀ h₂^-1 + a₅ h₃^-1 + ... + a_-N-2 h_N^-1 = 0

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

a_N h₀^-1 + a_N-1 h₁^-1 + a_N-2 h₂^-1 + a_N-3 h₃^-1 +... +a₀ h_-N ^-1 = 0

В случае вещественных фильтров, когда а_i = a_-i и h₀* = h₀, уравнение (7.3.9) идентично уравнению (7.3.5).

Уравнение Левинсона.Практический способ расчета оптимальных инверсных фильтров по уравнению (7.3.9) предложен в 1947 году Н.Левинсоном.

Перепишем уравнение (7.3.9) в матричной форме:

(7.3.10)

Так как коэффициенты инверсного фильтра достаточно определить с точностью до произвольного масштабного множителя, приведем h_o^-1 к 1, a функцию автоковариации переведем в функцию коэффициентов корреляции делением обеих частей уравнения на а_o. Обозначая А_i = а_i/a_o, W_i = h_i^-1/h_o^-1 и V = h_o*/(h_o^-1a_o) = h_oh_o*/a_o , получаем:

(7.3.11)

где для значений W и V введен индекс j номеров предстоящих итераций по циклу вычисления коэффициентов фильтра.

При нулевой итерации (N=0, j=0) имеем только одно уравнение:

. (7.3.12)

Благодаря проведенной нормировке решения уравнения (7.3.12) не требуется:

А₀= 1, V₀= 1, W₀₀= 1.

Составим уравнение для двучленного фильтра (N=1, j=1):

(7.3.13)

Перепишем уравнение (7.3.12) в прямой форме:

А₀ W₀₀ = V₀. (7.3.14)

Запишем вспомогательную систему, для чего к уравнению (7.3.14) добавим вторую строку с новой постоянной Е_j:

A₀ W₀₀ + A₁·0 = V_0,

A₁ W₀₀ + A₀ ·0 = E₁.

В матричной форме:

(7.3.15)

Реверсируем уравнение (7.3.15):

(7.3.16)

Вычтем (7.3.16) из (7.3.15) с неопределенным множителем R_j:

(7.3.17)

Из верхней строки уравнения (7.3.16) можно получить значение Е₁:

Е₁= A₁W₀₀. (7.3.18)

Уравнение (7.3.13) можно сделать равнозначным уравнению (7.3.17), если правую часть нижней строки уравнения (7.3.17) приравнять правой части нижней строки уравнения (7.3.13):

E₁ - R₁V₀ = 0

R₁ = E₁/V₀. (7.3.19)

При этом из правых частей верхних строк уравнений (7.3.13,7.3.17) будем иметь:

V₁ = V₀ - R₁E₁. (7.3.20)

Приравнивая друг другу левые части уравнений (7.3.13,7.3.17), получаем:

W₀₁ = W₀₀ - R₁·0 = W₀₀ = 1.

W₁₁ = 0 - R₁W₀₀ = -R₁W₀₀. (7.3.21)

Этим заканчивается первая итерация. Аналогично, для второй итерации:

(7.3.22)

(7.3.23)

(7.3.24)

(7.3.25)

Из верхней строки уравнения (7.3.24):

Е₂ = A₁W₁₁+A₂W₀₁.

Из правых частей нижней и верхней строк уравнений (7.3.22,7.3.25):

R₂ = E₂/V₁,

V₂ = V₁ - R₂E₂.

Новые коэффициенты из левых частей уравнений (7.3.22,7.3.25):

W₀₂ = W₀₁ - R₂ 0= 1,

W₁₂ = W₁₁ - R₂W₁₁,

W₂₂ = 0 - R₂W₀₁.

Анализ расчетов позволяет вывести следующие рекуррентные формулы:

E_j = A_iW_j-i,_j , j = 1,2,...,M. (7.3.26)

R_j = E_j/V_j-1,

V_j = V_j-1 - R_jE_j,

W_i,j = W_i,j-1 - R_jW_j-1,_j-1, i = 0,1,.., j.

Подпрограммы решения уравнений для ЭВМ приведены в литературе /12,22/.

7.4. Рекурсивная деконволюция /л22/.

Запишем уравнение (7.1.3) для инверсного фильтра в развернутой форме:

H^-1(z) = 1/(h₀+h₁z+h₂z²+ ...). (7.4.1)

Так как для минимально-фазового оператора всегда выполняется условие h₀ 0, приведем (7.4.1) к виду:

H^-1(z) = (1/h₀)/(1+h₁z/h₀+h₂z²/h₀+...) = G/(1+q₁z+q₂z²+ ...), (7.4.2)

где: G = 1/h₀, q₁ = h₁/h₀, q₂ = h₂/h₀ и т.д. Но уравнение (7.4.2) есть не что иное, как уравнение передаточной функции рекурсивного фильтра, где цепь обратной связи образована коэффициентами нормированного оператора h(n). Алгоритм вычислений:

y_k = G·x_k – q_n·y_k-n.

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