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

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

Процедура рекурсивно продолжается по достижениивершины дерева.Рис. 2.18Рис. 2.18 иллюстрирует наилучшее дерево (слева) по критерию энтропии, выведенное на экран при анализе звукового сигнала (см. пример 2.4):Function ss_treeload mtlb; x = mtlb(1:200); wpt = wpdec(x,3,'db1');wpt = wpsplt(wpt, [3 0]); plot(wpt)bst = besttree(wpt); plot(bst);Это дерево действительно короче полного. Справа на рисункедан исходный сигнал – временная зависимость в узле (0,0).

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

Известны и другие алгоритмы ВП с использованием вейвлет-пакетов [1, 3, 7, 8].Пакетные вейвлет-алгоритмы встроены в Wavelet Toolbox.Функции одно– и двумерного пакетного вейвлет-анализа и синтеза, вычисления энтропии, определения наилучшего дерева покритерию энтропии и другим приведены в прил.

2 (П.2.5).572.8. Удаление шумови компрессия сигналовТрадиционно для решения этих задач применяется известный из практики фильтрации метод подавления высокочастотных составляющих спектра. Этот метод был использован впримерах 2.2 и 2.3.Кроме того, с использованием вейвлетов есть еще один метод – ограничение уровня детализирующих коэффициентов.

Задав определенный порог для их уровня и «отсекая» коэффициенты ниже этого порога, можно значительно снизить уровеньшума и сжать сигнал или изображение. Это равносильно заданию оптимального пути по дереву ВП (см. п. 2.7). Возможныразличные типы порогов ограничения: мягкий или гибкий итвердый или жесткий (рис. 2.19). При этом устанавливаютсяразличные правила выбора порога: адаптивный порог, эвристический, минимаксный и др.Но самое главное состоит в том, что пороговый уровеньможно устанавливать для каждого коэффициента отдельно. Этопозволяет строить адаптивные к изменениям сигнала (изображения) способы очистки от шума и компрессии.Подробно познакомиться с этими проблемами и приобрестиопределенные навыки читатель может с помощью демонстрационных примеров, работая с подразделами De-noise и Compression интерфейса GUI.Для практических применений Wavelet Toolbox имеет рядсоответствующих функций (см.

прил. П.2.7). Использование некоторых из них проиллюстрировано в приводимых ниже примерах.Пример 2.5. Установка мягкого или жесткого пороговОсуществляется с помощью функции y = wtresh( X , SORN , T ) . Она возвращает жесткий ythard ( SORN = ' h ') или мягкий ytsoft ( SORN = ' s ') порог( thersold ) T для входного вектора или матрицы X .y = linspace(-1,1,100); thr = 0.4;ythard = wthresh(y,'h',thr); ytsoft = wthresh(y,'s',thr);subplot(131), plot(y); title('No thersold')subplot(132), plot(ythard); title('thard thersold')subplot(133), plot(ytsoft); title('tsoft thersold')На рис.

2.19 приведен вид зависимости y = f ( x) , полученной с помощьюфункции y = wtresh( X , SORN , T ) , как при отсутствии порога, так и при жестком и мягком порогах.58Рис. 2.19Пример 2.6. Очистка от шума тестового сигнала с использованиемфункции wdencmpВ памяти компьютера имеется шесть тестовых зашумленных сигналов. Загрузим один из них из файла noismima и очистим его с помощью функцииудаления шума и сжатия wdencmp (прил.

П.2.7). Функция [xd, cxd, lxd, perf0,perfl2] = wdencmp ('lvd',c,l,wname,lev,thr,'h') возвращает очищенный и сжатыйсигнал XC , полученный из исходного сигнала X (как одномерного, так идвумерного) с использованием глобального порога THR . Выходные аргументы [CXD, LXD] – структура вейвлет-разложения вектора XC . PERFO иPERFL 2 –L2 – нормы восстановления и сжатия в процентах.PERFL 2 = 100 * ( norm(C × C ) / norm(C ))2 , где norm – норма вектора; для од22номерного сигнала PERFL 2 = 100 XC / X . N – уровень вейвлет-разложения.

SORH ( ' s ' или ' h ' ) – установка гибкого или жесткого порога.Процедура удаления шума и сжатия включает в себя три шага:1) Разложение для уровня lev. Выполняется функцией wavedec(x, lev, wname);при этом выбираются тип вейвлета и уровень декомпозиции.2) Детализация. Выбирается определенный порог для детализирующих коэффициентов. Например, используем функцию [thr,nkeep] = wdcbm(c,l,alpha,m),возвращающую порог thr относительно установленного уровня и число сохраненных коэффициентов nkeep. Параметр alpha обычно устанавливается равным 1.5 для сжатия и 3.0 для удаления шума.3) Вейвлет-восстановление.

Выполняется функцией wdencmp(.).Ниже дан листинг загрузки сигнала, его очистка от шума с компрессией ипостроение исходного и очищенного сигналов (рис. 2.20 ):function noismimaload noismima; x = noismima;wname = 'db4'; lev = 5; [c,l] = wavedec(x,lev,wname), alpha = 2; m = 2*l(1);59[thr,nkeep] = wdcbm(c,l,alpha,m);[xd,cxd,lxd,perf0,perfl2] = wdencmp('lvd',c,l,wname,lev,thr,'h');subplot(211), plot(x), title('Original signal'); axis([0,500,-10,10])subplot(212), plot(xd), title('Compressed signal'); axis([0,500,-10,10])xlab1 = ['2-norm rec.: ',num2str(perfl2)];xlab2 = ['% -– zero cfs: ',num2str(perf0),'%'];xlabel([xlab1 xlab2])endРис. 2.20С примером желательно поэкспериментировать, меняя вид тестового сигнала, уровень декомпозиции lev, параметры alpha и m, типы порогов и вейвлетов.Пример 2.7. Очистка от шума бигармонического сигналаВид сигнала и его параметры те же, что и в примере 2.2.1) Наилучшее дерево.

