Материал: Применение вейвлет-преобразований

Пример просмотра процесса итерации вейвлета 'sym2'

iter=10; wav='db2';i=1: iter

[phi,psi,x]=wavefun(wav,i); plot(x,psi); hold on; hold off

Рис. 3.1

Примечание: Текст программы можно копировать и подставлять непосредствен но в командную строку MATLAB после знака ввода '>>'.

Средняя частота вейвлета определяет положение максимума пика его Фурье-образа на оси частот (в герцах). Формат функции:

Q=centfrq('wname'[,ITER]),

где ITER - количество итераций, если вейвлет 'wname' не задан аналитической функцией и вычисляется методом итераций.

Пример определения средней частоты и построения сравнительных графиков вейвлета db2 и гармоники с частотой, равной средней частоте вейвлета, приведен на рис. 3.2.

>> f=centfrq('db2',16,'plot'); grid=0.6667

Рис. 3.2.

3.2 Различные примеры функций вейвлет-анализа в MATLAB.

Функция intwave - первообразная вейвлет-функции .

Вычисляет первообразную функции  как интеграл функции  вейвлета ‘wname’ от  до переменной точки .

[integ,X]=intwave(‘wname’,N),

[integ,X]=intwave(‘wname’).

Значения параметров:

 - количество итераций(по умолчанию );

- массив значений переменной . Он состоит из точек, равномерно расположенных на носителе вейвлета с шагом . Если носитель имеет длину , то длина массива  равна  точка сетки , т.е. концевые точки носителя включаются; integ - значения первообразной функции  в точках . Функция  аппроксимируется на  точках сетки . Выходной массив integ является действительным или комплексным в зависимости от типа вейвлета.

Для биортогональных вейвлетов функция [intdec,x,intrec]=

=intwave(‘wname’,N) вычисляет интегралы intdec и intrec вейвлет-функций разложения  и восстановления . Выражение [integ,X]=intwave(‘wname’) эквивалентно [integ,X]= intwave(‘wname’,8).

Пример 1. Для нахождения первообразной вейвлет-функции Добеши  запускаем программу MATLAB. После чего перед нами появляется главное меню программы MATLAB, показанный на рис. 3.3.

Рис. 3.3. Главное меню программы MATLAB

Далее в командной строке вводим следующий алгоритм для нахождения первообразной данной вейвлет-функции:

[phi,psi,xval]=wavefun('db4',7);(211);plot(xval,psi);title('Wavelet');

[integ,xval]=intwave('db4',7);(212);plot(xval,integ);(['Wavelet integrals over [-Inf x] for each value of xval']);

После того как ввод алгоритма закончен, нажимая на клавишу Enter мы вызываем графическое решение примера. Появляться окно с двумя графиками, показанные на рис. 3.4. Первый график - это график самой вейвлет-функции Добеши, а второй - график первообразной этой функции.

Рис. 3.4. График вейвлет-функции Добеши и его первообразной.

Одноуровневое дискретное одномерное вейвлет-преобразование dwt.

В одноуровневом дискретном одномерном вейвлет-преобразовании используются следующие функции:

·    dwt - одноуровневое дискретное одномерное вейвлет-преобразование;

·        idwt - обратное одноуровневое дискретное одномерное вейвлет-преобразование;

·        dwtmode - мода расширения дискретного вейвлет-преобразования.

Дискретное вейвлет-преобразование. функция вычисляет аппроксимирующие коэффициенты cA и детализирующие коэффициенты cD первого уровня вейвлет-преобразования сигнала, заданного вектором . Коэффициенты получаются сверткой сигнала  с фильтром разложения нижних частот Lo_D для аппроксимации с фильтром разложения нижних частот Hi_D для детализации, а затем сопровождаются двоичной децимацией.

Команда dwt употребляются в следующих формах:

[cA,cD]=dwt(s,’wname’),[cA,cD]=dwt(s,’wname’,’mode’,MODE),

Функция [cA,cD]=dwt(s,’wname’) вычисляет аппроксимирующие коэффициенты cA и детализирующие коэффициенты cD вектора s. Строка wname - имя вейвлета, для которого возможны быстрые алгоритмы. Это - haar, db, sym, coif, bior, rbio, dmey. Функия [cA,cD]=dwt(s,Lo_D,Hi_D) находит вейвлет разложение, как выше, используя вместо имени вейвлета его низкочастотный фильтр разложения Lo_D и высокочастотный фильтр разложения Hi_D.

