Яковлев А.Н. Введение в вейвлет-преобразования (2003), страница 5

Формально обобщенный ряд Фурье (2.3) отличается от традиционноготем, что суммирование проводится не по одному, а по двуминдексам. Однако это несущественно, так как обе системы индексации принадлежат одному классу бесконечных счетныхмножеств.Диадное ВП часто называют дискретным. Однако, по мнению ряда авторов, например В.П. Дьяконова [8], такая подменаформулировки не совсем корректна: правильнее называть егодиадным, представляющим особую разновидность непрерывного ВП и позволяющим устранить избыточность последнего.Примечаниe. В некоторых публикациях параметрыфункции задаются в видеa , b и базисныеa = 1/ 2 j , b = k / 2 j , ψ jk (t ) = 2 j ψ (2 j t − k ) ,(2.1’)т.е.

с ростом j параметр a уменьшается, что соответствует сжатию функцииψ jk (t ) . Согласно формулам (2.1) с ростом m увеличивается и коэффициентa , т.е. функция ψ mk (t ) растягивается.Фреймы. Это особый вид вейвлетов, занимающих промежуточное положение между непрерывным и диадным ВП. Вейвлет-фреймы используюткратное двум масштабирование ( a = 2m ), но непрерывный сдвиг. Следовательно, они сохраняют избыточность, которая присуща непрерывному ВП, нов гораздо меньшей мере по сравнению с ним. Они не входят в пакеты расширения систем компьютерной математики (СКМ). Но если необходимо, то соответствующие им инструментальные средства легко получить незначительноймодификацией средств непрерывного ВП.342.2.

Дискретное преобразованиеСтатьи, касающиеся практического использования ВП,содержат в основной своей массе результаты компьютерныхрасчетов, в которых использовано дискретное вейвлет-преобразование (ДВП или DWT). При этом не только параметры a иb , но и сигналы также дискретизируются во времени.На основании теоремы Котельникова (теоремы отсчетов)непрерывный сигнал S (t ) , спектр которого не содержит частотвыше f m , полностью определяется дискретной последовательностью своих мгновенных значений {Si } , i = 0,1, ..., N − 1 , отсчитываемых через интервалы времени Δt :Δt = 1/ 2 f m , f д = 1/ Δt = 2 f m ,(2.5)где Δt и f д – интервал (шаг) и частота дискретизации.Таким образом, дискретизированный с шагом Δt сигналможно определить выражением:Sд (t ) = {Si } =N −1∑ S (iΔt )δ(t − iΔt ),(2.6)i =1где δ(t ) – дельта-функция.Если число отсчетов составляет N = 2n0 , то максимальноезначение m в формулах (2.1) будет равно n0 − 1 .

Наибольшеезначение k для текущего m определяется: k = 2n0 − m − 1 . В частности, для m = 0 (т.е. a = 1) число сдвигов k базисного вейв-лета составит 2n0 − 1 = N − 1 ; с каждым последующим значением m (1, 2, …) вейвлет ψ mk (t ) расширяется в два раза, а числосдвигов k уменьшается в два раза. Для максимального значения m = mmax , равного n0 − 1 , k = 0 , т.е. один вейвлет ψ mmax 0 (t )«накрывает» весь интервал сигнала (рис. 2.1; N = 8 ).Вейвлет-коэффициенты cmk (или c jk ) можно вычислить спомощью итерационной процедуры, известной под названиембыстрого вейвлет-преобразования БВП [1–3, 19, 25, 29]. Алгоритм БВП приведен в п. 2.4.

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

В ядро систем(начиная с версии Mathcad 8) встроен единственный вейвлет –Добеши db4 (или DB4). При этом реализация ВП происходит сбольшой скоростью (т.е. эффективностью) и можно осуществлять практическое исследование различных сигналов и временных рядов на выявление как их свойств, так и свойств ВП.Ядро систем Mathcad содержит две следующие функции ВП:wave(x) – вектор прямого ВП;iwave(w) – вектор обратного ВП.Вектор данных x и вектор вейвлет-спектра w должны иметьровно N = 2n0 элементов ( n0 – целое число).

