Модифицированные вейвлеты в обработке данных аналитических приборов. Ii. Алгоритмы обработки Текст научной статьи по специальности «Компьютерные и информационные науки»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2006, том 16, № 2, c. 78-91

ОРИГИНАЛЬНЫЕ СТАТЬИ

УДК 621.391.26 © Л. В. Новиков

МОДИФИЦИРОВАННЫЕ ВЕЙВЛЕТЫ В ОБРАБОТКЕ ДАННЫХ АНАЛИТИЧЕСКИХ ПРИБОРОВ. II. АЛГОРИТМЫ ОБРАБОТКИ

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

ВВЕДЕНИЕ

В работе [1] рассмотрены вейвлеты, полученные модификацией известных ортонормирован-ных вейвлетных базисов (вейвлетов-"праро-

дителей") {(pjk (t),Wj,k (t); je Z,ke N} некоторой функцией G (t). Новые базисные функции, обозначаемые как Qjk (t),дJk (t) — масштабирующие функции, T]jk (t),ñj k (t) — вейвлетные функции, позволяют, в частности, вычислить оценку полезного сигнала x (t), искаженного аппаратной функцией H (t) и шумом u (t) в наблюдаемом сигнале y (t)

y(t) = ¡H(t-t) x(r)dr + u(t). (1)

Эта оценка имеет вид

x(t) = E^ (k) j (t)+E Edj (k)Wjk (t), (2а)

k j =jo k

где j0 и J — самое тонкое и самое грубое значения масштаба соответственно. Показано, что получение оценки по формуле (2а) эквивалентно обработке наблюдаемого сигнала фильтром с частотным откликом G (ю)

x (a) = G (ю)у (a)

С0 (k)=Jу (t(t)dt, (k )=J y (t )ñjk (t )d^

(за)

(зб)

где функции в0.к () и Ц0.к () должны удовлетворять следующим масштабирующим уравнениям:

или

я, ,k (t )=Еh (n -2k )§j+hn (),

n

n, ,k (t )=E g (n -2k )§,+1,n (t),

n

§j (tЬЕ^ +1 (n)<Pj+1,n (t) , n

n j(t )=^Zß,+i(n M+i,n(t),

(4а) (4б)

(5а) (5б)

где к и g — масштабирующие и вейвлетные коэффициенты вейвлетов-"прародителей", а и /3 — масштабирующие и вейвлетные коэффициенты новых вейвлетов, которые определяются по формулам:

а

j+1

(n ) =

, 2 J+ln

1— J G(a)h (2-J-1rn)e'2-J-'anda, (6а)

2 j+1 2n

-2 jt1l

(2б)

ßj+1 (n ) =

и вейвлетным подавлением шума по алгоритму, предложенному в работе [3].

Коэффициенты вейвлет-разложения наблюдаемого сигнала у () — '¿^ (к) и (к) определяются как

, 2 J+ln

1— J G (a)g (2-J-1a)e' 2-J-'anda. (6б)

2 j+1 2п

-2 П

Обозначим

ä, +1 (ю)= G (a)h (2-j-1a),

(7а)

, (ю) = О (<о)ё (2--1т). (7б)

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

ложения с, (к) и ё, (к). При удачном выборе

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

1. БЫСТРЫЕ ВЫЧИСЛИТЕЛЬНЫЕ АЛГОРИТМЫ И БЛОКИ ФИЛЬТРОВ

Анализ

Чтобы получить быстрые алгоритмы для вычисления коэффициентов вейвлет-разложения

оценки полезного сигнала х (), подставим в (3 а, 3б) вместо в, к и г\,к их выражения из (4а) и (4б) соответственно. Получим:

С, (к) = Ек(т - 2к)С,+1 (к), (8а)

т

(к) = Е 8 (т - 2к), (к). (8б)

т

Так же как и в случае традиционного КМА, мы получили рекурсивный алгоритм для вычисления коэффициентов, который требует знания начальных значений с} (к) и (к) на самом "тонком" масштабе J:

С] (к) = /у (в (0^, (9а)

dJ (к) = /у (()dt. (9б)

Чтобы вычислить эти коэффициенты, так же как и в случае традиционного КМА, необходимо воспользоваться отсчетами наблюдаемого сигнала

у (). Подставляя в (9а, 9б) выражения для вв,к и П,к из (5а, 5б) и выполняя преобразования, получим

С] (к) = +1 (т - 2к)/у (+1,т ()dt =

т

= Ху (т)й]+1 (т - 2к) (10а)

т

и

ё] (к) = £у (т)&+1 (т - 2к). (10б)

Рис. 1. Вейвлетное подавление шума. ДВП — дискретное вейвлетное преобразование; ПО — пороговая обработка; ОДВП — обратное ДВП; ОП — оценка параметров; у ^) — наблюдаемый сигнал;

с, (к), dj (к) — коэффициенты вейвлет-разложения наблюдаемого сигнала; с, (к), ё, (к) — коэффициенты вейвлет-разложения наблюдаемого сигнала после пороговой обработки; х ) — оценка полезного сигнала

(X,

(-п)

у(п)

с/-1

к(-п) к(-п)

> С/-3

&+1(-п) \ 8 (-п)

(-п)

к (-п)

С,0 +1 -►

8(-п)

к (п)

(п)

'(п)

к(п) к(п)

:(п) (п)

к (п)

С,0 +1

б

Рис. 2. Схема 4-ступенчатого анализирующего (а) и синтезирующего (б) блоков фильтров.

а+1, к

и в/+1, 8 — коэффициенты низкочастотного и высокочастотного фильтров вейвлетов

► -2

/-1

/-2

/-3

а

/-3

/-1

/-2

/-3

/-2

/ -1

По формулам (8*) и (10*) можно построить вычислительную схему алгоритма анализа наблюдаемых данных, которая показана на рис. 2, а. В отличие от схемы [1, рис. 3, а] отсчеты наблюдаемого сигнала в первом каскаде подвергаются обработке низкочастотным и высокочастотным фильтрами с другими коэффициентами, а именно +1 и ¡3/+1 соответственно. Коэффициенты фильтров на последующих ступенях анализа остаются прежними.

Синтез

Учитывая что О,+1 = О, ©М, (см. [1, п. 3]), базисную функцию вв,+1 к е и,+1 можно представить в виде

в, +1, т ^ ) =

= Ш+1твк()+1(в,+1т,,к ()

к т

или с учетом формул (26) и (29) в [1]:

в +1,т ^ ) =

= Xк (т - 2к)вк () + X8 (т - 2кл () (11)

к т

Умножая далее (11) скалярно слева и справа на у (), получим

+1(т ) =

= Xк (т - 2к)с, (к) + X8 (т - 2к(к). (12)

к т

Мы получили, таким образом, рекурсивную формулу восстановления коэффициентов вейвлет-разложения искомой оценки х () "от грубого к тонкому разрешению" с использованием, как и в анализе, фильтров-"прародителей" к и 8 . Вычислительная схема восстановления отсчетов оценки х () ничем не отличается от схемы традиционного КМА (ср. [1, рис. 3, б] и рис. 2, б).

Покажем теперь, что схема анализа—синтеза, показанная на рис 2, обеспечивает получение оценки полезного сигнала х () по наблюдаемому сигналу у^), если значение масштаба / выбрано равным -1, т. е. при / = -1. Чтобы положить коэффициенты с0 (к), равными отсчетам наблюдаемого сигнала у (к), т. е. с0 (к) = у (к), как уже отмечалось в [1], спектр масштабирующей функции ф0 (ю) в пределах спектральной полосы у (о)

должен быть максимально плоским. Это означает, что при использовании в качестве "прародителей" вейвлетов Добеши, имеющих плоскую частотную характеристику вплоть до частоты Найквиста, частота отсчетов наблюдаемого сигнала должна быть меньше п . Тогда из формулы (10) на первом шаге анализа необходимо выполнить:

С-1 (к) = ^<х0 (т - 2к)у (т),

т

ё-1 (к) = ^А) (т - 2к)у (т).

(12а) (12б)

■(ю)=

1

а0 (ю)к (ю)+ ¡30 (ю)к (ю)

1

+— 2

а0 (ю + п)к(ю)+ ¡30 (ю + п)к (о)

у (ю) + у (ю+п).

Вторая сумма в квадратных скобках равна нулю, и, следовательно, равна нулю "элайзинговая" составляющая синтезируемого сигнала. Действительно, имеем с учетом (7*) при 0 = -1

а0 (ю + п)к (ю) + /30 (ю + п)к (ю) = = О (ю + п)к (ю + п)к (ю) + +О(ю + п)к (ю + п)к(ю) =

= О (ю + п)х

X

к (ю)к (ю + п) + ехр(т)к (ю)к (ю + п) Первая сумма в квадратных скобках равна

= 0.

а0 (ю)к (ю)+ ¡30 (ю)к (ю) =

Рассмотрим работу первого каскада анализа и последнего каскада синтеза блока фильтров, который должен восстанавливать искомый сигнал в его дискретных точках х (к) (рис. 3). Можно показать [2], что вход и выход такого фильтра связаны выражением

= О(ю)к(ю)к(ю) + О(ю)к (ю)к (ю) = = О(ю) \к(ю)|2 + \к(ю) = 20(ю).

Следовательно,

х (ю) = О (ю)у (ю).

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

2. ОПИСАНИЕ АЛГОРИТМОВ

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

а) выбор или синтез функции О () в зависимости от круга задач решаемых прибором;