Параметр MODE определяет метод расширения данных для вейвлет-разложения. Например, [cA,cD]=dwt(s,’db2’,’mode’’sym’).

Пример 2. Зададим сигнал . Найдем его разложение, используя вейвлет Хаара, затем - Добеши .

n=1:20;s=sin(0.5*n);

[cal,cd1]=dwt(s,'haar');(311);plot(s);title('Original signal');(323);plot(cal);title('Approx. coef. for haar');(324);plot(cd1);title('Detail coef. for haar');

[ca2,cd2]=dwt(s,'db4');(325);plot(ca2);title('Approx. coef. for db2');(326);plot(cd2);title('Detail coef. for db2');

Ниже на рис. 3.5 приведено графическое решение.

Рис. 3.5 Разложение данного сигнала используя вейвлеты Хаара и Добеши

Непрерывное вейвлет-преобразование cwt.

Команда cwt вычисляет вейвлет-коэффициенты одномерного сигнала, используя непрерывное вейвлет-преобразование.

Пусть  - сигнал и  - вейвлет. При непрерывном преобразовании вейвлет коэффициенты сигнала , соответствующие масштабному коэффициенту  и положению , определяются формулой:

.

Если сигнал  является дискретным, заданным выборкой значений , то в качестве  можно взять кусочно-постоянную функцию, определенную следующим образом:  , при . Параметр  пробегает значения с шагом , соответствующие длине сигнала , т.е. . Полагаем, что значения масштабного параметра  также дискретны и составляют вектор, который далее обозначается SCALES.

Основное отличие команды cwt от дискретного многоуровневого вейвлет-преобразования wavedec заключается в том, что в cwt можно брать любые значения масштабного коэффициента , тогда как в wavedec используются степени двойки.

Для команды cwt предусмотрена возможность графического изображения вейвлет-коэффициентов сигнала , называемое вейвлет-спектрограммой сигнала.

Синтаксис cwt:

=cwt(s,SCALES,’wname’).=cwt(s,SCALES,’wname’,’plot’),=cwt(s,SCAKES,’wname’,PLOIMODE),=cwt(s,SCALES,’wname’,PLOTMODE,XLIM).

Функция COEFS=cwt(s,SCALES,’wname’) вычисляет вейвлет-коэффициенты вектора , когда  принимает положительные значения из вектора SCALES, используя вейвлет ’wname’. Сигнал  вещественный, вейвлет может быть вещественным или комплексным.

Для каждого значения масштабного коэффициента  в пределах вектора SCALES вейвлет-коэффициенты  вычисляются для и и сохраняются в виде вектора-строки COEFS(i,:), соответствующей -му значению параметра a=SCALES(i).

Результат COEFS является  на  матрицей, где  является длиной вектора SCALES, . Результат COEFS является вещественной или комплексной матрицей в зависимости от типа вейвлета.

Функция COEFS=cwt(s,SCALES,’wname’,’plot’) вычисляет и создает вейвлет-спектрограмму сигнала .

Таблица 3.1

По умолчанию принят вариант окраски absglb, т.е.

COEFS=cwt(…,’plot’) эквивалентно COEFS=cwt(…,’absglb’).

Функция COEFS=cwt(s,SCALES,’wname’,PLOTMODE,XLIM) вычисляет и создает вейвлет-спектрограмму сигнала  с указанием пределов изображения XLIM=[x1 x2] по горизонтали, .

Пример 3. Варианты написания команды.

=linspace(-1,1,512);s=1-abs(t);

c=cwt(s,1:32,'cgau4');=cwt(s,[64 32 16:-2:2 1.5],'morl');=cwt(s,1:64,'sym4','abslvl',[100 400]);

Результаты показаны на рис. 3.6.

Рис. 3.6 Вейвлет коэффициенты одномерного сигнала

Функции в MATLAB для удаления шума.

Функция wden автоматическое удаление шума(1-D).

Функция wden автоматически удаляет шум одномерного сигнала , используя вейвлеты. Применяются в виде:

[XD,CXD,LXD]=wden(X,TPTR,SOHR,SCAL,N,’wname’),

[XD,CXD,LXD]=wden(C,TPTR,SOHR,SCAL,N,’wname’).

Первая версия функции вычисляет очищенную от шума версию XD выходного сигнала , используя пороговую обработку вейвлет-коэффициентов. Дополнительно находит вейвлет-коэффицинеты [CXD,LXD] разложения очищенного сигнала XD.

