Тема 5: Функции вейвлет-преобразований В MATLAB
При изучении наук примеры не менее поучительны, нежели правила.
Исаак Ньютон (1643-1727). Английский ученый, лорд, XVII-XVII в.
В теории вейвлетов есть правило, которое поучительней самой теории - попробуй, посмотри, и делай выводы.
Евгений Прокопчук. Иркутский геофизик Уральской школы, казак, XX в.
Содержание: Введение. 5.1. Непрерывное одномерное вейвлет-преобразование. Основные функции. Интерфейс GUI. 5.2. Дискретное одномерное вейвлетное преобразование. Многоуровневое вейвлет-разложение. Обратное многоуровневое вейвлетное преобразование. Обратное одноуровневое восстановление сигнала. Функции коэффициентов аппроксимации. Задание граничных условий вейвлет-преобразования. Функции детализирующих коэффициентов вейвлетного преобразования. Вейвлет-преобразование одного уровня. Обратное вейвлет-преобразование одного уровня. Прямое восстановление сигнала из коэффициентов. Интерфейс GUI. 5.3. Дискретное двумерное вейвлет-преобразование. Функции вейвлет-преобразования. Интерфейс GUI. 5.4. Дискретное стационарное вейвлет-преобразование. Одномерное вейвлет-преобразование. Обратное одномерное вейвлет-преобразование. Интерфейс GUI. Двумерное стационарное вейвлет-преобразование. Обратное двумерное стационарное вейвлет-преобразование. Интерфейс GUI. Литература.
Введение
Пакет расширения Wavelet Toolbox системы Matlab позволяет использовать вейвлетный анализ и преобразование данных в самых различных областях науки и техники.
Программное обеспечение пакета позволяет выполнять вейвлет-преобразования как в командном режиме (и готовить специализированные программы), так и в диалоговом режиме по интерфейсу GUI (включение командой ""wavemenu"" или из окна редактора, Wavelet Toolbox ® Main Menu).
Пакет имеет демонстрационные примеры вейвлетных преобразований, окно которых включаются командой "wavedemo". Первой кнопкой окна (Command line mode) включается довольно обширное меню примеров работы в командном режиме с одно - и двумерными вейвлетами всех типов в обычном и в пакетном исполнении (непрерывные и дискретные вейвлет-преобразования с декомпозицией и реконструкцией сигналов, сжатие сигналов, очистка от шумов и пр.). Слайды примеров сопровождаются соответствующими листингами программных фрагментов, которые можно переносить в буфер (при нажатии клавиш Shift + Delete) и затем использовать в командной строке Matlab.
Аналогично второй и третьей кнопками включается доступ к демонстрационным примерам в интерфейсе GUI.
5.1. Непрерывное одномерное преобразование /3/.
Основные функции. Непрерывное одномерное вейвлет-преобразование (НВП-1D) уже само по себе, без реконструкции сигналов, используется для анализа формы сигналов и выявления их локальных особенностей. Преобразование выполняется функцией cwt в следующих форматах:
● C = cwt (S, SCALES, "wname") – возвращает коэффициенты "c" прямого НВП вещественного или комплексного сигнала S вейвлетом "wname" в шкале масштабирования SCALES. Стандартное задание SCALES = начало: шаг: конец (по значениям коэффициента масштабирования "a").
● C = cwt (S, SCALES, "wname", "plot") – то же, плюс строит график коэффициентов.
● C = cwt (S, SCALES, "wname", plotmode) – то же, с заданием для построения графика коэффициентов следующих настроек цвета plotmode:
- "lvl" – окраска шаг за шагом,
- "glb" – окраска с учетом всех коэффициентов,
- "abslvl" или "lvlabs" – окраска шаг за шагом с использованием абсолютных значений коэффициентов,
- "absglb" или "glbabs" - окраска с масштабированием и с использованием абсолютных значений коэффициентов.
● C = cwt (S, SCALES,"wname", plotmode, xlim) – то же, с дополнительной настройкой цвета xlim = , где значениями в квадратных скобках устанавливаются номера точек векторов коэффициентов С (значения сдвига bmin и bmax), по интервалу которых определяются значения Cmin и Cmax интервала изменения цветовой окраски.
§ Пример непрерывного вейвлет-преобразования синусоиды двумя типами вейвлетов и с двумя видами окраски результатов (рис. 5.1.1).
t=linspace(-3,3,2048); s=sin(t).^17; subplot(331); plot(t, s);
subplot(334); = mexihat(-4,4,100); plot(X, psi);
subplot(337); = wavefun("bior1.5",8); plot(X, psi);
subplot(335); c1=cwt(s,1:32, "mexh","lvl",);
subplot(336); c2=cwt(s,1:32, "bior1.5","lvl",);
subplot(338);c1=cwt(s,1:32, "mexh","abslvl",);
subplot(339);c2=cwt(s,1:32, "bior1.5","abslvl",);
subplot(332); plot(t, c1(10,:)); subplot(333); plot(t, c2(10,:));
Как видно на рисунке, симметричность или несимметричность (нечетность) функции вейвлета существенно влияет на форму вейвлет-спектра, равно как и тип закраски, что позволяет целенаправленно выявлять особенности локальных сигналов, например – точки перегибов сигнальных функций, которые отмечаются вертикальными полосами с вейвлетом "mtxh" при plotmode=abslvl.
§ Пример НВП гамма-каротажной кривой (рис. 5.1.2).
A=dlmread("c:\MATLAB6p1\work\Zag3f. prn"," ",120);
xn=1; xk=2000; xd=xk-xn; GK=A(,5); da=1; dk=32;
subplot(311); plot(GK); grid; axis();
subplot(312); c1=cwt(GK,1:da:dk, "mexh","abslvl",);
subplot(313); c1=cwt(GK,1:da:dk, "sym4","abslvl",);
На рис. 5.1.2 приведены вейвлет-спектры гамма-каротажной кривой. В силу природы ядерных излучений диаграммы ГК всегда сильно осложнены статистическими флюктуациями, на фоне которых достаточно трудно выделять границы пластов горных пород с разной гамма-активностью. На вейвлет-спектрах с plotmode=abslvl эти границы фиксируются по перегибам функций ГК вертикальными полосами минимумов (темные цветовые полосы), что позволяет по масштабным сечениям среднего и верхнего уровня детализировать сигнальную функцию на локальные блоки.
Интерфейс GUI удобен для анализа данных в диалоговом режиме. На рис. 5.1.3 приведено окно, которое включается из "Wavelet Toolbox Main Menu" кнопкой "Continuous Wavelet 1-D". Демонстрационные сигналы загружаются в окно из меню File ® Example Analysis, устанавливается тип вейвлета, параметры анализа и нажимается кнопка "Analyze", после чего в графической части окна появляется сигнал и результаты его разложения в трех представлениях (полное, сечение по среднему уровню разложения и линии локальных максимумов). Графическое представление можно изменять ниже расположенными кнопками и переключателями, а также используя типовые возможности оконного меню. Коэффициенты разложения можно записать на диск в виде mat-файла (File ® Save Coefficients) и считать затем в рабочую область Matlab для детального изучения.
При выполнении непрерывного разложения комплексными вейвлетами в "Wavelet Toolbox Main Menu" используется кнопка "Complex Continuous Wavelet 1-D".
При обработке данных, записанных в других (не. mat) форматах, следует сначала перевести данные в mat-формат, что можно выполнить из основного окна Matlab (Файл ® Импорт данных), или из окна команд. При обработке одномерных сигналов второе предпочтительнее, так как одновременно дает возможность подготовить для GUI однострочные векторные массивы. Ниже приведен пример считывания каротажных геофизических данных из файла las-формата, предварительно переименованного в prn-формат для использования функции dlmread, и перевода столбцов 1 и 5 считанного массива (сетка глубин и диаграмма ГК) в строковые векторы с последующей записью в файлы mat-формата.
fprn="c:\MATLAB6p1\work\MainData\Заг3\Zag3f. prn";
fn=72; A=dlmread(fprn," ",fn); rows=size(A,1); cols=size(A,2);
xn=1; xk=rows; Zag3f_9_2473=A(xn:xk,1:5);
DEPT=A(xn:xk,1); GK=A(xn:xk,5);
save "c:\MATLAB6p1\work\MainData\Заг3\gk3f. mat" GK;
save "c:\MATLAB6p1\work\MainData\Заг3\dept3f. mat" DEPT;
Запись файлов в mat-формате может выполняться и непосредственно из окна рабочей области (клик правой кнопкой мыши на выбранном для записи массиве ® "сохранить выбранное как...).
5.2. Дискретное одномерное преобразование /3/.
Главным достоинством дискретного вейвлет-преобразования (ДВП) является возможность быстрого преобразования (БВП) с пирамидальным алгоритмом вычислений, что позволяет выполнять анализ больших выборок данных. Однако возможности БВП реализуются не для всех типов вейвлетов. Тем не менее, при обработке данных БВП используется весьма интенсивно и в пакете Wavelet Toolbox представлено большим количеством специальных функций.
Многоуровневое вейвлет-разложение сигналов (декомпозиция) выполняется функцией wavedec, которая используется в двух формах:
● = wavedec (S, N, "wname").
● = wavedec (S, N, LD, HD).
Функция возвращает векторы разложения сигнала S на уровне N с использованием вейвлета "wname" или его низкочастотного (LD) и высокочастотного (HD) фильтров декомпозиции. N должно быть целым числом и определяет длину вектора L (length(L)=N+2). Значение N при К-точках сигнала должно быть не более K/2 > 2N ≤ K. Состав векторов С и L поясняет нижеследующий пример.
§ Пример вычисления коэффициентов декомпозиции сигнала (рис. 5.2.1).
dwtmode("save","per");
load sumsin; S=sumsin(1:200); =wavedec(S,3,"db4");
subplot(131); plot(S); grid; title("signal S");
subplot(132); plot(C); grid; title("signal C");
subplot(133); plot(L); grid; title("signal L");
Первой строкой в приведенном примере стоит включение метода задания граничных условий свертки. В данном случае включен метод "per" - периодизации продолжения сигналов. Возможные методы задания граничных условий будут рассмотрены ниже.
Как можно видеть на рис. 5.2.1, при N=3 и K=200 точек сигнала, вектор С представляет собой массив значений аппроксимирующих коэффициентов декомпозиции уровня N=3 (отсчеты 1-25), и массивы детализирующих коэффициентов декомпозиции уровней N=3 (отсчеты 26-50), (N-1)=2 (отсчеты 51-100) и (N-2)=1 (отсчеты 101-200). Номера стыков коэффициентов фиксируются вектором L, в L(N+2) – число отсчетов сигнала. По такому же принципу вычисляются коэффициенты любого уровня декомпозиции, до максимального значения N, при котором значение 2N не превышает количества точек сигнала.
Обратное многоуровневое преобразование (реконструкция сигнала S по структуре ) выполняется функцией waverec:
● S = waverec (C, L, "wname").
● S = waverec (C, L, LD, HD).
§ Пример декомпозиции и реконструкции сигнала (рис. 5.2.2).
load sumsin; S=sumsin(1:200); =wavedec(S,6,"db4");
s=waverec(C, L,"db4"); err=norm(S-s)
subplot(131); plot(s); grid; title("signal s");
C(L(7)+1:L(8))=0; s1=waverec(C, L,"db4");
subplot(132); plot(s1); grid; title("signal s1");
C(L(3)+1:L(8))=0; s2=waverec(C, L,"db4");
subplot(133); plot(s2); grid; title("signal s2"));
На рис. 5.2.2 представлены три реконструкции сигнала: сигнал s полной реконструкции, сигнал s1 – реконструкции с обнулением детализирующих коэффициентов 1-го уровня, и сигнал s2 – с обнулением всех детализирующих коэффициентов. Вычисленная метрика полной реконструкции имеет порядок -12 степени, т. е. по существу представляет собой точность машинных вычислений.
Восстановление одиночной ветви коэффициентов из структуры может быть выполнено функцией wrcoef:
● X = wrcoef ("type",C, L,"wname",N) – возвращает вектор восстановленных коэффициентов уровня N, приведенных к исходным координатам сигнала. Уровень N может быть любым, но не более уровня Nmax = length(L)-2 структуры . Параметром "type" задается тип коэффициентов: "a" – аппроксимации, "d" – детализации.
● X = wrcoef ("type",C, L,LR, HR, N) – то же, с использованием фильтров реконструкции данного вейвлета.
Без указания параметра N восстановление производится с максимального уровня.
§ Пример восстановления одиночной ветви коэффициентов аппроксимации с использованием двух типов вейвлетов (рис. 5.2.3).
load sumsin; S=sumsin; subplot(311); plot(S); grid; axis();
Wavedec(S,5,"haar"); =wavedec(S,5,"db10");
subplot(323); plot(C1); grid; axis();
subplot(324); plot(C2); grid; axis();
a5h=wrcoef("a",C1,L1,"haar",5); a5d=wrcoef("a",C2,L2,"db10",5);
subplot(325); plot(a5h); grid; axis();
subplot(326); plot(a5d); grid; axis();
На рис. 5.2.3 достаточно ясно видно, что детализационные коэффициенты С2 отсутствуют на уровнях 2 и 5. Восстановление детализационных коэффициентов на уровнях 1, 3, 4 и сумма коэффициентов на уровнях 3 и 4 позволяют выделить высокочастотные составляющие сигнала, показанные на рис. 5.2.4.
§ Пример восстановления одиночных ветвей коэффициентов детализации вейвлета "db10" (продолжение предыдущего примера). Для наглядности результаты на рис. 5.2.4 приведены только в ограниченной части интервала восстановления.
d1=wrcoef("d",C2,L2,"db10",1); d3=wrcoef("d",C2,L2,"db10",3);
d4=wrcoef("d",C2,L2,"db10",4); a=d3+d4; subplot(141); plot(d1); grid;
subplot(142); plot(d3); grid; subplot(143); plot(d4); grid; subplot(144); plot(a); grid;
Обратное одноуровневое восстановление (уровня n с уровня n+1) выполняется функцией upwlew в форме:
● = upwlev (C, L, "wname"), где C, L – структура разложения сигнала на уровне length(L)-2 = n+1, сА – вектор аппроксимационных коэффициентов уровня n+1.
● = upwlev (C, L, LR, HR) – то же, с использованием низкочастотного (LR) и высокочастотного (HR) фильтров реконструкции данного вейвлета.
§ Пример декомпозиции и обратного одноуровневого восстановления (рис. 5.2.5).
load sumsin; S=sumsin(1:200); S=dwtmode("per","nodisp");
Wavedec(S,4,"db4"); = upwlev(C, L, "db4");
subplot(221); plot(S); grid; subplot(222); plot(C); grid;
subplot(223); plot(cA); grid; subplot(224); plot(nC); grid;
Из структуры вейвлет-разложения функциями appcoef и detcoef извлекаются соответственно коэффициенты аппроксимации и детализации.
Функции коэффициентов аппроксимации :
● А = appcoef (C, L, "wname", N) – возвращает коэффициенты аппроксимации сигнала вейвлетом "wname" на уровне N без изменения масштаба их представления.
● А = appcoef (C, L, "wname") – возвращает коэффициенты аппроксимации сигнала на последнем уровне Nmax = length(L)-2.
Вместо имени вейвлета могут быть указаны непосредственно фильтры реконструкции, т. е.:
● А = appcoef (C, L, LR, HR, N) , где LR и HR – векторы низкочастотного и высокочастотного фильтров данного вейвлета.
§ Пример выделения сигнала аппроксимации (рис. 5.2.6).
load sumsin; S=sumsin(1:800); =wavedec(S,4,"sym8"); A=appcoef(C, L,"sym8",4);
subplot(131); plot(S); grid; subplot(132); plot(C); grid; subplot(133); plot(A); grid;
На рис. 5.2.6 приведен пример применения функции appcoef. По существу, эта функция выполняет извлечение коэффициентов аппроксимации соответствующего уровня из комбинированного массива всех коэффициентов разложения сигнала. В данном примере на 4-ом уровне декомпозиции выделяется практически чистая низкочастотная составляющая сигнала (без учета начальной части искажения за счет граничных условий).
Задание граничных условий. Искажения в начальной части функций коэффициентов (см. рис. 5.2.6А) определяются граничными условиями, которые устанавливаются для данного входного сигнала (продление сигнала влево и вправо от границ задания для выполнения операции свертки) специальной функцией dwtmode:
● S = dwtmode ("mode","nodisp"), где параметр моды (метода) может быть следующим:
"sym" – симметричное дополнение (используется по умолчанию),
"zpd" – дополнение нулями,
"sp0" или "sp1" – гладкое дополнение нулевого или первого порядка,
"ppd" – периодическое дополнение,
"per" – метод периодизации.
Параметр "nodisp" не обязателен, при его отсутствии на дисплей выводится информация о смене типа моды для данного сигнала. Метод по умолчанию загружается из файла dwtmode. def (или, при его отсутствии, из dwtmode. cfg).
Изменение коэффициентов А, приведенных на рис. 5.2.6, в зависимости от начальных условий можно видеть на рис. 5.2.7.
Метод продления по умолчанию может изменяться командой:
● dwtmode ("save" ,"mode").
Информация о текущем методе по умолчанию выводится на экран командой >> dwtmode, об установленном методе для данного сигнала S: >> S=dwtmode.
Функции детализирующих коэффициентов :
● D = detcoef (C, L, N) – возвращает детализирующие коэффициенты на уровне N из структуры разложения . Значение N может быть любым в интервале 1≤ N ≤ length(L)-2.
● D = detcoef (C, L) – то же, на последнем уровне Nmax = length(L)-2.
§ Пример выделения сигнала детализации 1-го уровня (рис. 5.2.8).
load sumsin; S=sumsin; =wavedec(S,3,"db3"); D1=detcoef(C, L,1);
subplot(131); plot(S); grid; subplot(132); plot(C); grid; subplot(133); plot(D1); grid;
Детализирующие коэффициенты могут быть записаны в массив раздельных массивов коэффициентов функциями:
● DJ = detcoef (C, L, N, "cells"), где N – вектор целых чисел из интервала 1≤ N(j) ≤ Nmax, при этом DJ образует массив ячеек DJ{j}, каждая из которых содержит детализирующие коэффициенты соответствующего N-уровня. При записываются коэффициенты всех уровней.
● DJ = detcoef (C, L,) – то же, с непосредственным указанием вектора чисел N перечислением в , или последовательным перечислением типа .
● DJ = detcoef (C, L,"cells") эквивалентна функции detcoef (C, L,).
● = detcoef (C, L,) – записывает детализирующие коэффициенты в отдельные массивы.
Вейвлет-преобразование одного уровня вейвлетом "wname" вектора (сигнала) S выполняется функцией dwt:
● = dwt (S, "wname") – возвращает векторы коэффициентов аппроксимации сА и детализации cD.
● = dwt (S, LD, HD) – то же, с использованием низкочастотного LD и высокочастотного HD фильтров декомпозиции.
● = dwt (…, "mode", MODE) – то же, с заданием метода расширения.
Длина векторов cA и cD при методе расширения "per" равна половине входного вектора S, при других методах расширения равна половине суммы векторов входного сигнала и фильтра декомпозиции.
Обратное преобразование одного уровня вейвлетом "wname" (или фильтрами реконструкции) выполняется функцией idwt:
|
● S = idwt (cA, cD, "wname").
● S = idwt (cA, cD, LR, HR).
● S = idwt (……, L) – возвращает центральный блок вектора S, независимый от условий расширения.
● S = idwt (……,"mode", MODE) – возвращает вектор S с заданием метода расширения.
● Sа = idwt (cA, ,…..) – возвращает аппроксимированный вектор S.
● Sd = idwt (, cD,…..) – возвращает функцию детализации вектора S.
§ Пример декомпозиции и реконструкции сигнала на 1-м уровне (рис. 5.2.9).
randn("seed",123456789) % Set randn
s=2+kron(ones(1,10),)+((1:20).^2)/32+0.5*randn(1,20);
Dwt(s,"db4");
subplot(221); plot(cA1); grid;
subplot(222); plot(cD1); grid;
ss=idwt(cA1,cD1,"db4");
err=norm(s-ss); t=kron(1:20,1);
subplot(212); plot(t, s,t, ss); grid; xlabel(["err=",num2str(err)])
Прямое восстановление сигнала из коэффициентов уровня N задается функцией upcoef в следующем формате:
|
● Y = upcoef (O, X, "wname", N), где О = "a" для коэффициентов аппроксимации, или "d" для детализационных коэффициентов. Вместо имени вейвлета можно задавать фильтры реконструкции LR и HR. Если значение N не указывать, то восстанавливаются коэффициенты первого уровня.
● Y = upcoef (O, X, "wname", N, L) – то же, с выделением центрального блока результата размером L.
§ Пример восстановления сигнала из аппроксимационных коэффициентов разложения импульса Кронекера с разных уровней в одном масштабе (рис. 5.2.10).
S=1; N=4; es=10;
reg=upcoef("a", S, "db8",i);
ax=subplot(N,1,i); h=plot(reg(1:es)); grid;
set(ax, "xlim", ); es=es*2;
Интерфейс GUI при дискретном преобразовании используется аналогично непрерывному преобразованию.
Окно преобразования включается из "Wavelet Toolbox Main Menu" кнопкой "Wavelet 1-D" и имеет несколько большие функциональные возможности. Пример окна приведен на рис. 5.2.11 в режиме графического вывода всех коэффициентов разложения (Separate Mode). Переключателем "Display mode" режимы вывода можно изменять, детали вывода можно устанавливать в отдельном подокне (кнопка "More Display Options").
Под кнопкой "Analyze" окно имеет 4 кнопки включения окон выполнения специальных операций над результатами разложения сигнала.
Окна "Statistics" и "Histograms" предназначены для анализа и графического вывода статистических характеристик сигнала и всех коэффициентов его разложения. В окнах "Compress" и "De-noise" устанавливаются режимы компрессии (сжатия) сигналов и очистки сигналов от шумов, и выполняются эти операции. Их применение будет рассмотрено в дальнейшем.
Коэффициенты декомпозиции в структуре разложения , а также синтезированный сигнал могут быть записаны на диск в mat-файлы через меню окна File. Точно так же в окно может загружаться как сигнал, так и коэффициенты его разложения.
5.3. Дискретное двумерное преобразование.
Функции двумерного вейвлет-преобразования по своему назначению аналогичны описанным выше функциям одномерного преобразования и идентифицируются цифрой 2 в конце одноименных функций одномерного преобразования, т. е.: wavedec2, wavereg2, wrcoef2, upwlev2, appcoef2, detcoef2, upcoef2, dwt2, idwt2.
§ Пример декомпозиции и аппроксимационной реконструкции двумерного сигнала (рис. 5.3.1).
randn("seed",12345);
S(1:100,1:100)=0; S(60:70,60:70)=10;
S=S+5*randn(100,100);
Wavedec2(S,4,"db2");
a3=wrcoef2("a",C, L,"db2",3);
a4=wrcoef2("a",C, L,"db2",4);
figure(1); subplot(132); contour(a3); subplot(131); contour(S); subplot(133); contour(a4);
figure(2); subplot(131); mesh(S); subplot(132); mesh(a3); subplot(133); mesh(a4);
На рис. 5.3.1 приведен пример обработки двумерного сигнала, информационная часть которого соизмерима с мощностью шумов. Результаты обработки приведены в виде контурных карт и карт поверхностей. При двумерной обработке входной сигнал S представляет собой двумерную матрицу, структура разложения состоит из вектора коэффициентов аппроксимации и детализации C и матрицы структуры записи коэффициентов L. Вектор С состоит из последовательностей коэффициентов аппроксимации и детализации горизонтальных, вертикальных и диагональных строк разложения сигнала.
Аналогичным образом работают и другие функции двумерного преобразования.
Интерфейс GUI двумерного преобразования из "Wavelet Toolbox Main Menu" включается кнопкой "Wavelet 2-D". Окно преобразования показано на рис. 5.2.2 и по использованию во многом подобно окну одномерного преобразования.
Построение картины разложения двумерного сигнала (правый нижний график) в режиме "View mode: Square" следующее. Диагональные коэффициенты детализации первого уровня разложения находятся в правом нижнем углу графика, горизонтальные коэффициенты над ним в правом верхнем углу, вертикальные – слева, в левом нижнем углу. Левая верхняя четверть графика – коэффициенты последующих уровней, в данном случае второго и третьего, которые располагаются по аналогичному принципу. Левый верхний квадрат графика - коэффициенты аппроксимации последнего уровня разложения. В правом верхнем графике окна можно наблюдать отмеченные коэффициенты разложения (включая кнопки Operations on selected image) в исходном масштабе входного сигнала (Visualize), в увеличенной форме (Full Size), и в режиме реконструкции. Кнопками 1,2,3 и 4 Full Size можно включать на просмотр в увеличенном масштабе любой из четырех графиков окна.
Вид графической части окна в режиме View mode: Tree приведен на рис. 5.3.3 и комментариев не требует.
5.4. Дискретное стационарное вейвлет-преобразование /34/.
Дискретное преобразование dwt исходит из предпосылки нестационарности сигнала. В дискретном стационарном вейвлет-преобразовании (SWT – Stationary Wavelet Transform) сигнал рассматривается как стационарный. Использование данного преобразования предпочтительно при сжатии сигналов и очистке от шумов.
Одномерное вейвлет-преобразование выполняют функции:
|
● S = swt (X, N, "wname") – возвращает разложение сигнала Х до уровня N вейвлетом "wname". N – только положительное число, размер Х должен быть кратным 2N. Выходная матрица S состоит из векторов – строк и содержит детализирующие коэффициенты S(n,:) уровней n строках 1≤n≤ N, и аппроксимирующие коэффициенты уровня N в строке S(N+1,:). Временной (координатный) масштаб коэффициентов соответствует масштабу входного сигнала, что очень удобно для анализа.
● S = swt (X, N, LD, HD) – то же, с использованием низкочастотного LD и высокочастотного HD фильтров.
● = swt (…) – возвращаются матрицы коэффициентов аппроксимации cA и детализации cD для всех n-уровней. Коэффициенты уровней располагаются построчно с нумерацией от 1 до N.
§ Пример декомпозиции сигнала (рис. 5.4.1).
load noisbloc; s=noisbloc(640:879);
subplot(311); plot(s); axis();
sc=swt(s,3,"db1"); =swt(s,3,"db1");
subplot(312); plot(ca"); axis();
subplot(313);plot(cd(1,:)"); axis();
Обратное одномерное вейвлет-преобразование матриц S и реализуется практически с абсолютной точностью следующими функциями:
● X= iswt (S, "wname"),
● X= iswt (cA, cD, "wname"),
● X= iswt (cA(end,:), cD, "wname"),
● X= iswt (S, LR, HR),
● X= iswt (cA, cD, LR, HR),
● X= iswt (cA(end,:), cD, LR, HR).
§ Пример проверки точности восстановления сигнала.
load noisbloc; s=noisbloc; sc=swt(s,3,"db1"); =swt(s,3,"db1");
s1=iswt(sc,"db1"); s2=iswt(ca, cd, "db1"); err1=norm(s-s1); err2=norm(s-s2);
err1m=max(max(abs(s-s1))); err2m=max(max(abs(s-s2)));
Вычисления данного примера дают:
err1 = err2 = 9.6566e-014, err1m = err2m = 1.0658e-014.
Двумерное стационарное вейвлет-преобразование сигнала Х до уровня N выполняется функциями:
● S = swt2 (X, N, "wname"), ● = swt2 (X, N, "wname").
● S = swt2 (X, N, LD, HD), ● = swt2 (X, N, LD, HD).
В этих функциях N – только положительное число, размеры Х по обеим осям (size(X,1) и size(X,2)) должны быть кратными 2N. Выходная матрица S – массивы содержат коэффициенты аппроксимации на уровне n в матрице А(:,:,n), и детализационные коэффициенты на уровне n в матрицах H(:,:,n) – горизонтальные, V(:,:,n) – вертикальные, и D(:,:,n) – диагональные.
|
§ Пример двумерного разложения сигнала.
load facets; m=256; N=3;
codX=wcodemat(X, m);
subplot(221); image(codX);
Swt2(X, N,"sym4");
coda=wcodemat(ca(:,:,k),m);
codhd=wcodemat(chd(:,:,k),4*m);
codvd=wcodemat(cvd(:,:,k),4*m);
coddd=wcodemat(cdd(:,:,k),4*m);
decl=;
subplot(2,2,k+1); image(decl);
Утилита wcodemat в приведенном примере возвращает кодированные версии входных матриц в градациях значений между максимальным и минимальным значением. Числом градаций может устанавливаться световая контрастность графических отображений матриц. Полный формат утилиты
● Y = wcodemat(X, M,OPT, ABS), где М – число первых чисел (градаций) кодирования, ABS – тип кодирования (ABS=0 – с учетом знаков, ABS>0 – по абсолютным значениям), OPT – строковый параметр = "row" или "r" кодирования по строкам, = "col" или "c" кодирования по столбцам, "mat" или "m" глобального кодирования.
● Y = wcodemat(X, M,OPT) эквивалентна Y = wcodemat(X, M,OPT,1),
● Y = wcodemat(X, M) эквивалентна Y = wcodemat(X, M,"mat",1),
● Y = wcodemat(X) эквивалентна Y = wcodemat(X,16,"mat",1).
Обратное двумерное стационарное вейвлет-преобразование сигнала Х по матрицам, вычисленным утилитой swt2, выполняется функциями:
● X=iswt2 (S,"wname"), ● X=iswt2 (A, H,V, D,"…"), ● X=iswt2 (A(:,:,end),H, V,D, "…").
● X=iswt2 (S, LR, HR), ● X=iswt2 (A, H,V, D,LR, HR), ● X=iswt2 (A(:,:,end),H, V,D, LR, HR).
Сигнал, как и при одномерном восстановлении, также реконструируется с очень высокой точностью, что рекомендуется проверить самостоятельно.
Интерфейс GUI не имеет специального окна SWT. Но с использованием SWT выполняется сжатие сигналов и очистка от шумов в окне "SWT De-noising 1-D" "Wavelet Toolbox Main Menu", а, соответственно, в этом окне при необходимости можно просматривать и вейвлет-разложение сигнала.
литература
3. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. Специальный справочник. – СПб.: Питер, 2002, 608 с.
1.2.2. MatLab Wavelet Toolbox
MatLab Wavelet Toolbox – это открытый, дружественный для пользователя пакет расширения MatLab, позволяющий синтезировать всевозможные алгоритмы обработки информации - данных, сигналов и изображений - с использованием вейвлет-функций /6/. В своей работе пакет широко использует возможности системы MatLab (матричные алгоритмы вычислений, стильную и в тоже время мощную графику) для решения задач анализа (шумоподавления, расфильтровки, сжатия и восстановления): это предоставляет в распоряжение как начинающего, так и профессионального пользователя исчерпывающий набор функций для реализации собственных алгоритмов обработки данных, т.е. написания собственного m-кода, а также средства графического интерфейса (GUI). Можно сказать, пакет Wavelet Toolbox оказывается превосходным средством для решения задач обработки одно- и двумерной информации: действительно, спектр задач, решаемых с использованием пакета, настолько широк, что упоминание таких проблем, как обработка звука, статических изображений и видеокартинок, не говоря уже о передаче данных, исследовании массивов геофизических, сейсмоакустических данных, биомедицинских сигналов и изображений, будет, естественно, далеко не полным.
MatLab Wavelet Toolbox включает обширную библиотека вейвлет-функций (континуальных неортогональных вейвлетов, в том числе комплексных; ортогональных семейств функций, функций Добеши, Койфмана, а также симлетов; биортогональных вейвлетов); широкий набор вейвлет-фильтров /7/.
Основные возможности:
1) всевозможные функции для реализации континуального анализа, дискретного одноуровневого и дискретного многоуровневого анализа;
2) функции анализа и синтеза данных с использованием вейвлет-пакетов;
3) функции для решения задач аппроксимации данных, статистических распределений и т.п.;
4) функции внедрения в пакет собственных вейвлет-функций и работы с ними;
5) набор средств визуализации результатов анализа и синтеза;
6) средства GUI.
1.2.3. Вывод по аналитическому обзору
Список программных продуктов, безусловно, может быть расширен, но все же самые характерные и популярные разработки в него включены.
Однако, несмотря на множество достоинств, они имеют следующие недостатки:
1) не реализуют метод структурной индексации исходных сигналов;
2) обладают высокими требованиями к аппаратному обеспечению;
3) имеют высокую стоимость;
4) понятие вейвлета в них строго детерминировано для реализации уже разработанных алгоритмов.
Этих недостатков лишена система МАДС. Кроме того, ограничения, накладываемые вышеупомянутыми системами на структуру вейвлета, в данной работе сняты: вейвлет по своей сущности здесь ничем не отличается от сигнала. Это открывает перед нами широкое поле для экспериментов, в том числе и по изучению фрактальных свойств сигнала.
Поэтому данная разработка является востребованной в современной индустрии компьютерной обработки сигналов.
1.3. Основные требования к системе 1.3.1. Основные цели создания системы и критерии эффективности ее функционирования
Создание системы многомасштабного анализа дискретных сигналов позволит получить новые возможности по выявлению структурных особенностей сигналов, подавлению в них шумов, сжатию данных.
Для оценки эффективности работы системы МАДС можно использовать оценку размера данных до и после сжатия.
1.3.2. Функциональное назначение системы
Автоматизация процесса многомасштабного анализа дискретных сигналов подразумевает реализацию в системе определенных средств и функций. Следует выделить ряд функциональных особенностей, которыми должна обладать система МАДС:
1) осуществление вейвлет-преобразования исходных сигналов;
2) осуществление структурной индексации исходных сигналов;
3) конвертация результатов структурной индексации для получения исходного сигнала;
4) визуализация данных вейвлет-преобразования и структурной индексации для наглядного отображения их результатов.
1.3.3. Особенности системы и условия ее эксплуатации
Система МАДС предназначена для работы с текстовыми файлами, содержащими данные о различных сигналах. Таким образом, объем информации, обрабатываемый системой, может быть достаточно велик и составлять десятки мегабайт. Эти особенности накладывают ограничения на использование непроизводительных и медленных алгоритмов.
1.3.4. Требования к функциональной структуре
Построение системы многомасштабного анализа дискретных сигналов предполагает модульную структуру. Общий интерфейс и возможность доступа ко всем модулям в составе системы должна обеспечивать оболочка. Из оболочки МАДС вызываются следующие модули: подсистема вейвлет-анализа, подсистема структурной индексации, подсистема конвертации данных структурной индексации в исходный сигнал, подсистема визуализации исходного сигнала и результатов вейвлет-преобразования и структурной индексации. Обмен данными между подсистемами происходит через проект в рамках общей оболочки.
Подсистема вейвлет-анализа служит для вейвлет-преобразования исходного сигнала.
Подсистема структурной индексации предназначена для реализации методов структурной индексации исходного сигнала.
Подсистема конвертации данных структурной индексации служит для преобразования результата структурной индексации, а также для получения из него вновь исходного сигнала.
Подсистема визуализации предназначена для отображения исходного сигнала, результатов работы подсистем вейвлет-анализа, структурной индексации и конвертации данных структурной индексации в виде графического изображения.
1.3.5. Требования к техническому обеспечению
Задача обработки дискретных сигналов в системе МАДС связана с автоматическим анализом больших массивов информации. Преобразования, проводимые в системе, должны проводиться в процессе интерактивного взаимодействия с пользователем, поэтому паузы на обработку не должны превышать нескольких минут. Исходя из этого, сформулированы требования к техническим характеристикам персонального компьютера, на котором будет функционировать система. Требования сведены в табл. 1.1.
Таблица 1.1
Технические характеристики персонального компьютера
| Наименование | Значение |
| Частота процессора, МГц | от 900 |
| Объем оперативной памяти, Мб | от 128 |
| Разрешение экрана монитора | не менее 1024x768 |
Основным видом информации, получаемым в системе МАДС, является графическая информация в растровом представлении. Такой вид данных воспринимается человеком непосредственно и целостно, поэтому необходимо обеспечить средства наглядной визуализации изображений на различных этапах обработки.
1.3.7. Требования к программному обеспечению
Систему МАДС целесообразно разрабатывать для функционирования под операционной системой семейства Windows, так как ОС данного класса наиболее широко распространены в современном мире. Платформой для разработки выбрана среда для разработки приложений Microsoft Visual Studio .NET. Эта среда поддерживает язык C# и обладает при этом возможностями быстрой разработки и проектирования визуальных интерфейсов, что особенно важно при работе с графической информацией.
1.4. Основные технические решения проекта системы 1.4.1. Решение по комплексу технических средств
Как уже отмечалось в п.п. 1.3.5, для достижения удобного пользователю режима функционирования системы необходимо следующая минимальная конфигурация персонального компьютера: частота процессора 900 МГц, объем оперативной памяти 128 Мб, монитор, поддерживающий разрешение 1024x768 точек. Также желательно наличие следующих периферийных технических средств: цветной струйный принтер для вывода на печать результатов обработки изображений.
1.4.2. Описание системы программного обеспечения
Для реализации и функционирования проекта необходимо общесистемное программное обеспечение ОС Windows XP, в основе которой лежит ядро, характеризуемое 32-разрядной вычислительной архитектурой и полностью защищенной моделью памяти, что обеспечивает надежную вычислительную среду.
Разработка системы МАДС и ее подсистем будет вестись с использованием среды для разработки приложений Microsoft Visual Studio .NET. Среда разработки включает в себя высокопроизводительный 32-битный компилятор, что позволяет оптимизировать создаваемый код. Microsoft Visual Studio .NET включает обширный набор средств, которые повышают производительность труда программистов и сокращают продолжительность цикла разработки. Многофункциональная интегрированная среда разработки Microsoft Visual Studio .NET включает компилятор, удовлетворяющий стандарт ANSI/ISO, встроенный дизайнер форм, богатый набор средств для работы с компонентами, менеджер проектов и отладчик. Удобство разработки и эффективность созданных в данной среде разработки программ делают Microsoft Visual Studio .NET оптимальным выбором для построения исследовательской системы, какой является система МАДС.
2. РАЗРАБОТКА ПОДСИСТЕМЫ ВЕЙВЛЕТ-АНАЛИЗА 2.1. Описание постановки задачи вейвлет-анализа 2.1.1. Характеристика задачи
Для того чтобы функция могла называться вейвлетом, должны выполняться два условия /8/:
1) ее среднее значение (т.е. интеграл по всей прямой) равно нулю: ;
2) функция быстро убывает при .
Теперь возьмем произвольный сигнал – некоторую функцию (переменную будем называть временем), и произведем ее вейвлет-анализ при помощи вейвлета .
Результатом вейвлет-анализа этого сигнала будет функция , которая зависит уже от двух переменных – от времени и от масштаба . Для каждой пары и рецепт вычисления значения следующий:
1) растянуть вейвлет в раз по горизонтали и в раз по вертикали;
2) сдвинуть его в точку , полученный вейвлет обозначим ;
3) "усреднить" значения сигнала в окрестности точки a при помощи :
, (2.1)
где 
– вейвлет со смещением и масштабом /3/.
Но это всё в теории. На практике мы имеем дело с дискретным входным сигналом и дискретным (или дисркетизированным) вейвлетом. Соответственно, результатом дискретного вейвлет-анализа должна явиться матрица, каждую точку которой можно сопоставить конкретному значению входного сигнала и конкретному масштабу вейвлета.
Таким образом, задача вейвлет-анализа может быть разбита на несколько подзадач:
1) передискретизация сигнала, куда входит нахождение и разработка алгоритма, выполняющего масштабирование (сжатие и разжатие) вейвлета, заданного в дискретном виде (т.е. в виде, аналогичном входному сигналу);
2) перемножение сигнала и вейвлетва, т.е. расчёт одной строки матрицы результата вейвлет-анализа, соответствующей одному масштабу вейвлета;
3) собственно вейвлет-анализ, производящий последовательное масштабирование вейвлета и его перемножение с сигналом и получающий результирующую матрицу целиком.
Результат вейвлет-анализа легко визуализируется в любой цветовой шкале и может быть использован для выявления нестационарных составляющих сигнала, что крайне полезно при подборе способов фильтрования сигнала с помощью структурной индексации.
В результате построения подсистемы вейвлет-анализа система многомасштабного анализа дискретных сигналов (МАДС) дополнит свои функциональные возможности способностью выделения из исходного сигнала наиболее четких его составляющих, что должно быть учтено при дальнейшей его очистке от шумов.
2.1.2. Входная информация
Входной информацией являются текстовые файлы с расширением «.dat» (от англ. data – данные), содержащие данные исходного сигнала.
Структура входного файла «.dat»:
где – количество данных;
, – значение сигнала, целое число.
2.1.3. Выходная информация
Выходной информацией для данной задачи являются текстовые файлы с расширением «.war» (от англ. wavelet analysis result – результат вейвлет-анализа), содержащие результаты вейвлет-анализа.
Структура выходного файла «.war»:
где – ширина растра;
– высота растра;
, , – результат вейвлет-анализа, вещественное число.
2.1.4. Математическая постановка задачи 2.1.4.1. Математическое описание задачи передискретизации сигнала
Исходный и результирующий сигналы представляют собой одномерные массивы чисел.
Целью передискретизации исходного сигнала размером является получение сигнала размером по следующему закону:
где 
– индекс элемента в исходном сигнале, участвующего в вычислении -го элемента результирующего сигнала;
, – исходный сигнал;
, – передискретизированный сигнал;
– модуль (длина) вектора;
Соглашение об использовании материалов сайта
Просим использовать работы, опубликованные на сайте , исключительно в личных целях. Публикация материалов на других сайтах запрещена.
Данная работа (и все другие) доступна для скачивания совершенно бесплатно. Мысленно можете поблагодарить ее автора и коллектив сайта.
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
Подобные документы
Плоскость частота-время для анализа и сравнения частотно-временных локализационных свойств различных базисов. Понятие базисных функций. Прямое и обратное преобразование Фурье. Сущность дискретного вейвлет-преобразования и примеры функции вейвлет.
курсовая работа , добавлен 21.11.2010
Идея и возможности вейвлет-преобразования. Свойства вейвлетов: непрерывное прямое и обратное образование. Понятие и оценка преимуществ, сферы применения дискретного вейвлет-преобразования. Поиск изображений по образцу. Многомасштабное редактирование.
курсовая работа , добавлен 27.04.2011
Преобразования Фурье, представление периодической функции суммой отдельных гармонических составляющих. Использование преобразований как для непрерывных функций времени, так и для дискретных. Программа и примеры реализации алгоритмов с прореживанием.
реферат , добавлен 25.05.2010
Методика нахождения различных решений геометрических задач на построение. Выбор и применение методов геометрических преобразований: параллельного переноса, симметрии, поворота (вращения), подобия, инверсии в зависимости от формы и свойств базовой фигуры.
курсовая работа , добавлен 13.08.2011
Уравнения Фредгольма и их свойства как классический пример интегральных уравнений с постоянными пределами интегрирования, их формы и степени, порядок формирования и решения. Некоторые приложения интегральных уравнений. Общая схема метода квадратур.
курсовая работа , добавлен 25.11.2011
Определение, свойства и примеры функциональных уравнений. Основные методы их решения, доказательство некоторых теорем. Понятие группы функций, применение их при решении функциональных уравнений с несколькими переменными. Класс уравнений типа Коши.
курсовая работа , добавлен 01.10.2011
Параллельные методы решения систем линейных уравнений с ленточными матрицами. Метод "встречной прогонки". Реализация метода циклической редукции. Применение метода Гаусса к системам с пятидиагональной матрицей. Результаты численного эксперимента.
курсовая работа , добавлен 21.10.2013
Аналитические свойства интегральных преобразований. Интеграл Коши на различных кривых. Аналитическая зависимость от параметра. Существование производных всех порядков у аналитической функции. Вывод формулы Коши и формулировка следствий из данной формулы.
курсовая работа , добавлен 10.04.2011
Пакет расширения системы Matlab 7.0.1 Wavelet Toolbox 2 - одно из мощных инструментальных средств для изучения и применения вейвлет-преобразования (ВП). Этот пакет представляет обширные и уникальные возможности для работы с вейвлетами в двух режимах: командного окна и графического интерфейса пользователя.
Вейвлеты в пакете Matlab принято классифицировать по имени ученого, впервые предложившего тот или иной тип вейвлета, или по отличительным особенностям. В данном пакете содержится 15 базовых типов вейвлетов, в частности, Хаара - ‘haar’, Морле - ‘шогГ, Добеши - ‘db’, мексиканская шляпа - ‘mexh’ и т.д. Для получения справки по типу вейвлета необходимо в командном окне набрать waveinfo (‘type’), указав тип вейвлета. Например, для вейвлета типа "сомбреро" имеем:
>> waveinfо(1 mexh 1)
MEXHINFO Information on Mexican Hat wavelet.
Mexican Hat Havelet
Definition: second derivative of the Gaussian probability density function
mexh(x) = с * ехр(-х л 2/2) * (1-х л 2) where c = 2/(sqrt(3)*pi A {1/4})
Family Mexican hat
Compact support no
Support width infinite
Effective support [-5 5]
Режим командного окна
В режиме командного окна пакет системы Matlab Wavelet Toolbox предоставляет:
- инструментальные средства вейвлет-анализа и синтеза;
- множество встроенных вейвлетов;
- возможность конструирования вейвлета с заданными свойствами;
- средства обработки сигналов и изображений;
- средства для непрерывного и дискретного вейвлет-анализа;
- средства очистки сигнала от шума и другие возможности. Достаточно подробное описание приемов и средств вейвлет-
анализа приведено в .
Рассмотрим некоторые примеры использования вейвлетов в режиме командного окна.
Пример 6.11. Вычислить значения вейвлета типа «сомбреро» в диапазоне [-6; 6J.
Функция =mexihat(lb,uto,n); plot(x,psi)
получим итоговый результат, показанный на рис. 6.31. Здесь представлена вейвлет-функция типа "сомбреро", вычисленная в диапазоне [-6; 6J.
Рис.6.31. График вейвлет-функции типа «сомбреро» Непрерывное вейвлет-преобразование
Этот вид преобразования лежит в основе применения
вейвлетов в технике обработки сигналов. НПВ реализуется функцией cwt, синтаксис которой имеет следующие формы:
COEFS = cwt(S, SCALES, ‘wname’) - возвращает коэффициенты НПВ вещественного или комплексного сигнала S в вещественном положительном SCALES; wname задает имя вейвлета;
- COEFS = cwt(S, SCALES, ‘wname’, ‘plot’) - возвращает коэффициенты и строит их график;
- COEFS = cwt(S, SCALES, ‘wname’, PLOTMODE) возвращает коэффициенты и строит их график с
использованием настроек цвета PLOTMODE.
В итоге НВП получаем вейвлет-спектрограмму (ВС), которая представляет значения коэффициентов вейвлетов в плоскости масштаб (номер коэффициента) - время. В нижней части ВС располагаются коэффициенты с малыми номерами, которые дают детальную картину сигнала, в верхней части - с большими номерами, определяющие огрубленную картину сигнала. При этом их значения характеризуют цвет конкретной (обычно малой) области ВС. Гармоническим сигналам на ВС соответствуют яркие горизонтальные полосы; локальным особенностям (нарушениям гладкости) - вертикальные полосы, выходящие из точки, где имеется особенность. Пики сигналов отображаются сгущением светлых областей ВС; впадины - сгущением темных областей. ?
Пример 6.12. Непрерывный вейвлет-анализ синусоиды.
В командной строке записываются следующие выражения: t=linspace(-6,б,2048); s=sin(t); subplot(211), plot(t,s); title(1 function s(t)") subplot(212), c =cwt(s,1:1:16, "sym4", "abslvl", ); title (1 Wavelet spectr s(t)")|
Запись в третьей строке - "abslvl" - указывает на
характеристику окрашивания: здесь используются абсолютные значения коэффициентов. Результат расчета показан на рис.6.32.
Рис.6.32. График сигнала sin(t) и его вейвлет-спектрограмма
Полученная вейвлет-спектрограмма представляет собой зависимость коэффициентов вейвлет-представления (масштаба) от времени. Для такой простой функции, как sin(t), полезна оценка лишь младших коэффициентов, заданных комплексом 1:1:16 (выводятся коэффициенты от 1 до 16 с шагом 1). Как видно из рис.6.32, в спектрограмме синусоиды нет каких-либо отличительных особенностей, кроме областей перехода сигнала через нуль и экстремальные точки: периодичность синусоиды выражается чередованием темных и светлых областей. Лишь по краям выявляются зоны, усложняющие общий вид спектрограммы, что вызывается ограниченностью во времени области существования сигнала. Напомним, что обычный спектр Фурье такой функции представляет собой одну вертикальную линию с абсциссой, равной частоте синусоиды, и ординатой, определяемой ее амплитудой. ? шума:
s(t) = sin(0,3t) + b(t),
где b(t) - белый шум, равномерно распределенный в интервале [-0,5; 0,5 ]. Выполнить для этого сигнала НВП.
В командной строке запишем:
load noissin; vc=cwt(noissin, 1:48,1db4",1 plot1);
Результат вейвлет-преобразования показан на рис.6.33.
Рис.6.33. Г рафик сигнала (а) и его вейвлет-спектрограмма (б, в)
Аргументы в записанном выражении для cwt определяют использованный вейвлет (Добеши 4 -го порядка) и масштабы анализа. На рис.6.33.б приведены результаты ВА для 48 масштабов, на рис.6.33.в - для 128 масштабов. По существу, возвращаемый аргумент содержит коэффициенты для различных масштабов. В первом случае (рис.6.33.б) результат расчета -
матрица размером 48*1000, в которой каждый ряд соответствует одному фиксированному масштабу, а график дает коэффициенты НВП.
Второй аргумент в командной строке для cwt позволяет управлять числом масштабов (шкал), на которых проводится НВП. При изменении этого параметра необходимо учитывать:
Все масштабы должны быть действительными
положительными числами;
Приращение масштабов должно быть положительным.
Повторим анализ, используя масштабы от 2 до 128, набрав
следующую команду:
с = cwt(noissin,2:2:128,"db4","plot");
Полученный график (рис.6.33.в) дает более ясную картину того, что происходит с сигналом, высвечивая периодичность. ?
ДВП выполняется с помощью команды dwt, которая выполняет одномерную вейвлетную декомпозицию по отношению к конкретному вейвлету ("wname ^ или фильтрам декомпозиции (.Lo_D и Hi_D ).
Синтаксис команды имеет следующий вид:
Dwt(X,’wname") - вычисляет вектор коэффициентов аппроксимации сА и вектор коэффициентов детализации cD, полученных вейвлетным разложением вектора X. Запись в строке "wname" определяет имя вейвлета.
Dwt(X,Lo_D,Hi_D) - вычисляет вейвлетную
декомпозицию, как и выше, но используя указанные здесь фильтры в качестве входов: Lo_D - НЧ-ый фильтр декомпозиции; Hi_D - ВЧ-ый фильтр декомпозиции.
Пример 6.14. Анализируемый сигнал отображает
потребляемую предприятием электроэнергию на протяжении пяти недель с поминутной регистрацией, что формирует выборку, равную 50400 наблюдениям. Часть этого ряда показана на рис. 6.34.а. Выполнить для этого ряда ДВП с использованием вейвлета Добеши.
В командной строке записываются следующие выражения: load leleccum;
s = leleccum(1:3920) ;| l_s = length(s);
Dwt(s, 1 dbl 1);
subplot(1,2,1); plot(Al); title(1 Approximation Al 1) subplot(1,2, 2); plot(Dl); title(1 Detail D1 1)
Результат показан на рис.6.34.
Рис.6.34.
(б) и детализации (в)
Как видно из рис.6.34, коэффициенты аппроксимации достаточно хорошо описывают исходный ряд, а коэффициенты детализации показывают тонкую структуру этого ряда. ?
Режим графического интерфейса
Перейдем от командного окна к режиму графического интерфейса (Graphical Interface ), для чего необходимо набрать в поле командного окна wavemenu.. В результате появляется окно с перечнем основных применений вейвлетов (рис.6.35).
Рис.6.35.
Вначале выполним одномерное НВП, нажав кнопку Continuous Wavelet 1-D. В появившемся окне (рис.6.36) необходимо провести ряд действий для загрузки исследуемого сигнала.
Рассмотрим следующий пример.
Пример 6.15. Выполним НВП для зашумленного синусоидального сигнала, который был использован ранее в примере 6.13.
Загрузим выбранный сигнал с помощью операций, показанных на рис.6.36. Для анализа воспользуемся вейвлетом Добеши 4-го порядка с числом масштабов от 1 до 48. Установив в правой части окна требуемые параметры анализа, получим результат, который показан на рис. 6.37.
Рис.6.36.
Рис.6.37,
На рис.6.37 приведены следующие графики (сверху вниз):
- анализируемый сигнал;
- полученные коэффициенты;
- график коэффициентов, соответствующих масштабу а = 24. Выделив с помощью мыши участок анализируемого сигнала
- (рис.6.38), можно выполнить вейвлет-анализ этой части сигнала (рис.6.39).
Рис.6.38.
Рис.6.39.
Информацию о структуре сигнала можно получить не только в координатах: «Положение (X)» и «Масштаб (Sea)», но и в терминах «Положение (X)» и «Частота (Frq) в герцах. Удерживая правую клавишу мыши на графике коэффициентов, получаем информацию о положении и масштабе (в нижней части окна, рис.6.40). Выполнив такую же операцию, но при отмеченной позиции «Frequency», имеем информацию о положении и частоте (рис.6.41).
Рис.6.40. Вейвлет-анализ в координатах: «Положение» (X = 507) и «Масштаб» (25)
Рис.6.41. Вейвлет-анализ в координатах: «Положение» (X = 600) и «Частота» (0,027)
Часто возникает вопрос о связи масштаба и частоты. Ответ на этот вопрос может быть дан только в широком смысле, поэтому лучше говорить о псевдочастоте
, соответствующей масштабу. Для того чтобы выявить такую связь, необходимо вычислить центральную частоту вейвлета F c ,
равную где а -
масштаб; А - период выборки; F c
- центральная частота вейвлета; F a
- псевдочастота, соответствующая масштабу а
, в Гц.
Сущность такой связи заключается в ассоциации с данным вейвлетом периодического сигнала частоты F c . В системе Matlab имеется функция centfrq, которая может быть использована для вычисления центральной частоты и построения
аппроксимирующего графика вейвлета, основанного на этой частоте. В качестве примера на рис.6.42 приведен график, показывающий периодический сигнал и вейвлет Добеши 7 -го порядка. Как видно из рис.6.42, аппроксимация, основанная на центральной частоте, охватывает основные колебания вейвлета.
Рис.6.42. Вейвлет Добеши 7 -го порядка и периодический сигнал (период: 1,44; частота: 0,7)
Дискретное вейвлет-преобразование
Рассмотрим с помощью графического интерфейса выполнение ДВП на следующем примере.
Пример 6.16. Проведем вейвлет-анализ для сигнала SeriesjG, который имеется во многих статистических пакетах и используется как классический пример сезонных временных рядов.
На рис.6.43 показаны результаты ДВП, выполненного
посредством вейвлета Хаара. Декомпозиция на первом уровне представляет собой сумму коэффициентов аппроксимации и детализации, т.е. 5 = а г + d lt что определяет моду full
decomposition. При ДВП можно более подробно проанализировать структуру сигнала. Так, на рис.6.44 приведены отдельные моды ВА {separate mode): аппроксимации и детализации вместе с
коэффициентами разложения.
Рис.6.43.
Рис.6.44. (а - сигнал и его аппроксимация; б - коэффициенты, сигнал и его
детализация)
Применение вейвлета Добеши 2-го порядка для этого же ряда дает следующие результаты. На рис.6.45 показаны декомпозиция ряда SeriesjG на сумму коэффициентов аппроксимации и детализации в виде
Рис.6.45.
го порядка)
На рис.6.46 приведены отдельные моды BA (separate mode): аппроксимации и детализации вместе с коэффициентами разложения. ?
Рис.6.46. Окно отдельных мод вейвлет-преобразования (а - сигнал и его аппроксимация; б - коэффициенты, сигнал и его
