А.Н. Яковлев - Введение в вейвлет преобразования (1275343), страница 9
Текст из файла (страница 9)
Для практических применений Wavelet Toolbox имеет ряд соответствующих функций (см. прил. П.2.7). Использование некоторых из них проиллюстрировано в приводимых ниже примерах.
Пример 2.5. Установка мягкого или жесткого порогов
Осуществляется с помощью функции . Она возвращает жесткий
или мягкий
порог (
)
для входного вектора или матрицы
.
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 приведен вид зависимости , полученной с помощью функции
, как при отсутствии порога, так и при жестком и мягком порогах.
Рис. 2.19
Пример 2.6. Очистка от шума тестового сигнала с использованием функции wdencmp
В памяти компьютера имеется шесть тестовых зашумленных сигналов. Загрузим один из них из файла noismima и очистим его с помощью функции удаления шума и сжатия (прил. П.2.7). Функция [xd, cxd, lxd, perf0, perfl2] = wdencmp ('lvd',c,l,wname,lev,thr,'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 noismima
load noismima; x = noismima;
wname = 'db4'; lev = 5; [c,l] = wavedec(x,lev,wname), alpha = 2; m = 2*l(1);
[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) Наилучшее дерево. Как уже отмечалось, пакетные вейвлет-алгоритмы могут быть использованы для очистки от шума и сжатия сигнала. Пример применения функции (наилучшее дерево по критерию энтропии) представлен ниже:
function binar_tree
t = 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;
wpt = 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 раз).
Рис. 2.21
Для восстановления (реконструкции) очищенного сигнала (в первоначальном масштабе времени) необходимо выполнить обратное ДВП, отсекая при этом высокочастотные компоненты, т.е. осуществляя пороговую обработку.
2) Глобальный порог. Для удаления шума при реконструкции сигнала используем функцию (см. прил. П.2.6) с применением глобального порога
, который описывается функцией
и получается по «штрафному» методу Бирге-Массарта [13]:
function binar_de_noise
t = 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.22
С этим примером можно также поэкспериментировать, меняя среднеквадратическое значение шума , уровень декомпозиции lev, параметры alpha и sigma, тип и/или порядок вейвлета.
В табл. 2.2 в качестве примера приведены результаты исследований влияния уровня декомпозиции lev на среднеквадратическое значение шума после ВП, на погрешности
и
измерения фазового сдвига
(см. пример 2.2).
Т а б л и ц а 2.2
1 | 2 | 3 | 4 | 5 | |
0.307 | 0.208 | 0.128 | 0.108 | 0.087 | |
–5.346 | –5.050 | –4.901 | –4.657 | –9.449 | |
14.81 | 14.79 | 14.09 | 13.78 | 18.25 |
Очевидно, что существует оптимальное значение (
), при котором отклонение
и среднеквадратическая ошибка
в измерении фазового сдвига эхосигнала минимальны. При
наряду с дальнейшим подавлением шума происходит искажение формы бигармонического импульсного сигнала, что приводит к росту
и
.
Пример 2.8. Очистка от шума звукового сигнала
Рассмотрим случай использования функции wden(.). Эта функция осуществляет автоматическое одномерное подавление шума (см. прил. П.2.7). Правило выбора порога определяется аргументами (гибкий
или жесткий
порог) и
:
– использует алгоритм Штейна несмещенной оценки риска,
– эвристический вариант предыдущего алгоритма,
– инверсный порог,
– минимаксный порог.
function mtlb1_noise
load 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
Рис. 2.23 дает наглядное представление о результатах очистки с помощью порогов minimax и rigrsure. Рекомендуется поэкспериментировать с раз-
Рис. 2.23
личными: среднеквадратичскими значениями шума, уровнями декомпозиции
, типами порога (‘s’ или ‘h’) и правилами выбора порога (
,
,
,
).
Глава 3