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

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2006, том 16, № 1, c. 3-15

ОБЗОРЫ= =

УДК 621. 391. 26 © Л. В. Новиков

МОДИФИЦИРОВАННЫЕ ВЕЙВЛЕТЫ В ОБРАБОТКЕ ДАННЫХ

АНАЛИТИЧЕСКИХ ПРИБОРОВ. I. ОСНОВЫ ТЕОРИИ

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

ВВЕДЕНИЕ

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

Существующие подходы решения этих задач включают применение различных методов сглаживания (с помощью фильтров Савицкого—Го-лэя, Калмана, Винера, Фурье и др.), дифференцирующих фильтров и методов деконволюции [1]. В последнее десятилетие в практику обработки данных аналитических приборов стали внедряться подходы, основанные на время-масштабном преобразовании сигнала в вейвлетном базисе [2], включающем масштабирующие qjk (t ) и вейвлет-

ные функции yjk (t) ( jе Z определяет масштаб

(ширину), а k е N — смещение этих функций).

Уже опубликовано более 370 статей по применению вейвлет-преобразования (ВП) в аналитической химии, химической физике и др. сходных областях. ВП применяется для сжатия данных, сглаживания и удаления шума (denoising), коррекции базовой линии и фона, разрешения наложив-шихся пиков (деконволюции) и др. Существенные положительные результаты применения ВП достигнуты в проточном инжекционном анализе (FIA), высокоэффективном капиллярном электрофорезе (CE), высокоэффективной жидкостной хроматографии (HPLC), инфракрасной спектроскопии (IR), масс-спектрометрии (MC), спектрометрах ядерного магнитного резонанса (NMR), электроаналитической химии и др. [3].

В борьбе с шумами применяются вейвлетные подходы подавления шума, основанные на работах Donoho D.L. и Johnstone I.M. [4]. Они предложили

метод подавления шума для случая аддитивной смеси сигнала с шумом (искажения, вызванные приборной функцией не учитываются), основанный на вычислении коэффициентов c■ (к) и dj (к)

вейвлет-разложения наблюдаемого сигнала у ()

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

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

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

