Опрос

Какой архиватор наиболее эффективный?:

Новички

Виктор Васильев
Юрий Антонов
Сергей Андреевич
Генадий
Avanasy

Преобразование DWT

Информация, которая производится и анализируется в повседневной жизни, является дискретной. Она чаще поступает в виде чисел, а не в форме каких-то непрерывных функций. Поэтому на практике чаще применяются дискретные вейвлетные преобразования (DWT).

Конечно, непрерывные вейвлетные преобразования (CWT, см., например, [Lewalle 95] и [Rao, Bopardikar 98]) также интенсивно изучаются, поскольку это позволяет лучше понять действие DWT.

Преобразование DWT использует свертку, однако из опыта известно, что качество преобразований такого типа сильно зависит от двух вещей: от выбора масштабирующих множителей и временных сдвигов, а также от выбора вейвлета.

На практике преобразование DWT вычисляется с помощью масштабирующих множителей, которые равны отрицательным степеням двойки, и временных сдвигов, которые равны положительным степеням числа 2. На рис. 4.24 показана так называемая двухчленная решетка, которая иллюстрирует такой выбор. Используемые вейвлеты порождают ортонормальные (или биортонормальные) базисы.

Основное направление исследования веивлетов состоит в поисках семейств веивлетов, которые образуют ортогональный базис. Среди этих веивлетов предпочтение отдается веивлетам с компактным носителем, поскольку они позволяют делать преобразования DWT с конечным импульсным откликом (FIR, finite impulse response).

Самый простой способ описания вейвлетных преобразований использует произведение матриц. Этот путь уже был продемонстрирован в § 4.2.1. Преобразование Хаара зависит от двух коэффициентов фильтра со и с\, которые равны 1/\/2 « 0.7071. Наименьшая матрица, которую можно построить в этом случае, равна результату вычисления.

Эта матрица имеет размер 2x2. С ее помощью порождаются два коэффициента преобразования: среднее и разность. (Заметим, что эти среднее и разность не равны в точности полусумме и полуразности, поскольку вместо 2 используется знаменатель у/2. Более точными терминами были бы, соответственно, выражения «грубые детали» и «тонкие детали»). В общем случае, DWT может использовать любое число фильтров, но все они вычисляются с помощью этого метода независимо от вида фильтров.

Сначала мы рассмотрим один из самых популярных вейвлетов, а именно вейвлет Добеши, который принято обозначать D4. Из этого обозначения видно, что он основан на четырех коэффициентах фильтра со, с\, C2 и сз, значения которых приведены в (4.12). Матрица W преобразования равна (ср. с матрицей А\ из (4.1))

Если эту матрицу применить к вектору-столбцу исходных данных (х\,Х2,...,хп), то ее верхняя строка даст взвешенную сумму s\ = = cqX\ +с\Х2 + С2#з + сз#4- Третья строка матрицы определит сумму 52 = со#з + с\х± + С2Х5 + с^х^, и встроки с нечетными номерами зададут аналогичные взвешенные суммы s^. Такие суммы совпадают со свертками исходного вектора Х{ и четырех коэффициентов фильтра. На языке вейвлетов все они называются гладкими коэффициентами, а вместе они именуются сглаживающим фильтром Н.