в (-П) -► |2 а_1 -► 12 -► g(n)

Рис. 3. Первая ступень блока фильтров анализа-синтеза: у(т)— наблюдаемый сигнал; х(т) — оценка сигнала; |2 и |2 — прореживание и интерполяция с коэффициентом 2; с-1 и ё-1 — коэффициенты анализа сигнала первой ступени

б) выбор вейвлета-"прародителя" (коэффициен- Интегрирование. При оценке текущего значе-тов к (п) и 8 (п)); t

ния интеграла J х (т^т с учетом искажений,

в) вычисление коэффициентов а0 (п) и (30 (п); 1-Т

г) ДВП; вносимых прибором и шумами, функция О(ю)

д) пороговая обработка вейвлет-коэффициен- будет иметь вид

тов;

е) ОДВП; О(ю) = ( ю)-1 _ Н (ю)-. (17)

ж) оценивание параметров пиков. |Н (ю) | + в Я (ю)

Рассмотрим более подробно алгоритмы, реализуемые на каждом этапе обработки. (Карунена—Лоэва)-подобные разложения

(декорреляция). Известно, что разложение Кару. . нена—Лоэва приводит к некоррелированным ко-Выбор функции О () эффициентам разложения [6]. Однако этот же результат можно получить с помощью предлагаемых Деконволюция. Как уже отмечалось, выбор вейвлетов. Предположим, что на входе системы

