Лекции по "Технологии цифровой обработки сигналов"

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

3.1. Понятие деконволюции.

При достаточно сложном физическом представлении во временной (координатной) области деконволюция проста для понимания в частотном представлении. Допустим, что в какой-либо регистрирующей системе происходит резонансное поглощение энергии и сдвиг по фазе определенного гармонического колебания в составе входного сигнала, например, преобразование гармоники sin w₀t ® 0.5 sin (w₀-p/4). Соответственно, для восстановления первоначальной формы сигнала операция деконволюции должна выполнить усиление данной гармоники в выходном сигнале в 2 раза и осуществить обратный сдвиг фазы на p/4. Для одной гармоники выполнение такой операции труда не представляет. Но практические задачи деконволюции значительно сложнее, так как требуется, как правило, восстановление полного спектра исходного сигнала, имеющего непрерывный характер. 

Определение деконволюции. Если для прямой свертки дискретного сигнала x(k) c импульсным откликом h(n) линейной системы (фильтра) имеем уравнение:

y(k) = h(n) ③ x(k) ó H(z)X(z) = Y(z),   Y(z) =

y_k z^k,  z = exp(-jj),

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

X(z) = Y(z)/H(z) = Y(z)H^-1(z) ó y(k) ③ h^-1(n) = x(k),                 (3.1.1)

где индексом "-1" символически обозначена передаточная функция оператора  обратного фильтра H^-1(z) = 1/H(z), т.е. оператора деконволюции, инверсного прямому оператору свертки (импульсному отклику системы). Соответственно, при обратном z-преобразовании получаем оператор деконволюции:

H^-1(z) = 1/H(z) ó h^-1(n).                                        (3.1.2)

Очевидно, что если имеет  место H(z)H^-1(z) = 1, то при обратном z-преобразовании этого выражения мы должны иметь:

H(z)H^-1(z) = 1 ó h(n) ③ h^-1(n) = d_o(n),                             (3.1.3)

где d_o(n) - импульс Кронекера. При этом последовательная свертка сигнала x(k) с оператором системы h(k) и оператором деконволюции h^-1(k) будет эквивалентна свертке сигнала x(k) с импульсом Кронекера, что не должно изменить форму сигнала x(t).

При z = exp(-jj) все изложенное действительно и для спектрального представления операторов. Пример инверсии оператора через спектральное представление приведен на рис. 3.1.1 (исходный оператор h_n ® спектральная плотность H(ω) ® инверсная спектральная плотность H^-1(ω) ® инверсный оператор h^-1_n на начальном интервале отсчетов).

Рис. 3.1.1.

Особенности деконволюции. Выражение (3.1.2) позволяет сделать некоторые выводы об особенностях выполнения деконволюции.

При ограниченной импульсной реакции h(n) инверсный оператор h^-1(n) в общем случае не ограничен. Так, например, если импульсная реакция представлена нормированным диполем h(n) = {1,a} ó (1+az) = h(z), то имеем:

H^-1(z) = 1/(1+az) = 1-az+a²z²-a³z³+ ....

h^-1(n) = {1, -a, a², -a³,....}.

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

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

Рис. 3.1.2.

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

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

Устойчивость  фильтров деконволюции. Функция H(z) в выражении (3.1.2) имеет особые точки - нули функции, которые становятся полюсами функции H^-1(z) = 1/H(z) и определяют устойчивость инверсного фильтра. Для того чтобы фильтр деконволюции был устойчивым, ряд 1/H(z) должен сходиться, т.е. полюса функции должны находиться вне единичного круга на z-плоскости (внутри круга при использовании символики z^-1).     

Многочлен H(z) порядка N может быть разложен на N простых сомножителей - двучленов (диполей):

H(z) = (a-z)(b-z)(c-z)....,                                       (3.1.4)

где а, b, с,… - корни полинома. Обращение передаточной функции:

H^-1(z) =                                      (3.1.5)

Если каждый из диполей  функции (3.1.4) является минимально-фазовым  диракоидом, т.е. корни диполей находится  вне единичного круга на z-плоскости и модули нулевых членов диполей всегда больше следующих за ними первых членов (в данном случае: |а|>1, |b|>1, |с|>1), то функция H(z) также является минимально-фазовым диракоидом. Максимум энергии импульсного отклика минимально-фазового диракоида сосредоточен в его начальной части и последовательность отсчетов представляет собой затухающий ряд. При этом функция 1/H(z) также будет представлять собой передаточную функцию минимально-фазового оператора, обеспечивающую выполнение условия (3.1.3), а ее обратное преобразование - сходящийся ряд устойчивого фильтра. Так, например, фильтр, реализующий передаточную функцию (3.1.5), в самой общей форме может быть выполнен в виде включенных последовательно фильтров, каждый из которых имеет передаточную функцию следующего типа (для первого фильтра):

H₁^-1(z) = 1/(a-z) = (1+z/a+z²/a²+...)/a.

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

    Пример 1.  Оценить возможность инверсии оператора свертки h_n = {0.131, 0.229, 0.268, 0.211, 0.111, 0.039, 0.009, 0.001}, N = 7.

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

    На рисунке  показан также оператор деконволюции, который получен обратным преобразованием  Фурье передаточной функции Hi(w). Оператор бесконечен, но затухает достаточно быстро, что обеспечивает высокую точность деконволюции при использовании ограниченного числа членов оператора фильтра (устанавливается пользователем по заданной точности).

