ЦИФРОВАЯ
ОБРАБОТКА
СИГНАЛОВ
Digital signals processing. Deconvolution of signals
Тема 13. ДЕКОНВОЛЮЦИЯ ЦИФРОВЫХ СИГНАЛОВ
Если дом красив, то мы понимаем, что он был выстроен для хозяев, а не для мышей.
Хрисипп. Древнегреческий философ. III в.д.н.э.
Деконволюция сходна с археологией. Задача – восстановить дом из развалин. Обнадеживает, если обломки на месте. Но если только песок, то понять, был ли тут дом для людей или кладовка для мышей, не так-то просто.
Григорий Старцев. Уральский геофизик. ХХ в.
Содержание:
Введение.
1.
Понятие деконволюции.
Определение деконволюции. Особенности деконволюции. Устойчивость фильтров деконволюции. Обращение недиракоидных функций.
2. Инверсия импульсного отклика фильтра. Вычисление коэффициентов инверсного фильтра. Пример инверсии оператора фильтра.
3. Оптимальные фильтры деконволюции. Принцип оптимизации. Пример расчета оптимального фильтра деконволюции.
Уравнение оптимальной инверсии. Уравнение Левинсона.
4. Рекурсивная деконволюция. Уравнение фильтра рекурсивной деконволюции. Пример рекурсивной деконволюции.
5. Фильтры сжатия сигналов. Передаточная функция фильтра. Операторы оптимальных фильтров.
введение
Основное назначение деконволюции (deconvolution
) – восстановление истинной формы сигнала, несущего информацию об исследуемом физическом или технологическом процессе, явлении природы и т.п., после его искажения при регистрации какой-либо линейной системой - измерительным трактом прибора (аппаратной или приборной функцией) или каналом связи. Естественно, что для восстановления необходимы сведения о характеристиках искажающей системы, и в первую очередь, об импульсном отклике системы или его частотной передаточной функции. Для выполнения деконволюции реализуются фильтры, частотные характеристики которого обратны частотной характеристике искажающих систем. Построение таких фильтров не всегда возможно. Так, невозможно в принципе восстановить в сигнале частоты, которые были полностью подавлены, а при восстановлении частотных составляющих, ослабленных до уровня шумов, одновременно происходит значительное усиление дисперсии шумов, в которых полезный сигнал может полностью затеряться.
Вместе с тем, деконволюция или обратная свертка используется и для решения других задач обработки данных. Так, в геофизике она применяется для сжатия сигналов с целью повышения временного или пространственного разрешения результатов измерений. В грави- и магниторазведке с использованием деконволюции производятся перерасчеты аномальных полей вниз. В ядерной геофизике методы деконволюции являются основными при количественной интерпретации результатов измерений, чему способствует принцип суперпозиции ядерно-физических полей.
При неизвестных частотных характеристиках искажающей системы речь может идти о слепой деконволюции
(blind deconvolution). Слепая деконволюция - намного более сложная проблема, которая не имеет общего решения и разрабатывается с учетом определенных априорных сведений для специфических приложений.
13.1. Понятие деконволюции.
При достаточно сложном физическом представлении во временной (координатной) области деконволюция проста для понимания в частотном представлении. Допустим, что в какой-либо регистрирующей системе происходит резонансное поглощение энергии и сдвиг по фазе определенного гармонического колебания в составе входного сигнала, например, преобразование гармоники sin w0
t ® 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-1
n
на начальном интервале отсчетов).
Рис. 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+a2
z2
-a3
z3
+ ....
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), следует также, что чем больше значения модулей корней фильтра (дальше от единичного круга полюс фильтра), тем быстрее затухает инверсный оператор.
Пример 1.
Оценить возможность инверсии оператора свертки hn
= {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.
Сдвинем вышеприведенный оператор на одну позицию вправо и для сохранения той же энергии оператора примем h0
= 0.001. Требуется оценить устойчивость инверсного оператора свертки hn
= {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) не выполняется, о чем наглядно свидетельствует результат свертки hk
*
h-1
k
.
Обращение недиракоидных функций.
Если H(z) - реверсоид, т.е. корни составляющих его диполей находятся внутри и на единичном круге в z-плоскости, то устойчивое обращение H(z) является антиимпульсом (с отрицательными степенями z) и для его использования необходимо располагать "будущими" значениями входного сигнала.
Если диполи функции (13.1.4) представляют собой и диракоиды, и реверсоиды, то обращение будет центроидом и определяется полным рядом Лорана:
H-1
(z) = ...+h-2
z-2
+h-1
z-1
+h0
+ h1
z1
+h2
z2
+ ...,
т.е. оператор инверсного фильтра является двусторонним и не обязательно симметричным, как мы обычно рассматривали ранее двусторонние операторы.
Пример 3.
Передаточная функция фильтра: H(z) = 1-2z. Инверсная функция H-1
(z) = 1/(1-2z). Частотные спектры функций приведены на рис. 13.1.5.
|
Рис. 13.1.5. |
Полюс функции zp
= 1/2 и находится внутри единичного круга на z-плоскости.
Перепишем выражение для инверсного фильтра в следующем виде:
H-1
(z) = -(1/2z) [1+1/2z+1/(2z)2
+...].
Это выражение является разложением в ряд по степеням 1/z и сходится к кругу радиусом 1/2 при z → ¥. Коэффициенты при степенях 1/z являются, соответственно, коэффициентами инверсного фильтра с координатами (-n), т.е. фильтр оперирует с "будущими" отсчетами входного сигнала (см. рис. 13.1.5).
13.2. Инверсия импульсного отклика фильтра.
Вычисление коэффициентов инверсного фильтра
по значениям каузального (одностороннего) оператора h(n) может быть проведено на основе выражения (13.1.3):
h-1
(k)h(n-k) = do
(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/h0
)h-1
(k)h(n-k). (13.2.2)
Если фильтр деконволюции устойчив и ряд h-1
(n) сходится, то появляется возможность ограничения количества членов ряда с определенной ошибкой восстановления исходного сигнала. Метрика приближения Е (квадратичная норма разности) определяется выражением:
Е2
=[do
(n) - h(n) ③ h-1
(n)]2
. (13.2.3)
Ошибка восстановления исходного сигнала появляется со сдвигом на длину оператора фильтра деконволюции и проявляется на интервале длины прямого оператора.
Пример инверсии оператора фильтра
hn
= {0.219, 0.29, 0.257, 0.153, 0.061, 0.016, 0.003}.
|
Рис. 13.2.1. |
Полином по zn
: H(z) = Sn
hn
zn
.
Модули корней полинома: zn
= {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 инверсного фильтра и метрика приближения: sn
={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/.
Можно рассчитать оптимальные фильтры деконволюции, метрика приближения которых меньше, чем у усеченных фильтров деконволюции. Для получения общего уравнения оптимальной деконволюции будем считать, что число коэффициентов оператора hn
равно M+1, a число коэффициентов инверсного оператора hn
-1
равно N+1.
Принцип оптимизации.
Выходная функция приближения при использовании уравнения свертки (13.1.2) с ограничением числа членов оператора фильтра:
F = Е2
= [do
(k)-sk
]2
. sk
=hn
-1
hk-n
. (13.3.1)
Чтобы определить минимум функции, приравняем нулю частные производные от Е по неизвестным коэффициентам фильтра:
dF/dhj
-1
= -2hk-j
[do
(k) -hn
-1
hk-n
] = 0. (13.3.2)
hk-j
hn
-1
hk-n
= hk-j
do
(k) = h-j
. (13.3.3)
hn
-1
hk-n
hk-j
= hn
-1
aj-n
= h-j
, j = 0,1,2, ..., N, (13.3.4)
где aj-n
- функция автоковариации импульсной реакции h(n). Учитывая также, что hn
= 0 при n<0 и aj
= a-j
(функция автоковариации является четной функцией), окончательное решение определяется следующей системой линейных уравнений:
a0
h0
-1
+ a1
h1
-1
+ a2
h2
-1
+ a3
h3
-1
+ ... +
hN
-1
= h0
(13.3.5)
a1
h0
-1
+ a0
h1
-1
+ a1
h2
-1
+ a2
h3
-1
+ ... + aN-1
hN
-1
= 0
a2
h0
-1
+ a1
h1
-1
+ a0
h2
-1
+ a1
h3
-1
+ ... + aN-2
hN
-1
= 0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
aN
h0
-1
+ aN-1
h1
-1
+ aN-2
h2
-1
+ aN-3
h3
-1
+... +a0
hN
-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(А): sn
= {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(В): sn
= {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-отображение спектральной плотности мощности) представляет собой произведение:
A(z) = H(z)H*(z), (13.3.6)
где H*(z)- функция, комплексно сопряженная с H(z). Переносим H(z) в левую часть формулы и для функций диракоидного типа заменяем выражение 1/H(z) = H-1
(z):
А(z)H-1
(z) = H*(z). (13.3.7)
Запишем последнее равенство в развернутом виде:
(a-N
z-N
+ ... +a-1
z-1
+a0
+a1
z1
+ ... +aN
zN
)(h0
-1
+h1
-1
z1
+h2
-1
z2
+ ... +hN
-1
zN
) =
= h0
*+h1
*z-1
+h2
*z-2
+ ... +hN
*zN
. (13.3.8)
В выражении (13.3.8) сумма коэффициентов при одинаковых степенях z в левой части равенства должна быть равна коэффициенту при соответствующей степени z в правой части равенства, что позволяет составить следующую систему из N уравнений для коэффициентов при степенях z0
, z1
, z2
, ... , zN
:
a0
h0
-1
+ a-1
h1
-1
+ a-2
h2
-1
+ a-3
h3
-1
+ ... + a-N
hN
-1
= h0
* (13.3.9)
a1
h0
-1
+ a0
h1
-1
+ a-1
h2
-1
+ a-2
h3
-1
+ ... + a-N-1
hN
-1
= 0
a2
h0
-1
+ a1
h1
-1
+ a0
h2
-1
+ a5
h3
-1
+ ... + a-N-2
hN
-1
= 0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
aN
h0
-1
+ aN-1
h1
-1
+ aN-2
h2
-1
+ aN-3
h3
-1
+... +a0
h-N
-1
= 0
В случае вещественных фильтров, когда ai
= a-i
и h0
* = h0
, уравнение (13.3.9) идентично уравнению (13.3.5).
Уравнение Левинсона.
Практический способ расчета оптимальных инверсных фильтров по уравнению (13.3.9) предложен в 1947 году Н.Левинсоном.
Перепишем уравнение (13.3.9) в матричной форме:
(13.3.10)
Так как коэффициенты инверсного фильтра достаточно определить с точностью до произвольного масштабного множителя, приведем ho
-1
к 1, a функцию автоковариации переведем в функцию коэффициентов корреляции делением обеих частей уравнения на ао
. Обозначая Ai
= ai
/ao
, Wi
= hi
-1
/ho
-1
и V = ho
*/(ho
-1
ao
) = ho
ho
*/ao
, получаем:
(13.3.11)
где для значений W и V введен индекс j номеров предстоящих итераций по циклу вычисления коэффициентов фильтра.
При нулевой итерации (N=0, j=0) имеем только одно уравнение:
. (13.3.12)
Благодаря проведенной нормировке решения уравнения (13.3.12) не требуется:
А0
= 1, V0
= 1, W00
= 1.
Составим уравнение для двучленного фильтра (N=1, j=1):
(13.3.13)
Перепишем уравнение (13.3.12) в прямой форме:
А0
W00
= V0
. (13.3.14)
Запишем вспомогательную систему, для чего к уравнению (13.3.14) добавим вторую строку с новой постоянной Ej
:
A0
W00
+ A1
·0 = V0,
A1
W00
+ A0
·0 = E1
.
В матричной форме:
(13.3.15)
Реверсируем уравнение (13.3.15):
(13.3.16)
Вычтем (13.3.16) из (13.3.15) с неопределенным множителем Rj
:
(13.3.17)
Из верхней строки уравнения (13.3.16) можно получить значение Е1
:
Е1
= A1
W00
. (13.3.18)
Уравнение (13.3.13) можно сделать равнозначным уравнению (13.3.17), если правую часть нижней строки уравнения (13.3.17) приравнять к правой части нижней строки уравнения (13.3.13):
E1
- R1
V0
= 0, R1
= E1
/V0
. (13.3.19)
При этом из правых частей верхних строк уравнений (13.3.13, 13.3.17) будем иметь:
V1
= V0
- R1
E1
. (13.3.20)
Приравнивая друг другу левые части уравнений (13.3.13,13.3.17), получаем:
W01
= W00
- R1
·0 = W00
= 1.
W11
= 0 - R1
W00
= -R1
W00
. (13.3.21)
Этим заканчивается первая итерация. Аналогично, для второй итерации:
(13.3.22)
(13.3.23)
(13.3.24)
(13.3.25)
Из верхней строки уравнения (13.3.24):
Е2
= A1
W11
+A2
W01
.
Из правых частей нижней и верхней строк уравнений (13.3.22, 13.3.25):
R2
= E2
/V1
,
V2
= V1
- R2
E2
.
Новые коэффициенты из левых частей уравнений (13.3.22, 13.3.25):
W02
= W01
- R2
0= 1,
W12
= W11
- R2
W11
,
W22
= 0 - R2
W01
.
Анализ расчетов позволяет вывести следующие рекуррентные формулы:
Ej
=Ai
Wj-i
,j
, j = 1,2,...,M. (13.3.26)
Rj
= Ej
/Vj-1
,
Vj
= Vj-1
- Rj
Ej
,
Wi,j
= Wi,j-1
- Rj
Wj-1
,j-1
, i = 0,1,.., j.
Подпрограммы решения уравнений для ЭВМ приведены в литературе /12,22/.
Курсовая работа 21-07.
Разработка программы расчета оптимального фильтра деконволюции сигналов и фильтрации цифровых данных.
13.4. Рекурсивная деконволюция /22/.
Уравнение фильтра рекурсивной деконволюции
. Запишем уравнение (13.1.2) для инверсного фильтра в развернутой форме:
H-1
(z) = 1/(h0
+h1
z+h2
z2
+ ...). (13.4.1)
Так как для минимально-фазового оператора всегда выполняется условие h0
0, приведем (13.4.1) к виду:
H-1
(z) = (1/h0
)/(1+h1
z/h0
+h2
z2
/h0
+...) = G/(1+q1
z+q2
z2
+ ...), (13.4.2)
где: G = 1/h0
, q1
= h1
/h0
, q2
= h2
/h0
и т.д. Но уравнение (13.4.2) есть не что иное, как уравнение передаточной функции рекурсивного фильтра, где цепь обратной связи образована коэффициентами нормированного оператора h(n). Алгоритм вычислений:
yk
= G·xk
–qn
·yk-n
.
Выражение (13.4.2) уникально по своим возможностям. В принципе, оно может реализовать оператор инверсной фильтрации с бесконечным импульсным откликом. На практике оно может использоваться вместо медленно затухающих инверсных операторов, модуль одного из полюсов которого очень близок к 1 (менее 1.1) при высоких требованиях к метрике приближения.
Пример рекурсивной деконволюции.
|
Рис.13.4.1. |
Оператор hn
= {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. Форма операторов приведена на рис. 13.4.1.
|
Рис.13.4.2. |
2. При использовании оптимального инверсного оператора с N=100 значение погрешности приближения уменьшается более чем в 20 раз, что позволяет уменьшить длину оператора до N=35 при погрешности приближения порядка 0.1 (рис 13.4.2(А)), при этом абсолютные значения погрешностей приближения не превышают 0.03 (рис. 13.4.2(В)).
3. Расчет коэффициентов фильтра рекурсивной деконволюции:
G = 1/ho
= 2.441
gn
= hn
·G. gn
= {1.932, 0.978, -0.472, -0.896, -0.405, 0.077, 0.165, 0.065, -0.003}, n=1, …, 9.
|
Рис. 13.4.3. |
На рис. 13.4.3 приведен результат рекурсивной деконволюции оператора hn
. Деконволюция абсолютно точно, с нулевой метрикой, восстанавливает импульс Кронекера, хотя собственный импульсный отклик рекурсивного оператора повторяет оператор h-1
n
при его вычислении по формуле (13.2.2) и длительность его значимой части близка к 200. Коэффициент усиления дисперсии шумов при данной операции вычисляется по значениям импульсного отклика оператора рекурсивной деконволюции и весьма существенен, как и для всех инверсных операторов.
Курсовая работа 22-07.
Разработка программы расчета рекурсивного фильтра деконволюции сигналов и фильтрации цифровых данных.
13.5. Фильтры СЖАТИЯ СИГНАЛОВ
Рассмотренные выше методы относятся к расчету так называемых идеальных инверсных фильтров, т.е. фильтров полной инверсии системных операторов. Однако использование идеальных инверсных фильтров на практике не всегда возможно, т.к. регистрируемые данные обычно осложнены влиянием помех (шумов), а инверсные фильтры обычно имеют коэффициент усиления дисперсии шумов значительно больше 1. В этом случае задача точной деконволюции (восстановления истинной первоначальной формы сигнала), как правило, не ставится, а инверсные фильтры считаются оптимальными с точки зрения максимального приближения к форме полезного сигнала с определенным допустимым коэффициентом усиления дисперсии помех. Такие фильтры называются фильтрами неполной (частичной, ограниченной) деконволюции или фильтрами сжатия сигнала. При проектировании фильтров неполной деконволюции учитываются статистические характеристики помех во входном сигнале и их соотношение со статистическими характеристиками самого входного сигнала.
Передаточная функция фильтра
неполной деконволюции с учетом помех во входном сигнале определяется выражением:
H-1
(z) = H*(z)/[|H(z)|2
+g2
], (13.5.1)
где g2
= k·sh
2
- дисперсия шумов в единицах дисперсии оператора hn
, sh
2
– дисперсия значений оператора hn
, (при условии суммы значений оператора, равной 1), k - отношение дисперсии шумов к дисперсии оператора hn
. Коэффициент g2
играет роль регуляризирующего фактора при выполнении операции деконволюции информации.
|
Рис. 13.5.1. |
На рис. 13.5.1 пример формы оператора hn
и спектральных функций (13.5.1) при разных значениях параметра g. При g = 0 выражение (13.5.1) обращается в идеальный инверсный фильтр 1/H(z). Во втором крайнем случае, при g2
>>|H(z)|2
, фильтр (13.5.1) переходит в фильтр, согласованный с сигналом по частотному спектру: H-1
(z) = H*(z)/g2
, который только максимизирует отношение сигнал/помеха.
На рис. 13.5.2 приведена форма инверсных операторов, соответствующая их частотным характеристикам на рис. 13.5.1(В), и результаты свертки инверсных операторов с прямым (для лучшего просмотра графики прямой оператор при свертке сдвинут вправо на 2 значения Dt). При g=0 коэффициент усиления дисперсии шумов равен 11, при g=0.4sh
2
равен 4.6. Однако снижение усиления дисперсии шумов сопровождается увеличением погрешности приближения, что можно видеть на рис. 13.5.2(В), при этом уменьшается амплитуда восстановления импульса Кронекера и появляются осцилляции после импульса. Но при наличии шумов и правильном выборе параметра g общее отношение амплитудных значений сигнал/ шум для оператора по (13.5.1) больше, чем для прямой инверсии по (13.1.2), что объясняется более существенным уменьшением коэффициента усиления дисперсии шумов при увеличении параметра g, чем увеличением погрешности приближения.
|
Рис. 13.5.2. |
Операторы оптимальных фильтров
сжатия сигналов также могут вычисляться с учетом помех. Если сигнал s(k) и помеха статистически независимы, то функция автоковариации сигнала на входе фильтра:
ai
= asi
+ bi
, (13.5.2)
где asi
и bi
- функции автоковариации сигнала и помех. При помехе типа белого шума функция автоковариации помех представляет собой весовую дельта-функцию в точке 0:
bi
= c2
di
, (13.5.3)
где с2
- дисперсия помех. С учетом этого фактора расчет оптимальных инверсных фильтров может проводиться по вышеприведенным формулам (13.3.5, 13.3.9) с изменением значения коэффициента ао
:
ao
= ao
+ c2
. (13.5.4)
На рис. 13.5.3(А) приведены примеры операторов оптимальных инверсных фильтров, вычисленные по прямому оператору, приведенному на рис. 13.5.1(А). Значения коэффициента с2
заданы в долях дисперсии прямого оператора. Ввод коэффициента с2
в функцию автоковариации резко уменьшает значения коэффициентов инверсного оператора и, соответственно, уменьшает коэффициент усиления дисперсии помех.
Рис. 13.5.3.
Для приведенного примера, при исходном значении коэффициента усиления дисперсии шумов порядка 12 для с2
=0, его значение уменьшается до 1.8 при с2
=0.1s2
и становится меньше 1 при с2
> 0.3s2
. Естественно, что общая погрешность приближения деконволюции при этом также существенно изменяется (см. рис. 13.5.3(В)), но амплитуда значения сигнала на месте импульса Кронекера (там, где он должен быть) изменяется много меньше, чем коэффициент усиления дисперсии шумов, а, следовательно, отношение сигнал/шум при введении коэффициента с2
существенно увеличивается.
литература
12. Канасевич Э.Р. Анализ временных последовательностей в геофизике. - М.: Недра, 1985.- 300 с.
22. Рапопорт М.Б. Вычислительная техника в полевой геофизике: Учебник для вузов. - М.: Недра, 1993.- 350 с.
А.В.Давыдов.
20.08.07.
Cайт автора ¨ Лекции ¨ Практикум
О замеченных опечатках, ошибках и предложениях по дополнению: davpro@yandex.ru.
Copyright ©2008 Davydov А.V.