у^) = |H(-т) х(т^т + u(). (1а)

(Здесь и в дальнейшем для простоты бесконечные пределы в интеграле опускаются).

Если прибор безынерционный, то H () = 5 (),

где 5 () — дельта-функция, и имеет место обычная смесь сигнала с шумом

у () = х(^ + u (1б)

В выражениях (1а, 1б) приборная функция H (), спектр частот и дисперсия шума с предполагаются известными.

Рис. 1. Информационная модель прибора: х (/)— полезный сигнал; Н (/) — импульсный отклик стационарной, инвариантной во времени системы (приборная функция); и () — стационарный

в широком смысле шум; у (/) — наблюдаемый сигнал

Развиваемый в настоящей работе вейвлетный подход направлен на решение задачи оценки сигнала, его производных или текущего значения интеграла, когда имеет место искажение полезного сигнала прибором вследствие его инерционности и шумом, т. е. для модели (1а).

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

х (') = (О * у)('),

или в образах Фурье

х (т) = О(т)у(т).

(2а)

(2б)

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

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

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

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

1. КРАТНОМАСШТАБНЫИ АНАЛИЗ И БЛОКИ ФИЛЬТРОВ

Так как изложение материала работы основано на теории кратномасштабного (мультиразрешаю-щего) анализа (КМА) сигналов, поясним основные ее положения, применив понятный инженерам спектральный подход к анализу сигналов. Пусть некоторое множество сигналов, имеющих частотный спектр в диапазоне от 0 до 21П, где П — некоторая круговая частота (рад/с), образуют подпространство пространства Ь2, которое мы обозначим как У1 (см. рис. 2). Пусть другое множество сигналов имеет частотный спектр (0,21 _ П), которое обозначим как У _1, т. е. при 1 = 0 У0 — (0,20П), У_ — (0,2-1 П), Р_2 — (0,2_2П) и т. д.

Очевидно, что все сигналы, которые принадлежат У1 _1, также принадлежат и У1, т. к. спектр сигналов из У _1 полностью входит в спектральный диапазон сигналов из У1. Условно это записывается как У1 _1 с У (говорят, У1 _1 вложено в У1). То же можно сказать и о подпространстве У0, т. е. У0 с у , и о других подпространствах У1. В результате подпространства, сконструированные приведенным выше способом, образуют цепочку вложенных друг в друга подпространств

...с У_2 сУ_1 сУ0 сУ сУ2 с....

(3)

Однако разбиение всей шкалы частот можно выполнить более рациональным способом, например нижние частоты сигнала сохранить в диапазоне частот (0,2_2П) и затем путем расширения спектрального диапазона добавлением частотных полос (2_2П,2_1 П), (2_1 П,20П), (20П,21 П) и т. д. добиться необходимой для данного сигнала общей ширины полосы 21 П (см. рис. 2). Множества сигналов, входящих в эти спектральные полосы, образуют подпространства вейвлетов и обозначаются как Ж_2,Ж_1,Ж0,Ж1 и т. д. Тогда подпространство,

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

например, У2 может быть представлено как сумма

У2 = У © Щ = У0 © щ0 © щ = = ... = у-2 ©Щ-2 ©Щ-1 ©щ0 ©щ

Или в общем виде

Уг = Уг , © щг , = У, © щ, © щ, , ©••• ©

J ° 1 J 1 у0 Jo Jo 1

©Wy © ... © WJ-2 © WJ

(4)

где J называют самым тонким, а у0 — самым грубым значениями масштаба. Подпространство Щу-1 называют ортогональным дополнением подпространства У._1 до У,. Это означает, что

У, ± Щу при всех у е 2.

Из выражения (4) следует, что любой сигнал из У}, например порожденный каким-либо изображением, можно рассматривать по частям: сначала его общий вид, содержащийся в низкочастотной части спектра У,0 , затем — добавив крупные детали, содержащиеся в Щ0, затем — все более мелкие детали по мере расширения спектра добавлением деталей, содержащихся в Щ,. Именно поэтому такой анализ сигнала получил название "математического микроскопа".

Заметим, что двукратное разбиение спектральных поддиапазонов — не единственный способ построения цепочки (3), хотя и весьма удобный для формирования теории и создания простых вычислительных алгоритмов [6]. В отличие от преобразования Фурье неравномерное разбиение частотной шкалы имеет определенный физический смысл: сигналы небольшой продолжительности

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

Определение: Последовательность вложенных друг в друга замкнутых подпространств (частотных поддиапазонов) (3)

...с V-2 cV-1 cV0 cV cV2 с....

называется КМА, если удовлетворяются следующие условия:

а) их объединение включает все сигналы из L2:

У Vj = L2, или F = L2;

yez

б) пересечение образует пустое множество

П Vj = {0} или V_={0};

jeZ

в) справедливо соотношение масштабируемо-

сти

f (t)e Vj « f (2t)e Vj+1;

г) найдется такая функция p e V0, которую называют масштабирующей функцией (scaling function), что множество ее сдвигов (p(t - к) образует базис Рисса или ортонормированный базис подпространства V0.

Поясним некоторые положения этого определения. Свойства а) и б) означают, что для любой функции xe L2 найдется такое подпространство Vj , что ее проекция на это подпространство образует сколь угодно малую ошибку приближения. Свойство в), по мнению И. Добеши [7], считается основной причиной, позволившей ввести термин кратномасштабность: оно означает (совместно со свойством (3)), что если x (t)e Vy, то одновременно x (t )e Vy+1, но в Vy+1 входит еще и масштабируемая версия x (t), а именно x (2t)e Vy+1, что

порождает масштабирующее уравнение, рассматриваемое ниже. Свойство г) совместно с в) означа-

ет, что для любого Vj существует ортонормиро-ванный базис