Как уже отмечалось, пакетные вейвлет-алгоритмымогут быть использованы для очистки от шума и сжатия сигнала. Примерприменения функции besttree (наилучшее дерево по критерию энтропии)представлен ниже:function binar_treet = 0:0.000001:0.001024;A1 = 1; A2 = 1; F1 = 10000; F2 = 2*F1; a = 90; b = 90;s1(1:200) = 0; a1 = a*0.0174533; a2 = b*0.0174533;t2 = 0.0002:0.000001:0.0008;s2 = A1*sin(2*pi*F1*t2-a1) + A2*sin(2*pi*F2*t2-a2);s3(1:224) = 0; s = [s1 s2 s3];randn('state',0); g = 0.5; n = g*randn(size(t)); x = s + n;60wpt = wpdec(x,3,'db4');wpt = wpsplt(wpt, [3 0]); plot(wpt)bst = besttree(wpt); plot(bst); endНаилучшее дерево показало на рис.

2.21, а слева; оно действительно короче полного дерева. Справа представлен исходный сигнал – временная зависимость в узле (0,0). Для получения временной зависимости в нужном узле надоустановить на него мышь и щелкнуть её левой клавишей. Рекомендуется просмотреть диаграммы сначала в первых узлах (1,0) и (1,1) дерева, чтобы «почувствовать» разделение сигнала на низкочастотную и высокочастотную составляющие. Затем, перемещаясь по дереву вниз, можно наблюдать формуотфильтрованного НЧ-сигнала (левые ветви) и подавляемой ВЧ составляющей.

На рис. 2.21, б, в даны временные диаграммы в узлах (2,0) и (4,0), т.е. влевых ветвях дерева. Это декомпозиция (аппроксимация) сигнала второго ичетвертого уровней. Очевидна очистка аппроксимированного сигнала и егосжатие (соответственно в 4 и 16 раз).(0,0)(1,1)(1,0)(2,0)(3,0)(2,1) (2,2) (2,3)(3,1)(4,0) (4,1)aвбРис. 2.21Для восстановления (реконструкции) очищенного сигнала (в первоначальном масштабе времени) необходимо выполнить обратное ДВП, отсекая при этомвысокочастотные компоненты, т.е. осуществляя пороговую обработку.612) Глобальный порог. Для удаления шума при реконструкции сигналаиспользуем функцию wpdencmp (см. прил. П.2.6) с применением глобальногопорога THR , который описывается функцией wpbmpen и получается по«штрафному» методу Бирге-Массарта [13]:function binar_de_noiset = 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;wname = 'db4'; lev = 4; tree = wpdec(x,lev,wname);det1 = wpcoef(tree,2); sigma = median(abs(det1))/0.6745;alpha = 2; thr = wpbmpen(tree,sigma,alpha);keepapp = 1; xd = wpdencmp(tree,'s','nobest',thr,keepapp);subplot(311), plot(s), title('Исходный сигнал'); axis([0,1000,-3,3])subplot(312), plot(x), title('Зашумленный сигнал'); axis([0,1000,-3,3])subplot(313), plot(xd), title('Очищенный сигнал'); axis([0,1000,-3,3])endГрафики исходного и очищенного сигналов приведены на рис.

2.22.Рис. 2.2262С этим примером можно также поэкспериментировать, меняя среднеквадратическое значение шума σ , уровень декомпозиции lev, параметры alpha иsigma, тип и/или порядок вейвлета.В табл. 2.2 в качестве примера приведены результаты исследований влияния уровня декомпозиции lev на среднеквадратическое значение шума σ yпосле ВП, на погрешности Δ ψ и σψ измерения фазового сдвига ψ% (см. пример 2.2).Т а б л и ц а 2.212345σy , В0.3070.2080.1280.1080.087Δ ψ , град–5.346–5.050–4.901–4.657–9.449σψ , град14.8114.7914.0913.7818.25levОчевидно, что существует оптимальное значение lev ( levopt = 4 ), при котором отклонение Δ ψ и среднеквадратическая ошибка σψ в измерении фазового сдвига эхосигнала минимальны.

При lev > levopt наряду с дальнейшимподавлением шума происходит искажение формы бигармонического импульсного сигнала, что приводит к росту Δ ψ и σψ .Пример 2.8. Очистка от шума звукового сигналаРассмотрим случай использования функции wden(.). Эта функция осуществляет автоматическое одномерное подавление шума (см. прил. П.2.7). Правило выбора порога определяется аргументами ' SORH ' (гибкий ' s ' или жесткий ' h ' порог) и TPTR : ' rigrsure ' – использует алгоритм Штейнанесмещенной оценки риска, ' heursure ' – эвристический вариант предыдущегоалгоритма, ' sqtwo log' – инверсный порог, 'min i max i ' – минимаксный порог.function mtlb1_noiseload mtlb; s = mtlb(501:1000);subplot(411), plot(s); title('Исходный звуковой сигнал');axis([0,500,-3,3]); t=(0:0.000001:0.000499);randn('state',0); g = 0.6; n = g*randn(size(t))'; x = s + n;subplot(412), plot(x); title('Зашумленный звуковой сигнал');axis([0,500,-3,3]); lev = 3;xdm = wden(x, 'minimaxi', 's', 'sln', lev, 'db4');subplot (413), plot(xdm), title('Minimax'); axis([0,500,-3,3])xdr = wden(x, 'rigrsure', 'h', 'sln', lev, 'db4');subplot (414), plot(xdr), title('Rigrsure'); axis([0,500,-3,3]); endРис.