Результатом функции wave(x) является вектор, скомпонованный из коэффициентов двухпараметрического вейвлет-спектра cmk .Пример 2.1. Прямоугольный импульс с шумомИсследуемый сигнал x(t ) представляет собой аддитивную смесьx(t ) = S (t ) + n(t )прямоугольного видеоимпульса S (t ) и белого нормального шума(рис.

2.2):U = 5 (В), t0 = 40 (мкс), τ = 60 (мкс)s(t):= U ift 0 ≤ t ≤ t 0 +τ0 otherwisen(t ).Такая модель может характеризовать в первом приближении сигнал в видеотракте приемника радара (радиолокатора), сонара (гидролокатора) и оптического локатора.367xi−276543210120020406080100 120 140 160 180 200 220 240 260i260Рис. 2.2Представление сигнала и шума в дискретном виде:n0 = 8 , N = 2n0 , N = 256 , i := 0..N − 1 , si := s (i )σ := 0,3ni := σ −2ln(rnd (1)) sin(2π rnd (1)) , xi := si + niВейвлет-анализ, т.е. прямое ДВП:i := 0..N − 1 y := x w := wave( y )levelcoeffs (level ) := submatrix( w, 2level +1,2z := n0 − 1 z := 7 m := 1, 2..z− 1,0,0) ci , z − m := coeffs (m)⎡⎤⎢⎥i ⎥⎢flor⎢⎛ N ⎞ ⎥⎢⎜ m ⎟ ⎥⎣⎢ ⎝ 2 ⎠ ⎦⎥Семейства коэффициентов вычисленного вейвлет-спектра показаны нарис. 2.3, а весь спектр – на рис.

2.4.mПримечание. У коэффициентов (c )i нижний индекс i означает номертекущего отсчета времени и принимает N значений от 0 до N–1, а верхний mимеет тот же смысл, что и у вейвлет-коэффициентов cmk , определяемых поформуле (2.2). Напомним, что параметры m и k (которым соответствуютиндексы вейвлет-коэффициентов) характеризуют дискретные изменения временного масштаба ( a = 2m ) вейвлета и его сдвига (b = k ⋅ 2m ) во времени. Длятекущего масштаба m параметр k имеет 2n0 − m значений от 0 до 2n0 − m − 1 .

Вчастности, для m = 0 ( a = 1 ) вейвлет ψ 0 k ( x) смещается N раз (включая0нуль), т.е. индекс k в cmk и индекс i в (c )i совпадают. При m = 1 вейвлетψ1k ( x) расширяется по сравнению с вейвлетом ψ 0 k ( x) в два раза и общеечисло сдвигов будет в два раза меньше; при этом значение k будет изменяться через два отсчета i .

Для наибольшего временного масштаба, когдаm = n0 − 1 (в данном случае 7), k = 0 и один вейвлет ψ 7,0 ( x) «накроет» весьвременной интервал; при этом значение (cc7,0 при всех значениях i от 0 до N – 1.377)i будет постоянным и равным2.5( c 〈0 〉 ) i( c 〈1 〉 ) i−01.5(c2〈2 〉005000500050100150200150200i2502565)i( c 〈3 〉 ) i0( c 〈4 〉 ) i−6.20(c )i〈5 〉( c 〈6 〉 ) i( c 〈7 〉 ) i5100i2502562010010−15100150iРис. 2.3cРис.

