Тема 13. Деконволюция цифровых сигналов если дом красив, то мы понимаем, что он был выстроен для хозяев, а не для мышей (стр. 2 из 4)

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

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

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

Рис. 13.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.

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

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

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

Если диполи функции (13.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). Частотные спектры функций приведены на рис. 13.1.5.

Рис. 13.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), т.е. фильтр оперирует с "будущими" отсчетами входного сигнала (см. рис. 13.1.5).

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

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

h^-1(k)h(n-k) = d_o(n), (13.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). (13.2.2)

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

Е² =

[d_o(n) - h(n) ③ h^-1(n)]². (13.2.3)

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

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

Рис. 13.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), быстро затухает.

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

Значения прямого и инверсного оператора фильтра приведены на рис. 13.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.

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

Рис. 13.2.2.

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

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

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

F = Е² =

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

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

dF/dh_j^-1 =

-2h_k-j [d_o(k) - h_n^-1 h_k-n] = 0. (13.3.2)

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

h_n^-1 h_k-n h_k-j = h_n^-1 a_j-n = h_-j, j = 0,1,2, ..., N, (13.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₀ (13.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

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

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

Рис. 13.3.1.

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

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

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

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

Оператор по (13.2.2) – рис. 13.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.

Оператор по (13.3.5) – рис. 13.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 раза меньше усеченного.

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

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