ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ
Digital signals processing. Deconvolution of signals
Тема 13. ДЕКОНВОЛЮЦИЯ ЦИФРОВЫХ СИГНАЛОВ
Если дом красив, то мы понимаем, что он был выстроен для хозяев, а не для мышей.
Хрисипп. Древнегреческий философ. III в.д.н.э.
Деконволюция сходна с археологией. Задача – восстановить дом из развалин. Обнадеживает, если обломки на месте. Но если только песок, то понять, был ли тут дом для людей или кладовка для мышей, не так-то просто.
Григорий Старцев. Уральский геофизик. ХХ в.
Содержание:
Введение.
1. Понятие деконволюции. Определение деконволюции. Особенности деконволюции. Устойчивость фильтров деконволюции. Обращение недиракоидных функций.
2. Инверсия импульсного отклика фильтра. Вычисление коэффициентов инверсного фильтра. Пример инверсии оператора фильтра.
3. Оптимальные фильтры деконволюции. Принцип оптимизации. Пример расчета оптимального фильтра деконволюции. Уравнение оптимальной инверсии. Уравнение Левинсона.
4. Рекурсивная деконволюция. Уравнение фильтра рекурсивной деконволюции. Пример рекурсивной деконволюции.
5. Фильтры сжатия сигналов. Передаточная функция фильтра. Операторы оптимальных фильтров.
введение
Основное назначение деконволюции (deconvolution) – восстановление истинной формы сигнала, несущего информацию об исследуемом физическом или технологическом процессе, явлении природы и т.п., после его искажения при регистрации какой-либо линейной системой - измерительным трактом прибора (аппаратной или приборной функцией) или каналом связи. Естественно, что для восстановления необходимы сведения о характеристиках искажающей системы, и в первую очередь, об импульсном отклике системы или его частотной передаточной функции. Для выполнения деконволюции реализуются фильтры, частотные характеристики которого обратны частотной характеристике искажающих систем. Построение таких фильтров не всегда возможно. Так, невозможно в принципе восстановить в сигнале частоты, которые были полностью подавлены, а при восстановлении частотных составляющих, ослабленных до уровня шумов, одновременно происходит значительное усиление дисперсии шумов, в которых полезный сигнал может полностью затеряться.
Вместе с тем, деконволюция или обратная свертка используется и для решения других задач обработки данных. Так, в геофизике она применяется для сжатия сигналов с целью повышения временного или пространственного разрешения результатов измерений. В грави- и магниторазведке с использованием деконволюции производятся перерасчеты аномальных полей вниз. В ядерной геофизике методы деконволюции являются основными при количественной интерпретации результатов измерений, чему способствует принцип суперпозиции ядерно-физических полей.
При неизвестных частотных характеристиках искажающей системы речь может идти о слепой деконволюции (blind deconvolution). Слепая деконволюция - намного более сложная проблема, которая не имеет общего решения и разрабатывается с учетом определенных априорных сведений для специфических приложений.
13.1. Понятие деконволюции.
При достаточно сложном физическом представлении во временной (координатной) области деконволюция проста для понимания в частотном представлении. Допустим, что в какой-либо регистрирующей системе происходит резонансное поглощение энергии и сдвиг по фазе определенного гармонического колебания в составе входного сигнала, например, преобразование гармоники sin w0t ® 0.5 sin (w0-p/4). Соответственно, для восстановления первоначальной формы сигнала операция деконволюции должна выполнить усиление данной гармоники в выходном сигнале в 2 раза и осуществить обратный сдвиг фазы на p/4. Для одной гармоники выполнение такой операции труда не представляет. Но практические задачи деконволюции значительно сложнее, так как требуется, как правило, восстановление полного спектра исходного сигнала, имеющего непрерывный характер.
Определение деконволюции. Если для прямой свертки дискретного сигнала x(k) c импульсным откликом h(n) линейной системы (фильтра) имеем уравнение:
y(k) = h(n) ③ x(k) - H(z)X(z) = Y(z), Y(z) =
yk zk, z = exp(-jj),то задача деконволюции в общей форме - определение сигнала на входе линейной системы по значениям выходного сигнала, т.е. снятие реакции (импульсного отклика) системы на сигнал и восстановление исходной формы сигнала, что весьма актуально, например, для регистрирующих систем:
X(z) = Y(z)/H(z) = Y(z)H-1(z) - y(k) ③ h-1(n) = x(k), (13.1.1)
где индексом "-1" символически обозначена передаточная функция оператора обратного фильтра H-1(z) = 1/H(z), т.е. оператора деконволюции, инверсного прямому оператору свертки (импульсному отклику системы). Соответственно, при обратном z-преобразовании получаем оператор деконволюции:
H-1(z) = 1/H(z) - h-1(n). (13.1.2)
Очевидно, что если имеет место H(z)H-1(z) = 1, то при обратном z-преобразовании этого выражения мы должны иметь:
H(z)H-1(z) = 1 - h(n) ③ h-1(n) = do(n), (13.1.3)
где do(n) - импульс Кронекера. При этом последовательная свертка сигнала x(k) с оператором системы h(k) и оператором деконволюции h-1(k) будет эквивалентна свертке сигнала x(k) с импульсом Кронекера, что не должно изменить форму сигнала x(t).
При z = exp(-jj) все изложенное действительно и для спектрального представления операторов. Пример инверсии оператора через спектральное представление приведен на рис. 13.1.1 (исходный оператор hn ® спектральная плотность H(ω) ® инверсная спектральная плотность H-1(ω) ® инверсный оператор h-1n на начальном интервале отсчетов).
Рис. 13.1.1.
Особенности деконволюции. Выражение (13.1.2) позволяет сделать некоторые выводы об особенностях выполнения деконволюции.
При ограниченной импульсной реакции h(n) инверсный оператор h-1(n) в общем случае не ограничен. Так, например, если импульсная реакция представлена нормированным диполем h(n) = {1,a} - (1+az) = h(z), то имеем:
H-1(z) = 1/(1+az) = 1-az+a2z2-a3z3+ ....
h-1(n) = {1, -a, a2, -a3,....}.
Это действительно практически для любых операторов свертки, энергия которых на каких-либо ограниченных участках главного частотного диапазона близка к нулевой. При инверсии спектральной функции таких операторов на этих участках возникают резкие спектральные пики, которые при обратном преобразовании Фурье дают медленно затухающие функции операторов. Пример такого явления приведен на рис. 13.1.2.
Отсюда следует, что для точного выполнения деконволюции необходимо располагать бесконечно длинным инверсным оператором фильтра. Практически, деконволюция выполняется, если инверсный оператор достаточно быстро затухает и может быть ограничен. Но использование усеченных операторов приводит к появлению определенной погрешности деконволюции, величину которой следует контролировать.
Рис. 13.1.2. |
Заметим также, что передаточные функции систем, как правило, имеют низкочастотный характер. Инверсия операторов таких систем всегда связана с усилением высоких частот, что приводит к высоким коэффициентам усиления инверсными фильтрами дисперсии помех, что может приводить к потере полезного сигнала среди усиленных шумовых флюктуаций.
Кроме того, быстрое затухание оператора деконволюции является необходимым, но еще не достаточным условием реализации деконволюции.
Устойчивость фильтров деконволюции. Функция H(z) в выражении (13.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)...., (13.1.4)
где а, b, с,… - корни полинома. Обращение передаточной функции:
H-1(z) =
(13.1.5)Если каждый из диполей функции (13.1.4) является минимально-фазовым диракоидом, т.е. корни диполей находится вне единичного круга на z-плоскости и модули нулевых членов диполей всегда больше следующих за ними первых членов (в данном случае: |а|>1, |b|>1, |с|>1), то функция H(z) также является минимально-фазовым диракоидом. Максимум энергии импульсного отклика минимально-фазового диракоида сосредоточен в его начальной части и последовательность отсчетов представляет собой затухающий ряд. При этом функция 1/H(z) также будет представлять собой передаточную функцию минимально-фазового оператора, обеспечивающую выполнение условия (13.1.3), а ее обратное преобразование - сходящийся ряд устойчивого фильтра. Так, например, фильтр, реализующий передаточную функцию (13.1.5), в самой общей форме может быть выполнен в виде включенных последовательно фильтров, каждый из которых имеет передаточную функцию следующего типа (для первого фильтра):
H1-1(z) = 1/(a-z) = (1+z/a+z2/a2+...)/a.
Отсюда, равно как и непосредственно из выражения (13.1.5), следует также, что чем больше значения модулей корней фильтра (дальше от единичного круга полюс фильтра), тем быстрее затухает инверсный оператор.