Рис. 3.1.3. Пример устойчивого  инверсного фильтра.

Пример 2.  Сдвинем вышеприведенный оператор на одну позицию вправо и для сохранения той же энергии оператора примем h₀ = 0.001. Требуется оценить устойчивость инверсного оператора свертки h_n = {0.001, 0.131, 0.229, 0.268, 0.211, 0.111, 0.039, 0.009, 0.001}, N = 7.

Рис. 3.1.4. Пример неустойчивого инверсного фильтра

Результаты проверки на устойчивость нового оператора деконволюции приведены на рис. 3.1.4. При сходной  форме операторов свертки модули спектров операторов практически не отличаются. Но за счет сдвига по фазе данного оператора относительно первого (на рис. 3.1.3) корни его z-полинома претерпели существенное изменение, а модуль одного из корней меньше 1. И хотя вычисленный оператор деконволюции также представляет собой ряд с конечной энергией, но условие (3.1.3) не выполняется, о чем наглядно свидетельствует результат свертки h_{k *} h^-1_k.   

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

Если диполи функции (3.1.4) представляют собой и диракоиды, и реверсоиды, то обращение будет  центроидом и определяется полным рядом  Лорана:

H^-1(z) = ...+h_-2z^-2+h_-1z^-1+h₀+ h₁z¹+h₂z²+ ...,

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

   Пример 3.  Передаточная функция фильтра: H(z) = 1-2z.   Инверсная функция H^-1(z) = 1/(1-2z). Частотные спектры функций приведены на рис. 3.1.5.

Рис. 3.1.5.

   Полюс функции  z_p = 1/2 и находится внутри единичного круга на z-плоскости.

   Перепишем выражение для  инверсного фильтра в следующем виде:

H^-1(z) = -(1/2z) [1+1/2z+1/(2z)²+...].

  Это выражение  является разложением в ряд  по степеням 1/z и сходится к  кругу радиусом 1/2 при z → ¥. Коэффициенты при степенях 1/z являются, соответственно, коэффициентами инверсного фильтра с координатами (-n), т.е. фильтр оперирует с "будущими" отсчетами входного сигнала (см. рис. 3.1.5).

3.2. Инверсия импульсного отклика  фильтра.

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

h^-1(k)h(n-k) = d_o(n),                                       (3.2.1)

для чего достаточно развернуть его в систему n-уравнений при 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) = (1/h₀) h^-1(k)h(n-k).                                (3.2.2)

Если фильтр деконволюции устойчив и ряд h^-1(n) сходится, то появляется возможность ограничения количества членов ряда с определенной ошибкой восстановления исходного сигнала. Метрика приближения Е (квадратичная норма разности) определяется выражением:

Е² = [d_o(n) - h(n) ③ h^-1(n)]².                                    (3.2.3)

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

   Пример инверсии оператора фильтра  h_n = {0.219, 0.29, 0.257, 0.153, 0.061, 0.016, 0.003}.

Рис. 3.2.1.

Полином по zⁿ:    H(z) = S_n h_n zⁿ.

Модули корней полинома:  z_n = {2.054, 2.054, 2.485, 2.485, 1.699, 1.699}.

Модули корней больше 1, инверсный фильтр устойчив и, судя по значениям корней (достаточно существенно отличаются от 1), быстро затухает.

Двенадцать значений оператора по (3.2.2): 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}.

Значения прямого и  инверсного оператора фильтра приведены  на рис. 3.2.1.

Значения свертки прямого  оператора с инверсным при  длине N=10 инверсного фильтра и метрика приближения:  s_n={1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.004, 0.006, 0.004, 0.002, <0.001, …}. E=0.0086.

На рис. 3.2.2 приведены  графики абсолютных значений концевой ошибки деконволюции при разной длине N оператора деконволюции.

         Рис. 3.2.2.

3.3. Оптимальные фильтры деконволюции  /12,22/.

Можно рассчитать оптимальные  фильтры деконволюции, метрика приближения  которых меньше, чем у усеченных  фильтров деконволюции. Для получения общего уравнения оптимальной деконволюции будем считать, что число коэффициентов оператора h_n равно M+1, a число коэффициентов инверсного оператора h_n^-1 равно N+1.

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

F = Е² = [d_o(k)-s_k]².      s_k = h_n^-1 h_k-n.                    (3.3.1)

Чтобы определить минимум  функции, приравняем нулю  частные  производные от Е по неизвестным коэффициентам фильтра:

dF/dh_j^-1 = -2h_k-j [d_o(k) - h_n^-1 h_k-n] = 0.                    (3.3.2)

h_k-j h_n^-1 h_k-n = h_k-j d_o(k) = h_-j.                       (3.3.3)

h_n^-1 h_k-n h_k-j  = h_n^-1 a_j-n = h_-j,    j = 0,1,2, ..., N,            (3.3.4)

где a_j-n - функция автоковариации импульсной реакции h(n). Учитывая также, что h_n = 0 при n<0 и a_j = a_-j (функция автоковариации является четной функцией), окончательное решение определяется следующей системой линейных уравнений:

a₀ h₀^-1 + a₁ h₁^-1 + a₂ h₂^-1 + a₃ h₃^-1 + ...  + a_N h_N^-1 = h₀                (3.3.5)

a₁ h₀^-1 + a₀ h₁^-1 + a₁ h₂^-1 + a₂ 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

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

  Повторим инверсию  оператора, приведенного на рис. 3.2.1, N=6.

Рис. 3.3.1.

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