Яковлев А.Н. Введение в вейвлет-преобразования (2003) (1245706), страница 6
Текст из файла (страница 6)
2.1.При m = 0 восстановление (синтез) сигнала происходит без подавленияшума и синтезируемый y1 и исследуемый xi полностью совпадают. С увеличением m (на рис. 2.7 m = 3 ) расширяется полоса подавления спектра, сужается полоса пропускания низкочастотной части спектра, что уменьшает уровень шума ( σ y ).Рис. 2.7mТ а б л и ц а 2.12301456σy , В0.4930.3270.1980.1250.1050.0780.05Δ ψ , град–5.73–5.35–5.01–4.24–4 .84–6.91120.σψ , град14.8714.8118.4714.0814.5418.8116Установлено, что существует оптимальное mopt , при котором отклонениеΔ ψ , или среднеквадратическая ошибка σψ , в измерении фазового сдвига эхо-41сигнала минимально. При m > mopt наряду с дальнейшим подавлением шумапроисходит искажение формы бигармонического импульсного сигнала, чтоприводит к росту Δ ψ и σψ . Отметим, что при отсутствии шума ( σ = 0 ) измерение фазового сдвига происходит без ошибок (т.е.
Δ ψ = 0 и σψ = 0 ).2.3.2. ДВП в MatlabИзучение ДВП, как и непрерывного ВП, лучше осуществлять с помощью графического интерфейса GUI (прил.1). Длявызова меню следует исполнить команду wavemenu. В появившемся окне со списком разделов ВП активизировать позициюWavelet 1-D, а далее установить File → Demo Analys и выбратьодин из 32 примеров применения вейвлет-технологии.Пакет расширения систем MATLAB 6.0/6.1 Wavelet Toolbox2/2.1 содержит несколько функций нахождения вейвлеткоэффициентов (прил.
2, П.2.2), например,coef = detcoef(C,L,M).Эта функция возвращает коэффициенты на уровне М из структуры wavelet разложения [C, L]; при этом уровень М долженбыть целым числом, таким, что 1 ≤ M ≤ MMAX , гдеMMAX = length( L) − 2 .Функция [С, L] = wavedec(S, M, 'wname')возвращает векторы wavelet разложения сигнала Х на уровне М,используя выбранный вейвлет (с именем ‘wname’).Пример 2.3. Бигармонический сигнал с шумомМодель такого сигнала приведена в примере 2.2. Найдем его непрерывныйи дискретный вейвлет-спектры:function binar_rauch_wav_1t = 0:0.000001:0.001; A1 = 1; A2 = 1; F1 = 10000; F2 = 2*F1; a1 = 0; a2 = 0;s1(1:200) = 0; t2 = 0.0002:0.000001:0.0008;s2 = A1*sin(2*pi*F1*t2 + a1) + A2*sin(2*pi*F2*t2+a2);s3(1:200) = 0; s = [s1 s2 s3]; randn('state',0); g = 0.5; n = g*randn(size(t));x = s + n; figure (1); subplot(311), plot(t,x,'k'); title('Сигнал x(t)'); grid on;gtext('F=10кГц, А1=А2=1В, g=0.5 B');subplot(312), c = cwt(x,1:64,'mexh','absglb',[0 400]);title('Вейвлет-спектр'); xlabel('Временной сдвиг, b');ylabel('Временной масштаб,a'); set(gca,'Xlim',[0 1000]);[c,l] = wavedec(s,6,'db4');for m = 1:6d = detcoef(c,l,m); d = d(ones(1,2^m),:);cfd(m,:) = wkeep(d(:)',1000);end42cfd = cfd(:); I = find(abs(cfd)<sqrt(eps));cfd(I) = zeros(size(I)); cfd = reshape(cfd,6,1000);subplot(313), colormap(pink(16));img = image(flipud(wcodemat(cfd,64,'row')));set(get(img,'parent'), 'YtickLabel',[]);title('Дискретное преобразование');ylabel('Уровень, m'); xlabel('Временной сдвиг, b'); endРис.
2.8На рис. 2.8 приведены диаграмма сигнала и его спектрограммы. Особенности спектрограммы непрерывного ВП этого сигнала обсуждены в примере1.6. Очевидно, что детали сигнала просматриваются и на спектрограмме дискретного ВП, но с худшим разрешением.Пример 2.4. Звуковой сигналЗагрузим звуковой сигнал из файла mtlb с выборкой в 200 отсчетов (см.пример 1.7) и построим его график и две спектрограммы – непрерывного идискретного ВП:function ss_cdload mtlb; v = mtlb(1:200)', lv = length(v);subplot(311), plot(v); title('Звуковой сигнал');set(gca, 'Xlim',[0 200]); [c,l] = wavedec(v,6,'sym2');cfd = zeros(6,lv); subplot(312), ccfs = cwt(v,4:127,'sym2','plot');title('Непрерывное преобразование'), colormap(pink(32));ylabel('Временной масштаб,a'); xlabel('Временной сдвиг, b');43for m = 1:6d = detcoef(c,l,m); d = d(ones(1,2^m),:);cfd(m,:) = wkeep(d(:)',lv);endcfd=cfd(:); I = find(abs(cfd)<sqrt(eps));cfd(I) = zeros(size(I)); cfd = reshape(cfd,6,lv);subplot(313), colormap(pink(32));img = image(flipud(wcodemat(cfd,64,'row')));set(get(img,'parent'), 'YtickLabel',[]);title('Дискретное преобразование');ylabel('Уровень'); xlabel('Временной сдвиг, b');endПолученные в результате выполнения программы графики приведены нарис.
2.9.Рис. 2.9Очевидно, что все мельчайшие детали сложной временной зависимостиS (t ) отчетливо просматриваются на спектрограмме как непрерывного, так идискретного ВП. Однако последнее выполняется значительно быстрее, хотя подетальности представления уступает непрерывному ВП.Примечание. Обычное дискретное ВП (DWT) осуществляется исходя изпредположения нестационарности сигнала. Если сигнал стационарный, то вэтом частном случае можно использовать стационарное ВП [8]. Это преобразование, применяемое для очистки сигнала от шума, имеет ряд форм записи(см. прил. 2.6).442.4. Быстрое вейвлет-преобразованиеПри исследовании сигналов полезно их представление ввиде совокупности последовательных приближений грубой (аппроксимирующей) Am (t ) и уточненной (детализирующей)Dm (t ) составляющихmS (t ) = Am ( t ) + ∑ D j ( t ) ,(2.7)j =1с последующим их уточнением итерационным методом.
Каждый шаг уточнения соответствует определенному масштабуa m (т.е. уровню m ) анализа (декомпозиции) и синтеза (реконструкции) сигнала. Такое представление каждой составляющейсигнала вейвлетами можно рассматривать как во временной, таки в частотной областях. В этом суть кратномасштабного анализа(КМА). В прил. 3 описан КМА для непрерывных сигналов.Как уже отмечалось в разд. 2.3, в практике ВП в большинстве случаев мы имеем дело с дискретными сигналами. Однакоформулы для ВП дискретных сигналов не могут быть полученыпростой дискретизацией формул диадного ВП для непрерывного сигнала.
Найдем их из предпосылок КМА.Пусть имеется непрерывный сигнал S (t ) ∈V0 . Дискретныйсигнал Sд интерпретируем как последовательность коэффициентов ak , полученную в ходе КМА сигнала S (t ) при масштабирующих функциях ϕ0 k (t ) :S (t ) = A0 (t ) = ∑ a0 k ϕ0 k (t ) ,(2.8)kгдеa0 k = ak = ( S (t ), ϕ0 k (t ))– коэффициенты аппроксимации на уровне m = 0 .По концепции КМА сигнал S (t ) декомпозируется на две составляющие (принадлежащие подпространствам V1 и W1 ):S (t ) = A1 (t ) + D1 (t ) = ∑ a1k ϕ1k (t ) + ∑ d1k ψ1k (t ) .(2.9)kkСледовательно, получены две новые последовательности a1k иd1k .
Отметим, что последовательности a1k и d1k имеют поло45винную длину по сравнению с a0 k . Далее процесс декомпозиции может быть продолжен по A1 (t ) (подпространства V2 иW2 ). Сигнал S (t ) на уровне декомпозиции m будет представлен совокупностью коэффициентов amk и d mk .Однако вычисления amk и d mk по-прежнему зависят от непрерывных базисных функций ϕ(t ) и ψ(t ) . Как показано вприл. 3, эти функции однозначно определяются коэффициентами hl :ϕ(t ) = 2∑ hl ϕ(2t − l ) ,(2.10)ψ (t ) = 2∑ (−1)l h1−l ϕ(2t − l ) = 2∑ gl ϕ(2t − l ) ,(2.11)hl = (ϕ(t ), ϕ(2t − l ) ,(2.12)gl = (−1)l h2 n −1−l ,(2.13)lllгде l = 0,1,..., lo = 2n − 1 , n – порядок вейвлета. Вейвлеты n-гопорядка существуют только на интервале длиной 2n − 1 и имеют 2n отличающихся от нуля коэффициентов hl .Из (2.10) и (2.11) можно получить следующие соотношения:amk = ( S (t ), ϕmk (t )) = ∑ hl − 2 k (ϕ(t ), ϕm −1,l (t )) = ∑ hl − 2 k al ,m−1 ,ll(2.14)d mk = ( S (t ), ψ mk (t )) = ∑ gl − 2k (ϕ(t ), ϕm −1,l (t )) = ∑ gl − 2k al ,m−1 .ll(2.15)Итерационная процедура быстрого вейвлет-анализа получила название анализа от «тонкого» к «грубому» масштабу.На практике наименьший возможный масштаб (наибольшийуровень разрешения n0 ) определяется числом N дискретныхзначений сигнала ( N = 2n0 ).
На самом «тонком» значении масштаба ( m = 0 , a = 2m = 1 ) за аппроксимирующие коэффициентыa0 k принимаются сами отсчеты Si сигнала S (t ) , т.е. a0 k = Si ,46k = i , i = 0, 1,..., N − 1 . При переходе от текущего масштаба m кследующему m + 1 число вейвлет-коэффициентов уменьшаетсяв два раза и они определяются по рекуррентным соотношениям:am +1,k = ∑ hl − 2k aml , d m +1,k = ∑ gl − 2 k aml .(2.16)llПроцесс останавливается после конечного числа уровнейm = MMAX , которое зависит от протяженности сигнала ( N ) ипорядка ( l ) фильтра hl .При восстановлении (реконструкции) сигнала по его вейвлет-коэффициентам процесс идет от крупных масштабов к мелким и на каждом шаге описывается выражениемam −1,k = ∑ (hk − 2l aml + g k − 2l d ml ) ,(2.17)lкоторое получается из соотношений (2.10) и (2.11).Число операций умножения при прямом быстром ВП (БВП)будет 2 LN , где L = 2n [29].
Столько же операций необходимои для реконструкции сигнала. Таким образом, для анализасинтеза сигнала в базисе вейвлетов необходимо выполнить4 LN операций, что не превышает ( и даже меньше) числа операций для быстрого преобразования Фурье ( N log 2 N ).Таким образом, в практических приложениях с применениемБВП используются только коэффициенты hl , сами же вейвлетыне вычисляются и в расчетах не используются.Пакет Wavelet Toolbox позволяет осуществлять КМА с использованием БВП.
При этом порядок следования коэффициентов – «дерево» коэффициентов приведено на рис. 2.10: декомпозиция сигнала – сверху-вниз и реконструкция – снизу-вверх.SA1A2D1D2m=MРис. 2.1047m=0S = {Si}m=1S = A1 + D1m=2S = A2 + D2+D1S = AM + DM+DM–1,…,+D1На рис. 2.11 приведен пример КМА, взятый из раздела демонстрационных примеров «wavedemo». Слева под сигналомпредставлены аппроксимирующие коэффициенты am , а справа– детализирующие d m ( m от 1 до 5). Очевидно, что коэффициенты аппроксимации являются грубыми копиями сигнала, а детализирующие коэффициенты выделяют локальные особенности и свойства сигнала.
Справа сверху приведен также вейвлетспектр сигнала – cfs.Функции для нахождения этих коэффициентов имеют рядформ и, в частности:A = appcoef (C , L,' wname ', M ) ,(2.18)D = det coef (C , L, M )– возвращают аппроксимирующие и детализирующие коэффициенты из структуры вейвлет-разложения [C , L] , ' wname ' –строка, содержащая имя вейвлета, уровень M должен бытьцелымчислом,таким,что1 ≤ M ≤ MMAX ,гдеMMAX = length( L) − 2 .Рис. 2.1148Структура вейвлет-разложения(2.19)[C , L] = wavedec( X , N , ' wname ')возвращает векторы вейвлет-разложения сигнала X на уровнеN ; выходная структура разложения содержит векторы C и L .Реконструкцию (восстановление) сигнала S с многоуровневой структурой разложения [C , L] осуществляет функция waverec:S = waverec(C , L,' wname ') .(2.20)2.5.
О вейвлетах для БВПБольшинство используемых вейвлетов не имеют, к сожалению, аналитического выражения. Однако из предыдущегорассмотрения следует, что для практических расчетов используются не сами вейвлеты, а их коэффициенты hl .Эти коэффициенты, однозначно определяющие отцовскийϕ(t ) и материнский ψ(t ) вейвлеты, могут быть найдены из решения уравнения (2.11).Следует отметить, что процесс определения коэффициентовhl , т.е. конструирования вейвлетов, достаточно сложен дляпользователя. Да в этом и нет особой необходимости, так какуже создано большое число вейвлетов, в том числе входящих впакет расширения Wavelet Toolbox, например, вейвлеты Добеши (dbN), Симплета (sumN), Койфлета (coifN), Хаара (haar) идр.; их подробное описание приведено в [7, 8].Особо следует отметить вейвлеты Добеши.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