Аналогично вторая строка матрицы W порождает величину d\ = = С3Х1 — С2Х2 + С1Х3 — С0Х4, а все остальные четные строки матрицы определяют подобные свертки. Каждое число d\ называется детальным коэффициентом, а все вместе они образуют фильтр G. Фильтр G не является сглаживающим. На самом деле, его коэффициенты подобраны так, чтобы фильтр G выдавал на выход маленькие числа, когда данные Х{ коррелированы. Все вместе, Н и G называются квадратурными зеркальными фильтрами (QMF, quadrature mirror filters).

Таким образом вейвлетное преобразование любого изображения представляет собой прохождение исходного образа через фильтр QMF, который состоит из низкочастотного фильтра (Н) и высокочастотного фильтра (G).

Если размер матрицы W равен п х п, то она порождает п/2 гладких коэффициентов S{ и п/2 детальных коэффициентов d. Можно показать, что матрица W будет ортогональной, если четыре порождающие ее коэффициента удовлетворяют соотношениям: Cq + с\ + с\ + Сз = 1, с2со + С3С1 = 0. Еще два уравнения для вычисления коэффициентов фильтра имеют вид: сз — с2 + ci — со = 0 и Осз — 1с2 + 2ci — Зсо = 0. Они получаются из условия равенства нулю первых двух моментов последовательности (сз, —c2,ci, —со). Решением этих четырех уравнений служат следующие величины:

со = (1 + л/3) / (4\/2) « 0.48296, ci = (3 + л/3) / (4д/2) « 0.8365, с2 = (3 - л/3) / (4л/2) « 0.2241, с3 = (1 - л/3) / (4л/2) « 0.1294.

(4.12) Умножение на матрицу W очень просто и наглядно. Однако этот метод не практичен, так как W должна иметь такой же размер, что и исходное изображение, которое обычно велико. Если взглянуть на матрицу W, то видна ее регулярная структура, поэтому нет необходимости строить ее целиком. Достаточно хранить лишь ее верхнюю строку. На самом деле, достаточно иметь массив, состоящий из четырех коэффициентов фильтра.

function wcl=fwtl(dat,coarse,filter)

У, Одномерное прямое дискретное вейвлетное преобразование

У, dat - это вектор-строка размера 2~п,

У, coarse - это самый грубый уровень преобразования

У, (заметим, что coarse должен быть <<п)

У, используется ортогональный квадратурный фильтр,

У, длина которого должна быть <2~(coarse+l)

n=length(dat); j=log2(n); wcl=zeros(l,n); beta^dat;

for i=j-l:-l:coarse

alfa=HiPass(beta,filter);

vcl((2~(i)+l):(2~(i+l)))=alfa;

beta=LoPass(beta,filter) ; end wcl(l:(2"coarse))=beta;

function d=HiPass(dt,filter) У, пропускает высокие частоты

d-iconv(mirror(filter),lshift(dt));

У, iconv - функция свертки из пакета Matlab

n=length(d); d=d(l:2:(n-1));

function d=LoPass(dt,filter) У, пропускает низкие частоты

d=aconv(filter,dt);

У, aconv - это функция свертки из Matlab с обращенным

У, во времени фильтром

n=length(d);

d=d(l:2:(n-D);

function sgn=mirror(filt)

У* возвращает коэффициенты фильтра с обратными знаками

sgn=-((-1). " (1:length(filt))).*f ilt;

Рис. 4.25. Одномерное прямое DWT (Matlab).

На рис. 4.25 показана программа Matlab, которая делает эти вычисления. Функция fwtl(dat,coarse,filter), аргументом dat которой является исходный вектор из 2П элементов, а аргументом filter служит фильтр, вычисляет первые грубые уровни дискретного веивлетного преобразования и записывает их в переменную coarse.

Простой тест для функции fwtl:

n=16; t=(l:n)./n;

dat=sin(2*pi*t)

filt= [0.4830 0.8365 0.2241 -0.1294];

wc=fwtl(dat,1,filt).

На рис. 4.26 приведена программа одномерного дискретного веивлетного преобразования, которое вычисляется с помощью функции iwtl(wc,coarse,filter). Ниже приведен простой тест для ее проверки.

function dat=iwtl(wc,coarse,filter)

'/• Одномерное обратное дискретное вейвлетное преобразование

dat=wc(l:2~coarse);

n=length(wc); j=log2(n);

for i=coarse:j-l

dat=ILoPass(dat,filter)+ . . .

IHiPass(wc((2~(i)+l):(2~(i+l))),filter); end

function f=ILoPass(dt,filter) f=iconv(filter,AltrntZro(dt));

function f=IHiPass(dt,filter)

f=aconv(mirror(filter),rshift(AltrntZro(dt)));

function sgn=mirror(filt)

У. возвращает коэффициенты фильтра с обратными знаками

sgn=-((-1).~(1:length(filt))).*filt;

function f=AltrntZro(dt)

7, возвращает вектор длины 2*п из нулей,

У, размещенных между последовательными значениями

n =length(dt)*2; f =zeros(l,n);

f(l:2:(n-l))=dt;

Рис. 4.26. Одномерное обратное DWT (Matlab).

Простой тест для функции iwtl:

n=16; t=(l:n)./n;

dat=sin(2*pi*t)

filt= [0.4830 0.8365 0.2241 -0.1294];

wc=fwtl(dat,l,filt)

rec=iwt1(wc,1,filt)

Для читателей, которые потрудились разобраться с одномерными функциями fwtl и iwtl (рис. 4.25 и 4.26), мы приводим двумерные аналоги этих функций fwt2 и iwt2 (см. рис. 4.27 и 4.28) для которых также приведена простая тестовая программа.

Кроме семейства фильтров Добеши (между прочим, вейвлет Ха-ара порождает фильтр Добеши степени 2) существует множество других семейств веивлетов, имеющие другие полезные свойства. Вот некоторые известные фильтры: фильтр Белкина, фильтр Койфмана, симметричный фильтр.

Семейство вейвлетов Добеши состоит из ортонормальных функций с компактным носителем, в котором каждая следующая функция имеет большую гладкость, чем предыдущая. В § 4.8 обсуждается вейвлет Добеши D4, а также его «строительный блок». Словосочетание «компактный носитель» означает, что эти функции равны нулю вне некоторого конечного отрезка числовой оси времени.

function wc-fwt2(dat,coarse,filter)

У, Двумерное прямое вейвлетное преобразование

У, dat - это матрица размера (2~п:2~п),

У, «coarse» - самый грубый уровень преобразования

У, (coarse должен быть <<п)

У, длина фильтра <2~(coarse+l)

q=size(dat); n = q(l); j=log2(n);

if q(l) =q(2),

disp('Неквадратное изображение!'), end;

wc = dat; nc = n; for i=j-l:-l:coarse, top = (nc/2+l):nc; bot = l:(nc/2); for ic=l:nc, row = wc(ic,l:nc); wc(ic,bot)=LoPass(row,filter); wc(ic,top)=HiPass(row,filter); end

for ir=l:nc,row = wc(l:nc,ir)'; wc(top,ir)=HiPass(row,filter)'; wc(bot,ir)=LoPass(row,filter)'; end

nc = nc/2; end

function d=HiPass(dt,filter) % пропускает высокие частоты d=iconv(mirror(filter),lshift(dt)); У, iconv - свертка из Mat lab n=length(d); d=d(l:2:(n-1));

function d=LoPass(dt,filter) % пропускает низкие частоты

d=aconv(filter,dt);

У, aconv - это функция свертки из Mat lab с обращенным

У, во времени фильтром

n=length(d);

d-d(l:2:(n-l));

function sgnssmirror(filt)

У» возвращает коэффициенты фильтра с обратными знаками

sgn=-((-1)."(1:length(filt))).*filt;

Рис. 4.27. Двумерное прямое DWT (Matlab).

function dat=iwt2(wc,coarse,filter) У. Двумерное обратное вейвлетное преобразование У, n=length(wc) ; j=log2(n); dat=wc;

nc=2~(coarse+l); for i=coarse:j-1, top=(nc/2+l):nc; bot=l:(nc/2); all=l:nc; for ic=l:nc, dat(all,ic)=ILoPass(dat(bot,ic)',filter)'

+IHiPass(dat(top,ic)',filter)'; end '/, ic for ir=l:nc,

dat(ir,all)=ILoPass(dat(ir,bot),filter)

+IHiPass(dat(ir,top),filter); end '/, ir nc=2*nc; end '/, i

function f-ILoPass(dt,filter) f=iconv(filter,AltrntZro(dt));

function f-IHiPass(dt,filter)

f=aconv(mirror(filter),rshift(AltrntZro(dt)));

function sgn=mirror(filt)

У, возвращает коэффициенты фильтра с обратными знаками

sgn=-((-l).~(l:length(filt))).*filt;

function f=AltrntZro(dt)

У, возвращает вектор длины 2*п из нулей

У, размещенных между последовательными значениями

n =length(dt)*2; f =zeros(l,n);

f(l:2:(n-l))=dt;

Рис. 4.28. Двумерное обратное DWT (Matlab).

Простой тест для функций fwt2 и iwt2:

filename=>house128'; dim=128;

fid=fopen(filename,'r');

if fid=--l dispCfile not found')

else img=fread(fid,[dim,dim])'; fclose(fid);

end

filt= [0.4830 0.8365 0.2241 -0.1294];

fwim=fwt2(img,4,filt) ;

figure(l), imagesc(fwim), axis off, axis square

rec=iwt2(fwim,4,filt);

figure(2), imagesc(rec), axis off, axis square

Вейвлет Добеши D4 строится по четырем коэффициентам, приведенным в (4.12). Аналогично, вейвлет D6 имеет шесть коэффициентов. Их можно найти, решив систему из шести уравнений, три из которых отражают свойство ортонормальности, а другие три получаются из условия равенства нулю первых трех моментов. Результат приведен в (4.13).

со = (l + >/1б + \Д + 2у/Тб) I (16v/2) « .3326,

Cl = 15 + л/ГО + Зл/5 + 2VTo) / (16л/2) ~ -8068,

с2 = (10 - 2>/10 + 2\/5 + 2ч/1о) / (16л/2) « .4598,

) ________/ ' (4 13)

с3 = f 10 - 2VT5 - 2^5 + 2х/Ш) / (16л/2) « -.1350,

с4 = (б + \/1б - 3\/5 + 2л/10) / (16л/2) « -.0854,

с5 = (l + VT0 - \/5 + 2л/Ш) / (16л/2) « .0352.

В каждой последовательности этого семейства число коэффициентов на два больше, чем в предыдущей, причем они являются более гладкими. Происхождение этих функций обсуждается в [Daubechies 88], [DeVore и др. 92] и [Vetterli, Kovacevic 95].