функции О(ю), равной Н-1 (ю), приводит к Н(ю) действует белый шум с единичной диспер-

большому уровню выходного шума. С целью по- сией х ^). Тогда на ее выходе он оказывается кор-

давления эффекта бесконечного усиления на прак- релированным с корреляционной функцией тике применяют окна Ж (ю) или используют восстанавливающий фильтр Винера—Тихонова [4], к(т) = -^ [К(ю)ехр( ЮT)dT , регуляризирующий усиление на верхних частотах, 2п ^

Н (ю) где К (ю) — спектральная плотность мощности и

О (ю) = --^-, (14) „ .2

|Н (ю) | +вЯ (ю) К (ю) = Н (ю) .

где в — регуляризирующий параметр, Я(ю) — Выберем О(ю) = К(ю)-1/2 и покажем, что ко-

некоторый четный полином. , , _ ^

, ч эффициенты разложения с, (к) и (к) при этом

Выбирая константу р и полином Я(ю), мож- л 1

будут некоррелированными. Действительно, с

но легк°з корректировать ход частотной характери- учетом (3а) получим

стики

О (ю).

Е[с, (т)с, (к)] =

Дифференцирование. Чтобы получить оценку

р-й производной сигнала х(), искаженного при- =JJЕ[у^)у(т)^,^ ()в,к (т^dт =

бором и шумами, функция О(ю) будет иметь вид =JJк(-т)§ (t)§, (т)dtdт =

О (ю) = ( ю)р-Н(ю)-. (15) = ^К(ю)(ю) ]\ехР((ю(т--ю =

Н (ю)| (ю) I . , , ч

= — I ^(ю)) ехр(ю(т-п))dю = 8(m-п),

Если же необходимо получить оценку произ- 2п •м 1

водной полезного сигнала без учета аппаратной

, 11 где символ Е обозначает статистическое усредне-

функции, то можно воспользоваться дифференци- ние

рующими фильтрами Савицкого—Голэя [5]. Тогда НИ Аналогично можно показать, что

функция О(ю) будет частотным откликом цифро- _ _

вого дифференцирующего фильтра, т. е. Е(т)а, (п)] = $(т - п).

Из формул [1, (19а, 21а)] при к = 0 можно получить интересный результат:

О (ю) = О (ехр ( ю)) = О0 ^), (16)

где 2 = ехр (/ю). К (юв {ю) = 0; (ю),

или

IК((т)^т = в] ()

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

Однако так же, как в случае с деконволюцией,

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

Коррекция частотной характеристики системы. Чтобы скорректировать частотную характеристику системы, необходимо умножить функцию О(ю) на подходящий полином от частоты ю . Например, чтобы поднять верхние частоты, достаточно умножить эту функцию на полином вида 1 + а0 ю2.

Выбор модифицируемых вейвлетов. Как было показано в [1], функции в,в,ц,ц порождают КМА, если в качестве модифицируемого выбран вейвлет с компактным носителем в частотной области (вейвлет Мейера или вейвлет Шеннона). Однако эти вейвлеты во временной области затухают очень медленно, определяются большим числом коэффициентов к и g и неудобны для приложений. Поэтому в качестве модифицируемых вейвлетов необходимо использовать другие вейвлеты, например Добеши высоких порядков, имеющих хорошую крутизну фронтов частотного спектра. Возникающая при этом систематическая ошибка может быть оценена численно в каждом конкретном случае.

Вычисление коэффициентов а0 и ¡в0. Эти коэффициенты определяются формулами (6а) и (6б) при 0 = -1:

1 п

а (п) =-1 О(ю)к(ю)еюМю,

2п -п

А (п ) = -2^1 О (юуё (ю)ею Мю.

(18а)

(18б)

эффициенты а0 и во вычисляются по формуле обратного дискретного преобразования Фурье (ДПФ), которая для а0 имеет вид

осп

1 N -1

(п) = ^ Е О (т)к (т)егп

1 т=-Ы1

для всех п = (М0: М1), где М0 и М1 выбираются такими, чтобы а0 (п) ~ 0 при М0 > п > М1.

По аналогичной формуле вычисляются коэффициенты в0 •

Если приборная функция Н (), которая в ряде задач определяет О (), задана на дискретной сетке t = (-N0 : Ы0), то ее ДПФ вычисляется по формуле

Н (т ) = Е Н (к )е

-гп mk/N1

к=0

для всех т = (-N1 : N -1).

Если же функция О(ю) определена по формуле (16) и так как к (п) и g (п) являются КИХ-

фильтрами, коэффициенты а0 (п) и (30 (п) могут быть получены путем умножения полиномов по

степеням г:

«0 (г) = (( * к)(г), Д, (г) = ( * к ))г).

(19а) (19б)

Для вычислений по этим формулам функции О (ю), к (ю) и g (ю) должны быть определены на дискретном множестве (сетке) частоты ю

ю = (N1 : N1 - 1)п/ N1,

где N выбирается в пределах 100^200, что дает вполне удовлетворительный результат. Тогда ко-

Тогда искомые коэффициенты будут коэффициентами полиномов а0 (г) и (30 (г).

Быстрые дискретные вейвлетные преобразования (ДВП). ДВП имеют целью вычисление коэффициентов с ■ (к) и (к) для значений масштабов 0 > 0 > 00. Вычислительные процедуры и реализующие их блоки фильтров рассмотрены в [1, раздел 1] и здесь в разделе 1. Из формул (8*) и (10*) следует, что в основе структуры ДВП лежит блок фильтров, показанный в [1, рис. 4]. Их последовательное соединение по схеме рис. 2, а реализует процедуру анализа наблюдаемого сигнала у ^). В отличие от традиционного КМА, как показано на рис. 2, а, для достижения требуемого результата первый каскад блока фильтров должен быть снабжен коэффициентами фильтров а0 (п) и

в0 (п). Напомним, что операция децимации на каждом шаге анализа понижает число коэффициентов вдвое. Поэтому, чтобы достигнуть наиболь-

шей глубины анализа по масштабу, выбирают число отсчетов наблюдаемого сигнала y (n) кратным степени двойки, т. е. N = 2J°. Тогда формально можно выполнить максимально J0 шагов анализа (J0 = log2 (N)), получив в результате N -1 коэффициентов d. и один коэффициент cjo .

На практике часто применяют ДВП с избыточным числом коэффициентов (избыточное (redundant) ДВП), при котором иногда достигается лучшее подавление шума. В этом преобразовании в результате децимации коэффициенты не отбрасываются, а сохраняются. Тогда на первом шаге анализа получим пару векторов dj (к): dj (к) и dj (к) длиной N/2 каждый — и пару векторов с. (к): cj (к) и cB (к) длиной N/2 каждый. Векторы cj (к) и cB (к) на втором шаге подвергаются аналогичной обработке, в результате которой получаем четыре вектора dj (к) длиной N/4 каждый и

четыре вектора с. (к) длиной N/4 каждый, и т. д.; выполнив J0 шагов, получим N log2 (N) коэффициентов d. и с ■ .

J J0

Программы, реализующие ДВП и избыточное ДВП, содержатся в последних версиях МАТЛАБ, а также их можно найти на различных сайтах ИНТЕРНЕТ, например [7].

Пороговая обработка коэффициентов вейв-лет-разложения. Пороговая обработка имеет целью удаление шума из наблюдаемых данных. Процедура пороговой обработки осуществляется только для коэффициентов dj (к), которые сравниваются с порогом Tr, обнуляются, если этот коэффициент ниже порога, или сохраняют свое значение в противоположном случае. Очевидно, что идеальное удаление шума будет в том случае, если энергия сигнала в результате ДВП оказывается сосредоточенной в минимальном числе коэффициентов, а энергия шума размыта по максимальному их числу. Этого можно достичь в том случае, если базис будет хорошо коррелирован с сигналом, что и происходит для сигналов, имеющих небольшую продолжительность (для пиков), при применении вейвлетного базиса. Очевидно, что наибольшая корреляция будет в том случае, если масштабирующие и вейвлетные функции совпадают с формой сигнала.

Различают жесткую и мягкую пороговые обработки. Жесткая пороговая обработка строится по алгоритму [8]

dj (к ) =

У. (к), если У (к)|> Тг; 0, если (к )|< Тг,. а мягкая пороговая обработка:

У. (к)- Тг.., если (к)|> Тг,. У. (к) + Тг., если

dj (к ) =

0,

если

(к )| (к ) ,(к )

< -ТГ..

< ТГ..

где порог для каждого значения масштаба . определяется как

Tr. = O. .J2ln N . (20)

Дисперсия шума о2 определяется по формуле а2 =о2 Xd0,j (n) , (21)

где — коэффициенты разложения функции

О (^) по вейвлетным функциям У0. = = |О ()ц.к () dt, а2 — дисперсия шума и ().

Эта дисперсия оценивается по наблюдаемым данным у (^) путем вычисления медианы коэффициентов У. (к) на самой тонкой шкале (. = 0). Если эта медиана равна МУ, то оценка а будет

о =-

0.6745

(22)

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

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

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

3. ПРИМЕРЫ ПРИМЕНЕНИЙ

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

| х (т)а

т при условии минимизации величины

соответствующего среднеквадратического отклонения (СКО):

х (г)-х (г )|<е0,

х (р)(г)- х (г )(р)

I I

| х(т)ёт- | х(т)ёт

< ет

где £, > 0 — некоторая заданная ошибка.

Чаще всего используют относительное СКО, определяемое как

1

ИМ8Б(х - х) = {е2/Лх (г)2 . (23)

1.4 1.2 1

0.8 0.6 0.4 0.2 0.

х()

1.

1.2 1

0.8 0.6 0.4 0.2-o^ -0.2,

у(г)

0

100 200 300 400 500 600

0

100 200 300 400 500 600

1.2 1

0.8 0.6 0.4 0.2 0 -0.2,

х(г (

0

100 200 300 400 500 600

1.2 1

0.8 0.6 0.4 0.2 0 -0.2

х(г( к\ \

0 100 200 300 400 500 600

Рис. 4. Пример деконволюции модельного сигнала.

а — оригинальный сигнал х (г); б — наблюдаемый сигнал у (г) с шумом при отношении сигнал/шум 10 (ИМ8Б (х - у ) = 0.1); в — оценка сигнала х (г) с шумом (ИМ8Б (х - х ) = 0.03, в = 0.02, Я(ю) = 1); г — то же, что и (в), но с усилением верхних частот с помощью полинома

1 + со

б

а

г

г

в

г

г

г

800

600

400

200

-200.

200 400 600 800 1000 1200

2000

1500

1000

500

-500,

X

()

200 400 600 800 1000 1200

Рис. 5. Фрагмент спектра инсулина времяпролетного масс-спектрометра. а — наблюдаемый сигнал, б — оценка сигнала, полученная методом обратной свертки

б

а

t

t

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

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

заданной формулой (14). Выберем модельный сигнал с белым шумом, определенный на множе-

стве точек N = 2

состоящий из шести пиков га-

( и *\2 \

уссовой формы ак ехр

( - Ч)

2^с2

где

ак =[0.1 0.251.0 0.71.0 0.35], tk = [100 138 150159 280 290], /л0 = 3.2.

Приборная функция t

Н () = ехр

2^0

принята равной Я(ю) = 1, и коэффициенты

а0 (п) и во (п) определялись по формулам (18*).

Результаты приведены на рис. 4. Как следует из рис. 4, г, определенный положительный эффект

дает усиление верхних частот умножением О(ю)

на полином 1 + ю 2 .

На рис. 5-8 приведены примеры обработки реальных сигналов времяпролетного масс-спектрометра, спектрометра Мессбауэра, электронного спектрометра и секвенатора НАНОФОР 3. Приборные функции в первых двух случаях определялись по форме одиночного пика этих приборов, во втором и третьем случаях выбиралась в виде гауссовой кривой. Фильтр О (ю) строился аналогично предыдущему примеру.

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

Сравним этот метод с существующим, например использующим дифференцирующие фильтры второго порядка Савицкого—Голэя (СГ) [5], на модельном сигнале рис. 4, а (см. рис. 9).

50 100 150 200 250 300

50 100 150 200 250 300

Рис. 6. Мессбауэровский спектр продуктов реакции образования 8и8 (олово—сера). а — наблюдаемый сигнал; б — оценка сигнала, полученная методом обратной свертки с усилением верхних частот с помощью полинома 1 + 10ю2

250

250

0 100 200 300 400 500 600

-50.

0 100 200 300 400 500 600

Рис. 7. Фотоэлектронный спектр.

а — наблюдаемый сигнал; б — оценка сигнала, полученная методом обратной свертки с усилением верхних частот с помощью полинома 1 + 150ю2

б

а

г

t

б

а

г

г

Ошибка дифференцирования определялась по величине СКО от истинного значения производной х", вычисленной по известным правилам дифференцирования. Если при отсутствии шума (рис. 9, б) этой ошибкой можно пренебречь (КМ8Е (х" - х ) ~ = 10-12), то в присутствии шума (рис. 9, в) метод СГ приводит к величине СКО ~28 %. Более того, производная от пика минимальной интенсивности (первый пик) оказалась полностью скрытой в шумах. Предлагаемый метод (рис. 9, г) приводит к величине СКО, равной ~8 %, и к четкому выделению производной минимального пика на фоне шумов.

Вычисление коэффициентов а0 (п) и в0 (п)

в этом случае производилось по формулам (19*).

(Карунена—Лоэва)-подобные разложения (декорреляция). В формуле (14) положим Я(ю) =

= 1 и в = 3.5 • 10-8. В качестве модельного сигнала х (г) выберем белый гауссовский шум, N = 210 и фильтр

Н (ю) = 1Г^=7

V! + У ю2

при у2 = 10 .

Тогда случайный процесс на выходе фильтра

будет описываться экспоненциальной корреляцион-

ной функцией К (т) = ехр < -

Коэффициенты

а0 (п) и в0 (п) вычисляются по формулам (19*). Результаты обработки приведены на рис. 10.

а

2600 2800 3000 3200 Время задержки, с

3400 3600

2800 3000 3200 Время задержки, с

б

3600

Рис. 8. Фрагмент сигнала ДНК, полученного на приборе НАНОФОР 3. а — наблюдаемый сигнал; б — оценка сигнала, полученная методом обратной свертки

1.2 1 0.8 0.6 0.4 0.2 0 -0.2.

у ()

0 100 200 300 400 500 600

а

0.04 0.02 0

-0.02 -0.04 -0.06 -0.08.

хТ (()|

г »

;

0 100 200 300 400 500 600 б

0.04 0.02 0

-0.02 -0.04 -0.06 -0.08,

х\1)

О 100 200 300 400 500 600

в

0.04г 0.020-

-0.08.

х" (() 1

I

О 100 200 300 400 500 600

г

Рис. 9. Пример вычисления второй производной модельного сигнала с использованием девятиточечного третьего порядка фильтра Савицкого—Голэя.

а — наблюдаемый сигнал у ) с шумом при отношении сигнал/шум 5 (ИМ8Е (х - у ) = 0.1); б — оценка второй производной без шума (КМ8Е(х"-X") = 2-10-12, в = 0.004, Я(ю) = 1); в — оценка второй производной без подавления шума (КМ8Е(х"-X") = 0.28, в = 0.04, Я(со) = 1); г — оценка второй производной с подавлением шума (ИМ8Е(х"-X") = 0.087, в = 0.04, Я(ю) = 1)

0.8

0.6

0.4

0.2

1 H (т) ■

т/п

0.2 pi 0.6 0.8 в

Рис. 10. Декорреляция (отбеливание) шума.

а — фрагмент тест-сигнала (N=1024); б — наблюдаемый сигнал у (г), в — частотная характеристика системы Н (ю ) =-1-

v 7

; г — декоррелированный сигнал (RMSE(x-x/)=0.028)

СПИСОК ЛИТЕРАТУРЫ

1. Новиков Л.В. Модифицированные вейвлеты в обработке данных аналитических приборов. I. Основы теории // Научное приборостроение. 2006. Т. 16, № 1. С. 3-15.

2. Добеши И. Десять лекций по вейвлетам. М.: РХД, 2002. 464 с.

3. Donoho D., Johnstoune I. Ideal spatial adaptation via wavelet shrinkage // Biometrika. 1994. V. 81. P.425-455.

4. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979. 288 с.

5. Savitzky A., Golay M.J.E. Smoothing and differentiation of data by simplified least squares procedures // Analytical Chemistry. 1964. V. 36,

N 8. P.1627-1639.

6. Френкс Л.Е. Теория сигналов. М.: Советское радио, 1974. 343 с.

7. (www-dsp.rice.edu/software/rwt.shtml); (www-stat.stanford.edu/~wavelab).

8. Mallat S. A wavelet tour of signal processing. San Diego: Academic Press, 1999. 637 p.

Институт аналитического приборостроения РАН, Санкт-Петербург

Материал поступил в редакцию 7.07.2005.

MODIFIED WAVELETS IN DATA PROCESSING FOR ANALYTICAL INSTRUMENTS.

II. DATA PROCESSING ALGORITHMS

L. V. Novikov

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The paper describes data processing algorithms using modified wavelets. The fast computational algorithms are implemented with filter banks. The application is illustrated by the examples of signal processing for various analytical instruments. It is shown that the proposed approach combines wavelet based noise rejection and de-convolution or other signal transformations, provides several times higher sensitivity and resolution of the instrument.