2.438200250256m := 07xix1i−2765432101276543210127x1i−27x1i−27x1i−200007654321012406080 100 120 140 160 180 200 220 240 260i260m := 120406080100 120 140 160 180 200 220 240 260i260m := 3007654321012200020406080100 120 140 160 180 200 220 240 260i260m := 420406080100 120 140 160 180 200 220 240 260i260Рис. 2.5Вейвлет-синтез, т.е. обратное ДВП. Синтезируемый сигнал:x1i := iwave( w) .39Осуществим синтезирование сигнала с подавлением коэффициентов cmkпри быстрых (высокочастотных) слагаемых обобщенного ряда (2.3):j := 2 z − m... N − 1w j := 0 .Результаты представлены на рис. 2.5. Очевидно, что при m = 0 синтез происходит без подавления составляющих и исследуемый xi и синтезируемый x1iсигналы полностью совпадают.С увеличением параметра m расширяется полоса подавления составляющих в вейвлет-спектре, что эквивалентно пропусканию сигнала через фильтрнизких частот с уменьшающейся полосой пропускания фильтра и, следовательно, росту подавления шума и относительно высокочастотных компонентовсигнала; последнее приводит к искажению (затягиванию) фронтов импульса.Пример 2.2.

Вейвлет-фильтрация бигармонического сигнала с шумомНа вход низкочастотного приемного тракта фазового параметрическогогидролокатора [38] поступает сигналx(t ) = n(t ) + sэ (t ) ,где n(t ) – помеха в виде белого гауссова шума с математическим ожиданиемmn = 0 и среднеквадратическим отклонением σ ; sэ (t ) – бигармонический(двухкомпонентный) сигналsэ (t ) = A1 sin[2πF1t − ϕ1(t )] + A2 sin[2πF2t − ϕ2 (t )] , t з < t < t з + τи ,где A1, A2 – амплитуды компонентов эхосигнала на частотах F1 = F иF2 = 2 F , ϕ1 (t ), ϕ2 (t ) – фазовые сдвиги, зависящие от акустической жесткости и структуры подводных объектов, tз и τи – время задержки и длительность эхоимпульса. Измеритель фазового сдвига приемника дает на своемвыходе: во время действия сигнала sэ (t ) напряжение, пропорциональное фазовому сдвигу Ψ (t ) = 2ϕ1 (t ) − ϕ2 (t ) , и равномерный шум вне интервалаtз < t < tз + τи .

Присутствие шума n(t ) привносит случайную ошибку в измерения Ψ (t ) .Задача ВП – это осуществление фильтрации sэ (t ) из шума n(t ) . Отфильтрованный сигнал подается на измеритель фазового сдвига (ИФС на рис. 2.6).y (t )x(t )ИФСВПΨ (t )Рис. 2.6Моделирование сигнала x(t ) , его ВП, алгоритма измерения фазовогосдвига ψ (t ) и оценки параметров входного и выходного сигналов осуществлено в пакете Mathcad (2001).Дискретное ВП осуществлено на основе встроенной в пакет Mathcad базисной функции Добеши db 4 : wave(x) – вектор прямого ВП, iwave(w) – вектор40обратного ВП.

Сигнал x(t ) был подвергнут (как и в примере 2.1) прямому ВП,и по найденному вектору коэффициентов cmk осуществлено обратное ВП сподавлением коэффициентов при быстрых (высокочастотных) слагаемых ряда(2.3): yi := iwave( w) , j := 2 z − m..N − 1 , w j := 0 .На рис. 2.7 представлены входной сигнал x(t ) и результат его вейвлетфильтрации y (t ) . Здесь A1 = A2 = 1 В, F = 10 кГц, ϕ1 (t ) = ϕ2 (t ) = ϕ = 90° ,σ = 0,5 В, t з = 0, 2 мс, τ = 0,6 мс, n0 = 10 (z = 9), N = 1024 .В ходе исследований изменялся параметр m и измерялись: среднеквадра% , отклонениетическое значение шума σ y после ВП, фазовый сдвиг Ψ%%Δ = Ψ − ϕ измеренного фазового сдвига Ψ от истинного Ψ = 2ϕ − ϕ = ϕ иψ12среднеквадратическое отклонение фазового сдвига σψ . Результаты сведены втабл.