(, (t ) = 2j2 р

(

t - 2-]к

jk м

V (ю)=;/22

ю

где

m.

=moi — т т

= h (ю). V2 V ;

te—'2 т i || =

ю

штабирующую функцию р(г), т. к. ее Фурье-

образ можно получить из выражения (6а) последовательной подстановкой в эту формулу выражений для р\ т I, рI т I, рI т I и т. д. [7].

Причем для подпространства V0 этот базис будет иметь вид

т0,к (t )=((t - к) •

В дальнейшем будем пользоваться также обозначениями

(,к (t) = 2j2 т(2Ч - к) = ( (t - 2-Jk)•

Масштабирующее уравнение. Из определения КМА следует, что простым масштабированием и сдвигами одной функции (p(t) можно построить базисы для всех подпространств Vj . Однако КМА

накладывает на сами функции р(t) определенные и достаточно жесткие требования: поскольку (е V0 и V0 с V1, то функция (p(t) является линейной комбинацией функций р1п (t):

p(t) = Хh(n )( (t) = £h(n)^2p(2t-n), (5)

neZ neZ

которая называется масштабирующим уравнением (scaling equation, two-scale equation, refinement equation). Этому уравнению должны удовлетворять все масштабирующие функции КМА. Коэффициенты h (n) называются масштабирующими.

Они используются в качестве импульсных откликов в блоках (банках) фильтров, реализующих быстрые вычислительные алгоритмы дискретных вейвлетных преобразований (digital wavelet transform, DWT) — ДВП.

В образах Фурье уравнение (5) имеет вид:

2 У '{ 4 ) {

Учитывая, что р. к (г)еУ;. с У3+1, масштабирующее уравнение (5) для произвольного значения масштаба можно представить в виде

р(2Зг _ к) = ^к (пр(2(237 _ к)_п) =

п

= £к(п)^2 р(+1г _(2к + п)),

п

которое после замены переменных 2к + п = т становится

р С г _ к) = ^к (т _ 2к)Лр (23+1 г _ т),

или

(к (t) = Еh (m - 2k)(+1,m (t) .

(7а)

Выполнив преобразование Фурье, получим

р. (т) = к (2_3_1 т +1 (т), (7б)

где

р. (т) = 2_32р(2_3т).

Отметим, что

р3кк (т)=р3 (тУ

-1'ю2 J к

(7в)

(6а)

Из (7а), умножив слева и справа скалярно на функцию р. п(г), получим

к(п _ 2к) = (23Ч/2р(( г_к),р(23+1 г _п)) = = ( л/2 р(г), р( _(п _ 2к)). (8)

Вейвлетные функции. Базисом для подпространства вейвлетов Ж. являются сдвиги и масштабирование одной вейвлетной функции у (г)

(

t - 2-jk

w

(6б)

Функция m0 (ю) полностью определяет мас-

ji¥], (t) = 2з/2¥ j е Z, k е N.

В дальнейшем будем пользоваться также обозначениями

¥jk (t) = 2j2i( t-k) = ¥j (t - 2-Jk).

Так как Щ0 е у , то функция ) может быть представлена в виде взвешенной суммы масштабирующих функций р^)

V(t) = ^g (п)р(2 _ п), (9)

п

где

8п =(_1)1-пй (1 _ п), или в образах Фурье

у ^ (п ^ 2 р (ю )= )( (ю).

где

т (ю) = ^2,8(пУ'пю=^8(ю) =

= в,ш 42т0 (ю + п).

Учитывая, что у,к (t)еЩу с VУ+1, выражение

(9) для произвольного значения масштаба можно представить в виде

у (2уt _ к) = 28 (т _ 2кр(2у+11 _ т),

п

или

У ^) = 2 8 (т _ 2к(0 . (10а)

т

Выполнив преобразование Фурье, получим

у/. (ю) = 8(2_у-1 ю )(,+1 (ю), (10б)

где

уу (ю) = 2_у2у (2_ую).

Отметим, что

уу (ю) = у (ю)е-,ю2-Ук. (10в)

Из (10а), умножив слева и справа скалярно на функцию 2у2 (р (2уЧ11 _ п), получим

8 (п _ 2к) = ^2]у[2у(2Ч _ к),р(2у +11 _ п) =

^ ^/2у(t), р( _(п _ 2к))). (11)

Анализ сигнала. Из (4) следует, что для любой функции х ()е УJ

х ()=2 сУ0(к ( () + 2 2 (к К* ^),

к У =У0 к

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

к (_п)

у(п) -^ CJ_1

(_п)

к(_п) к(_п)

(_п)

(_п)

а

к (п) к (п) к (п)

CJ _1 <-CJ _2 <-

'(п)

:(п) ^ч 8(п)

б

Ч +1

к (_п)

8 (_п)

к (п)

с,0 +1

Рис. 3. Схемы 4-ступенчатых анализирующего (а) и синтезирующего (б) блоков фильтров. к и 8 — коэффициенты низкочастотного и высокочастотного фильтров вейвлетов

J _2

J _3

J _1

J _2

J _3

J _3

J _1

J _3

J _2

J _1

сигнала по базисным функциям р.к (г) и у.к (г) соответственно

Сз (кИх,р3к) , ёз (кИХ,Узк) .

Умножив выражение (7а) скалярно слева и справа на х (г), получим:

(k ) = Х h (m - 2k )ji (k).

(12а)

2п

-Jx (ю)е~'т2 Jk dю = const• x(k2-J),

что при соответствующей нормировке обеспечивает выполнение равенства (13).

На практике вычисление с. (к) и ё. (к) по формулам (12а, 12б) осуществляется с помощью

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

(обозначим его как h (n))

y (k ) = S h (k - n )x (n).

(14)

Аналогично, умножив выражение (10а) скаляр-но слева и справа на х (г), получим:

ё3 (к) = Е 8 (т _ 2к )с.+1 (к). (12б)

т

Формулы (12) являются рекурсивными алгоритмами анализа сигнала "от тонкого к грубому разрешению", т. к. каждые последующие коэффициенты вычисляются из предыдущих с большим масштабным индексом, характеризующим более широкополосный процесс [7]. Вычислительная схема алгоритма приведена на рис. 3, а. Она реализуется с помощью многокаскадного блока последовательно соединенных фильтров, обеспечивающих быстрые вычисления ДВП. Начальные, значения коэффициентов с. (к) вычисляются при таком значении масштаба 1 , когда в пределах спектрального диапазона сигнала х (т) функция р1 (т) остается практически постоянной (базисная функция по отношению к х (г) имеет вид 5-функции). При этом значении масштаба коэффициенты с1 (к) можно положить равными отсчетам функции

х(г)| ;=к2- = х(к2_1 ) = С, (к). (13)

Действительно,

С, (к) = /х(гр к (г= = 2-/ х (т)р1 (т)е~'тТ 1 к ат = (0)

Для того чтобы этот фильтр осуществлял требуемую в соответствии с (12) свертку, необходимо инвертировать порядок коэффициентов h (n), т. е. положить h (k) = h (-k). Тогда будем иметь под знаком суммы в (14) h (k - n ) = h (n - k), что требуется на первом этапе вычислений. Далее в соответствии с формулами (12) необходимо убрать из выходной последовательности фильтра те отсчеты, которые имеют нечетные индексы, сохранив только четные, осуществив, таким образом, с помощью устройства, называемого дециматором, операцию понижения числа отсчетов вдвое (down-sampling). В результате получим, что для реализации алгоритмов (12) необходим блок фильтров, содержащий высокочастотный и низкочастотный каналы, как показано на рис. 4. Чтобы осуществить многомасштабный (мультиразрешающий) анализ необходимо применить многокаскадный блок фильтров, построенный по схеме рис. 3, а.

Синтез (восстановление сигнала по коэффициентам вейвлет-разложения). Получим далее рекурсивные формулы для восстановления функции

1

Рис. 4. Однокаскадный блок фильтров анализа, реализующий рекурсивную процедуру вычисления вейвлет-коэффициентов сигнала. к (_п) и 8 (_п) — низкочастотный и высокочастотный фильтры с обратной последовательностью коэффициентов масштабирующих уравнений для р и у ; 2 — дециматор, удаляющий вейвлет-коэффициенты с нечетными индексами на выходах фильтров

х () (точнее, ее отсчетов х (п)) по коэффициентам вейвлет-разложения су (к) и ёу (к). Учитывая, что Уу+1 = У у © Жу базисную функцию ру+1 к ()е У у+1 можно представить как

cj(m )=2h (m - 2k )cj (k)+

k

+ 2g(m - 2k)c: (k).

(16)

dJ

t:

s(n)

( +1, т ( ) = ^{Ру +1, т (у, к) Ру, к ( ) +

к

+ ^(Ру +1, т ,¥у, к )¥у, к ( ). к

С учетом формул для скалярных произведений под знаком суммы (8) и (11) получим

Ру+1,т () = Ек (т - 2к)Рук () +

к

+ Е8 (т - 2к)¥у,к (). (15)

к

Умножив далее выражение (15) скалярно слева и справа на х (), получим рекурсивное соотношение для вычисления коэффициентов вейвлет-раз-ложения "тонкого" масштаба у +1 как взвешенной

суммы с у (к) и ёу (к) более "грубого" масштаба у :

-» U -► h(n)

1

Вычислительная схема восстановления первоначальных отсчетов сигнала по коэффициентам его вейвлет-разложения показана на рис. 3, б. Чтобы выполнить вычисления по формуле (16) с помощью цифровых фильтров, осуществляющих свертку коэффициентов h и g с отсчетами Су (k),

необходимо эти фильтры дополнить еще одним устройством, называемым интерполятором. Он осуществляет повышение числа отсчетов вдвое (up-sampling) путем прореживания их нулями: все отсчеты с нечетными индексами на выходе интерполятора оказываются равными нулю. Тогда свертка коэффициентов фильтров h и g благодаря двойному сдвигу 2k в формуле (16) будет осуществляться со всеми значащими отсчетами Су (k)

и dy (k). Схема одного каскада блока фильтров

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

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

Рис. 5. Однокаскадный блок фильтров синтеза, реализующий рекурсивную процедуру вычисления вейвлет-коэффициентов сигнала. 2 — интерполятор, прореживающий вейвлет-ко-эффициенты нулями на входах фильтров

3. СИНТЕЗ МОДИФИЦИРОВАННЫХ ВЕЙВЛЕТОВ

Чтобы построить КМА для решения задач оценки сигнала, искаженного аппаратной функцией и шумом, и/или его линейных преобразований, необходимо синтезировать базисные функции соответствующих подпространств, ориентировав их на решение требуемой задачи. Теоретически этот синтез выполняется путем модификации известных вейвлетов, имеющих компактный носитель в частотной области, в частности вейвлетов Мейера, с помощью вещественной функции О (^ )е Ь2, достаточно быстро затухающей при ^^ . Практически в качестве вейвлетов-"прародителей" могут быть выбраны вейвлеты Добеши высоких порядков, характеризующиеся плоской частотной характеристикой в области пропускания и хорошей крутизной фронтов в области затухания.

Свертка функции О () с наблюдаемым сигналом у () дает требуемую оценку

х(1) = \о(1 -т) у (т) ёт, (17а)

или в образах Фурье

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

Представим выражение (1а) в образах Фурье:

у (ю) = Й (ю)х(ю) + и (ю). Подставляя это равенство в (17б), получим

х(ю) = О(о))Й (ю)х(ю) + О(ю)и (ю).

В частности, если выбрано О (ю) = Й-1 (ю), имеем

оценку полезного сигнала в виде

х(ю) = х(ю) + Н-1 (ю)и(ю).

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

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

^ )=Е (к К,к()+

к

+ ЕЕё!(к >()

] = ,0 к

(18)

вейвле-

где ф,о * (), () — модифицируемые ты; с,. (к) и ё, (к) — коэффициенты вейвлет-раз-ложения сигнала у ().

Для того чтобы найти оценку х (^), необходимо определить коэффициенты с^ (к) и ё, (к) по

наблюдаемому сигналу у ().

Выполним преобразование Фурье (18):

х (ю)=(ю) Е С0 (к) ехр (-ю2-л к)+

к=-"

" ^

+ Е (ю) Е ё, (к) ехр(-Ш2-1 к),

1 = ,0

где

у (ш) = О-1 (ф]о (ш) Е с,0 (к) ехР(-/°к)+

к=-"

„ " ^

+ Е О-1 (ю) (ю) Е ё, (к) ехр(-¡юк).

]=, к=-"

Введем функции в,к ()=в,( - к) и Пк (^)=п( - 2-'' к) (,,к ё2) такие, что

к (ю) = 0-1 (ю)фк (ю) =

= О-1 (ю)ф, (ю)ею2-к, (19а)

Пк (ю) = 0-1 (ю)^,,к (ю) =

= О4 (ю)> (ю)2-, (19б)

и выполним обратное преобразование Фурье у (ю), получим

у (' )= Е С0 (к )в0,к ()+

к=-"

/ "

+ ЕЕ ё, (к )п, к ().

,=,0 к=-"

(20)

Предположим, что существуют двойственные функции в, к () и г), к () такие, что

в:,к (ю)= 0(ю)Ф},к (ю) =

= О(ю)ф , (ю)еЧсо21 '

(21а)

П,к (ю)= С(ю)>. к (ю) = = О (ю)>; (ю^-'"2"

(21 б)

ф}0 (ю) = 2 2 ф (2-л ю) и (ю) = 2 2> (2-,ю).

Умножим последнее выражение слева и справа на функцию О- (ю), получим

Легко проверить, что функции в, к () и в, к () биортогональны по отношению к собственным сдвигам, а функции г],к () и г),* () биортого-

нальны по отношению к сдвигам и масштабу. Кроме того, биортогональными являются функции в!,к () и , (), а также в,,к () и Пк ().

Для того чтобы равенство (20) было справедливым, необходимо, чтобы функции в и ц , а также

в иг) являлись базисом соответствующих подпространств КМА. Заметим, что преобразования (19а, 19б) и (21а, 21 б) осуществляют линейное отображение подпространств КМА V , и Ж, в новые подпространства, как показано ниже:

V ——^ и

V ——>

Ж ——-и М

Ж ——> М

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

и, и и,, а функции ц и т) — образами вейвлет-

ной функции у в подпространствах М, и М..

и

°1

А-——2

1 1 \ з / \ 1 4 1

1 1 \ \ / \

Уч /Ч СО ! 71

-0.5

0.5

¿о. с«

О 2

б

7-4,0 (*)

0.05

"044

-0.05

-0.1

Рис. 6. Синтез модифицированных вейвлетов при двух значениях масштаба у = 0 (а, б, в) и у =-4 (г, д, е): 1 и 4 — 8(ю), 2 — О(ю), 3 — к (ю)

а

г

е

в

Согласно требованию г) определения КМА, сдвиги этих функций должны образовывать базис Рисса в соответствующих подпространствах. Следуя работе [8], можно показать, что функции в формируют базис Рисса подпространства и у, функции в — базис Рисса подпространства О у, функции п — базис Рисса подпространства Му и функции п — базис Рисса подпространства Му. Более того, функции п и п порождают две цепочки КМА:

-•си-2 си-1 си0 си1 си2 с•••

и

-с и-2 с и-1 с ио с й с и2 с -,

которые в свою очередь порождают вейвлетные подпространства {Му} и {Му } . При этом справедливы следующие соотношения биортогональности

иу 1 иу; иу 1 Му; и у 1 Му; Му 1М (22а)

и вложенности

V -1 с и с у+1; у -1 с О у с у,. (22б

В результате имеет место соотношение Ь2 (Я ) = и. © М © М +. ©... © М ©... © Мг, и

V / у0 у0 у0 +1 у

равенство (20) оказывается справедливым.

На рис. 6 показан пример построения модифицированных вейвлетов для некоторой типовой

функции О (ю) при двух значениях масштабов у = 0 и у = -4 . Как следует из рис. 6, д, е, уже при значении масштаба у = -4 форма масштабирующей и вейвлетной функций близки к соответствующим функциям-"прародителям".

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

- образуют биортогональную систему базисных функций, формирующих кратномасштабный анализ сигналов в Ь (Я);

- позволяют получить линейную оценку полезного сигнала;

- нестационарны по отношению к масштабу. Последнее вытекает из (19а, 19б) и (21а, 21 б), т. к. в зависимости от у изменяется носитель функций-

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

уже не является результатом растяжения/сжатия одной функции, как в случае традиционного КМА —

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

Для "всепропускающих" фильтров G (ra) = const

функции в,в,ц и f вырождаются в соответствующие функции-"прародители" р и у .

4. МАСШТАБИРУЮЩИЕ УРАВНЕНИЯ

МОДИФИЦИРОВАННЫХ ВЕЙВЛЕТОВ

Чтобы определить коэффициенты c у (k) и

dy (k), умножим скалярно слева и справа выражение (20) на двойственные базисные функции вук (t), fjyk (t) и с учетом биортогональности (22а) получим:

Со (kИy (t )! (tК (23а)

dy (k) = Jy(t)fyk (t)dt. (23б)

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

Так как подпространства {йу-} формируют КМА и ву+1 е иу+1 = иу © Му , то базисные функции ву. и Пу аналогично масштабирующим уравнениям (7а) и (10а) могут быть выражены в виде линейной комбинации ву+1. Для функции ву имеем

у (t) = Шk,ву 4®у +1,- (t) .

n

Определим скалярное произведение под знаком суммы:

{Оу, ,ву+и) = Jву ( - 2-k )+1 ( - 2-у-1 n )dt =

= 2~JС- (т)е~,ш21 k ! (m)eiaT1-n dra =

= J Р(2— Ю)(у+1 (m)e'm2'1 (n-2k) dra =

1 f 1 „, ч fra^ , | (n-2k)

= — —¡=р(га)р| — Ie 2 dra =

2П 42 ;2 J

= Jp(t)<Jlp(2t - (n - 2k))dt = h (n - 2k).

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

к (<) = Ек (п - 2кв+1, (х), (24а)

п

или

в, (ю) = к (2-]-1ю)§+1 (ю), (24б)

где в, (ю) определяется из (21а).

Аналогично для функции в, к (х) получим:

вк () = Ек (п - 2к)в+1,п (Г), (25а)

п

или

в , (ю) = к (2-1-1 ю))+1 (ю), (25б)

где в, (ю) определяется из (19а).

Из (25а), умножив скалярно слева и справа на функцию в ,+1 т (х) с учетом биортогональности, получим

к(т - 2к) = (в^тв*) . (26)

Далее построим масштабирующее уравнение, которому должна удовлетворять вейвлетная функция Л,,к (х)

п, к(х)=Е(п к ,п) +1,п(х).

п

Определим скалярное произведение под знаком суммы:

Пв = (ю)ею к ^(Щ)е'"2^-1 п (ю =

= 2п > (ю)Ф'+1 (ю)е'ю2'1(п-2к)аю =

= 2_\ (2-у ю)ф (2"l'"1ю)eiю2-1-1(n"2 к) аю 1 Г 1 V С 'юю(п-2к)

= — \—=у(ю)ф| —\е 2 аю = 2П у/2 У } \ 2 J

= \У (х) 72 < (2х - (п - 2к))х = ^ (п - 2к). Следовательно,

Пк () = Е£ (п - 2кв+1,п (х), (27а)

п

или

П (ю) = £ (^ю).+1 (ю). (27б)

Аналогично для функции п,к (х):

п,к (х) = Е£ (п - 2к)в+1,п (х), (28а)

п

или

П (ю) = £ (-ю)+1 (ю). (28б)

Из (28а), умножив скалярно слева и справа на функцию (х), с учетом биортогональности

получим

£(т - 2к) = (0в,1,т,п-к) . (29)

Справедливость уравнений (24а, 24б), (25а, 25б), (27а, 27б) и (28а, 28б) легко проверить, подставив

в них вместо функций в,6,ц и г) их выражения из (19а, 19б) и (21а, 21 б) и сравнив результаты соответственно с (7а) и (10а).

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

и, с У,.+1, в, е и, и ф,+1 е У,.+1, то существует такая 21+1 2п -периодическая функция а.+1 (ю), что

в, (ю^О,+1 (ю)<р,+1 (ю). (30)

Выполним обратное преобразование Фурье выражения (30). Так как

!+1

2 2 ф (2~-''~ 1ю), имеем

в ; (х) = Е, (п)2-¡ф^^'-п)((ю.

п

После замены переменной ю1 = 2"1"1ю

в, (х ) = Еаа,+1 (п )\ф (ю )е'ю (+1х-п ю1 =

п

= Еа+1 (п)2^+1 X - п)

п

и окончательно:

в, (х) = Еа+1 (п)ф,+1 (х - 2-'''-1 п) =

п

= Еа!+1 (п)Ф+1,п (х). (31)

п

Для произвольного значения сдвига к формула (31) имеет вид

в (ю) =

Е^^,+1(п )е

у+1

! (X) = Ъ&у +1 (п)2 2 р(2у +1 X - 2к -п),

п

а после замены переменных т = 2к + п становится в у к (X ) = Ха+1 (т - 2к )Ру+1,т (X). (32)

т

Получим формулу для вычисления йу+1 (п). Имеем в соответствии с (21а, 21б) с учетом (7б):

ву (ю) = О(ю)к (2-у-1ю)(ру+1 (ю),

где произведение первых двух членов образует 2у+1 2п -периодическую функцию. Следовательно,

и

а

'у +1

ау+1 (ra)= G(ra)h (2 у 1ra)

(n ) =

(33а)

, 2у+1п

-- J G(ra)h (2-у-1ra)ei2-1-1шndra. (33б)

1

2у+1 2n

Аналогично (30) для функции пу (t) буде иметь

м

пу (ю) = ву+1 (ю)ру+1 (ю) (34)

и, следовательно,

пу (X) = у +1 (п)Ру +1 (-2-у-1 п) =

п

= 2! (п )Ру +1,п (X). (35)

п

Аналогично (32) имеем

пу,к (X) = Х^ +1 (т - 2кР +1,т (X). (36)

т

Получим формулу для вычисления ¡Зу+1 (п). Имеем в соответствии с (21 б) с учетом (10б)

п\3 (ю) = О(ю)уг}. (ю) = О(ю)8(2-у-1 ю) РРу+1 (ю),

где произведение первых двух членов образует 2у+1 2п -периодическую функцию. Следовательно,

+1 (ra) = G (ra)g (2-у -1ra)

(37а)

и

Ру+1 (n ) =

1 у1 ж

J G(ra)g(2 --1ra)e"2 J 'ffl"dffl. (37б)

ЗАКЛЮЧЕНИЕ

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

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

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

1. Barclay V.J., Bonner R.F. Application of wavelet transforms to experimental spectra: smoothing, denoising, and data set compression //Anal. Chem. 1994. V. 69. P. 78-90.

2. Меркушева А.В. Классы преобразований нестационарного сигнала в информационно-измерительных системах. III. Время-масштабные (вейвлет-) преобразования для спектрально-временного анализа // Научное приборостроение. 2002. Т. 12, № 3. С. 68-81.

3. ShaoX.-G, Leung A.K.-M., Chau F.-T. Wavelet: a new trend in chemistry // Accounts of Chemical Research. 2003. V. 36, N 4. P. 276-283.

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

5. Mallat S. A wavelet tour of signal processing. Academic Press, 2000. 637 p.

6. Burrus C.S., Gopinath R.A., Haitao Guo. Introduction to wavelets and wavelet transform. Prentice Hall, 1997. 268 p.

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

8. Zhang J., Walter G. A wavelet-based KL-like expansion for wide-sense stationary random processes // IEEE Transaction on Signal Processing. 1994. V. 42, N 7. P. 1737-1745.

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

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

MODIFIED WAVELETS IN DATA PROCESSING FOR ANALYTICAL INSTRUMENTS.

I. BASIC THEORY

L. V. Novikov

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The theory and practice of wavelet transform based experimental data processing are considered. The mul-tiscale analysis is treated here as broadening band signal spectral analysis. Filter banks implementing fast computational algorithms are constructed, and a class of non-stationary wavelets obtained by modification of known orthonormal wavelets by the instrument apparatus function is studied.