Исходные параметры:

·        - сигнал;

·        TPTR - строка, определяющая выбор порога:

- ‘rigrsure’ - адаптивный выбор порога, используя принцип Штейна и несмещенной оценки риска (SURE);

?  ‘heursure’ - эвристический вариант первого выбора;

?      ‘qtwolog’ - универсальный порог sqrt(2*log(length(X)));

?      ‘minimaxi’ - мини-максный порог;

• SOHR=’s’ или ‘h’ - выбор мягкого или жесткого;

• SCAL - строка, которая определяет мультипликативное пороговое перемасштабирование. Если шум вне пределов  или не белый, то порог должен быть перемасштабирован с использованием оценки уровня шума. Возможны следующие варианты:

‘one’ - без перемасштабирования, когда используется базовый метод;

‘sln’ - для перемасштабирования порога с использованием оценки уровня шума на базе коэффициентов первого уровня;

‘mln’ - для перемасштабирования порога с использованием оценки уровня шума, зависящей от уровня разложения;

• N - уровень вейвлет-разложения и ’wname’ - имя ортогонального вейвлета.

Во втором случае функция wden производит тоже самое, но использует уже полученное вейвлет-разложение [C,L] исходного сигнала до уровня N и ортогональный вейвлет ’wname’.

Пример 4. Создаем зашумленный сигнал функцией wnoise. Для удаления шума используем мягкий эвристический метод SURE. Разложение  берем до уровня 5 по вейвлету sym8.

=3;init=2055615866;

[xref,x]=wnoise(3,11,snr,init);

xd=wden(x,'heursure','s','one',5,'sym8');(311),plot(xref),axis([1 2048 -10 10]);title('original signal');(312),plot(x),axis([1 2048 -10 10]);title('Noisy signal');(313),plot(xd),axis([1 2048 -10 10]);title('De-noised signal - heuristic SURE')

Результат показан на рис. 3.7.

Рис. 3.7. График данного сигнала, зашумленного и обесшумленного соответственно.

Функция wdencmp - удаление шума и сжатие при помощи вейвлетов(1-D,2-D).

[XC,CXC,LXC,PRO,PRL2]=(‘gbl’,X,’wname’,N,THR,SORH,KEEPAPP),

[XC,CXC,LXC,PRO,PRL2]=wdencmp(‘lvd’,X,’wname’,N,THR,SORH),

[XC,CXC,LXC,PRO,PRL2]=wdencmp(‘lvd’,C,L,’wname’,N,THR,SORH).

Описание. Случай параметра ‘gbl’. Функция производит:

·   XC - очищенную от шума версию входного сигнала X (1-D или 2-D);

·        [CXC,LXC] - структуру вейвлет-разложения сигнала XC;

·        PRO - оценку числа обнуленных коэффициентов сигнала X в процентах;

·        PRL2 - отношение -норм сжатого и первоначального сигнала в процентах. PRL2=100(norm(CXC)/norm(С))2, где [C,L] - вейвлет-разложение X. Если X - одномерный сигнал и ‘wname’ - ортогональный вейвлет, то PRL2=100.

Функция использует:

·   X - сигнал, или изображение;

·        ‘wname’ - имя вейвлета и N - уровень вейвлет разложения;

·        THR - значение порога;

·        SORH=’s’ или ‘h’ - выбор мягкого или жесткого порогового метода;

·        KEEPAPP=1, то коэффициенты аппроксимации не подвергаются пороговой обработке, иначе это возможно.

Функция wdencmp(‘gbl’,X,’wname’,N,THR,SORH,KEEPAPP) выполняет то же самое, используя вейвлет-разложение [C,L] исходного сигнала до уровня N.

В одномерном случае и опции ‘lvd’. Функция выполняет то же самое, но используя зависимые от уровня пороги, указанные в векторе THR (THR должен иметь длину N). Кроме того, коэффициенты аппроксимации сохраняются.

Для двумерного случая и опции ‘lvd’. В этом случае THR должен быть матрицей , содержащей зависимые от уровня пороги в трех направлениях: по горизонтали, по диагонали и по вертикали.

Отметим, что сравнительно с wden (автоматическое удаление шума), wdencmp дает больше возможностей и можно осуществить свою стратегию удаления шума. Пример 5. Загружаем изображение девушки. Делаем вейвлет-разложение  при . Используем wdencmp для компрессии изображения.