Этот пример показывает, как двойной древовидный комплексный дискретный вейвлет преобразовывает (DT-CWT), обеспечивает преимущества перед критически выбранным DWT для сигнала, изображения и обработки объема. Двойной древовидный DWT реализован, когда два разделяют двухканальные наборы фильтров. Получать преимущества описало в этом примере, вы не можете произвольно выбрать масштабирование и фильтры вейвлета, используемые в этих двух деревьях. lowpass (масштабирование) и highpass (вейвлет) фильтры одного дерева
, должен сгенерировать масштабирующуюся функцию и вейвлет, которые являются аппроксимированными Гильбертовыми преобразованиями масштабирующейся функции и вейвлета, сгенерированного lowpass и highpass фильтрами другого дерева
. Поэтому функции масштабирования с комплексным знаком и вейвлеты, сформированные из этих двух деревьев, приблизительно аналитичны.
В результате двойной древовидный DWT показывает меньше отклонения сдвига и больше направленной селективности, чем критически выбранный DWT только с
фактором сокращения для 
- размерные данные. Сокращение в двойном древовидном DWT является значительно меньше, чем сокращение в неподкошенном (стационарном) DWT.
Этот пример иллюстрирует, что аппроксимированная инвариантность сдвига двойного древовидного DWT, выборочная ориентация вейвлетов анализа двойного дерева в 2D и 3D, и использование двойного древовидного комплексного вейвлета преобразовывают в шумоподавление объема и изображение.
DWT страдает от отклонения сдвига, означая, что маленькие сдвиги во входном сигнале или изображении могут вызвать существенные изменения в распределении энергии сигнала/изображения через шкалы в коэффициентах DWT. DT-CWT является приблизительно инвариантом сдвига.
Чтобы продемонстрировать это на тестовом сигнале, создайте два переключенных импульса дискретного времени 128 выборок в длине. Один сигнал имеет модульный импульс в демонстрационных 60, в то время как другой сигнал имеет модульный импульс в демонстрационных 64. Оба сигнала ясно имеют модульную энергию
(норма).
kronDelta1 = zeros(128,1); kronDelta1(60) = 1; kronDelta2 = zeros(128,1); kronDelta2(64) = 1;
Получите DWT и двойной древовидный DWT двух сигналов вниз к уровню 3 с вейвлетом и масштабирующимися фильтрами длины 14. Извлеките коэффициенты детали уровня 3 для сравнения.
J = 3; dwt1 = dddtree('dwt',kronDelta1,J,'sym7'); dwt2 = dddtree('dwt',kronDelta2,J,'sym7'); dwt1Cfs = dwt1.cfs{J}; dwt2Cfs = dwt2.cfs{J}; dt1 = dddtree('cplxdt',kronDelta1,J,'dtf3'); dt2 = dddtree('cplxdt',kronDelta2,J,'dtf3'); dt1Cfs = dt1.cfs{J}(:,:,1)+1i*dt1.cfs{J}(:,:,2); dt2Cfs = dt2.cfs{J}(:,:,1)+1i*dt2.cfs{J}(:,:,2);
Постройте абсолютные значения DWT и коэффициентов DT-CWT для двух сигналов на уровне 3 и вычислите энергию (нормы в квадрате) коэффициентов.
figure; subplot(1,2,1) stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1]); title(['DWT Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2,3)],... 'fontsize',10); subplot(1,2,2) stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1]) title(['DWT Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)],... 'fontsize',10); figure; subplot(1,2,1) stem(abs(dt1Cfs),'markerfacecolor',[0 0 1]); title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)],... 'fontsize',10); subplot(1,2,2) stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1]) title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)],... 'fontsize',10);
Обратите внимание, что четыре выборки переключаются на нижний регистр, сигнал вызвал изменение на почти 6,5% в энергии уровня 3 коэффициенты вейвлета DWT. Однако двойные древовидные коэффициенты вейвлета показывают только изменение на 0,3% в энергии.
Чтобы продемонстрировать утилиту аппроксимированной инвариантности сдвига в действительных данных, мы анализируем электрокардиограмму (ECG) сигнал. Интервал выборки для сигнала ECG является 1/180 секундами. Данные взяты из Percival & Walden (2000), p.125 (данные, первоначально обеспеченные Уильямом Константином и На Reinhall, Вашингтонский университет). Для удобства мы берем данные, чтобы запуститься в t=0.
load wecg dt = 1/180; t = 0:dt:(length(wecg)*dt)-dt; figure; plot(t,wecg) xlabel('Seconds'); ylabel('Millivolts');
Большой положительный peaks на расстоянии приблизительно в 0,7 секунды является волнами R сердечного ритма. Во-первых, анализируйте сигнал с помощью критически выбранного DWT и постройте исходный сигнал наряду с коэффициентами вейвлета уровня 2 и уровня 3. Коэффициенты уровня 2 и уровня 3 были выбраны, потому что волны R изолируются наиболее заметно в тех шкалах для данной выборки уровня.
J = 6; dtDWT1 = dddtree('dwt',wecg,J,'farras'); details = zeros(2048,3); details(2:4:end,2) = dtDWT1.cfs{2}; details(4:8:end,3) = dtDWT1.cfs{3}; subplot(311) stem(t,details(:,2),'Marker','none','ShowBaseline','off') title('Level 2'); ylabel('mV'); subplot(312) stem(t,details(:,3),'Marker','none','ShowBaseline','off') title('Level 3'); ylabel('mV'); subplot(313) plot(t,wecg); title('Original Signal'); xlabel('Seconds'); ylabel('mV');
Повторите вышеупомянутый анализ для двойного древовидного преобразования. В этом случае только постройте действительную часть двойных древовидных коэффициентов на уровнях 2 и 3.
dtcplx1 = dddtree('cplxdt',wecg,J,'dtf3'); details = zeros(2048,3); details(2:4:end,2) = dtcplx1.cfs{2}(:,1,1)+1i*dtcplx1.cfs{2}(:,1,2); details(4:8:end,3) = dtcplx1.cfs{3}(:,1,1)+1i*dtcplx1.cfs{3}(:,1,2); subplot(311) stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off') title('Level 2'); ylabel('mV'); subplot(312) stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off') title('Level 3'); ylabel('mV'); subplot(313) plot(t,wecg); title('Original Signal'); xlabel('Seconds'); ylabel('mV');
И критически выбранный и двойной древовидный вейвлет преобразовывает, локализуют важную функцию формы волны ECG к подобным шкалам.
Важное приложение вейвлетов в 1D сигналах должно получить дисперсионный анализ шкалой. Это выдерживает обосновать, что этот дисперсионный анализ не должен быть чувствителен к циклическим сдвигам во входном сигнале. К сожалению, дело обстоит не так с критически выбранным DWT. Чтобы продемонстрировать это, мы циркулярный сдвиг сигнал ECG 4 выборками, анализирует непереключенные и переключенные сигналы с критически выбранным DWT и вычисляет распределение энергии через шкалы.
wecgShift = circshift(wecg,4); dtDWT2 = dddtree('dwt',wecgShift,J,'farras'); sigenrgy = norm(wecg,2)^2; enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT1.cfs,'uni',0)); enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT2.cfs,'uni',0)); levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'}; enr1 = enr1(:); enr2 = enr2(:); table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})
ans =
7x3 table
Level enr1 enr2
_____ ______ ______
'D1' 4.1994 4.1994
'D2' 8.425 8.425
'D3' 13.381 10.077
'D4' 7.0612 10.031
'D5' 5.4606 5.4436
'D6' 3.1273 3.4584
'A6' 58.345 58.366
Обратите внимание на то, что коэффициенты вейвлета на уровнях 3 и 4 показывают приблизительно 3%-е изменения в энергии между исходным и переключенным сигналом. Затем, мы повторяем, что этот анализ с помощью комплексного двойного древовидного вейвлета преобразовывает.
dtcplx2 = dddtree('cplxdt',wecgShift,J,'dtf3'); cfs1 = cellfun(@squeeze,dtcplx1.cfs,'uni',0); cfs2 = cellfun(@squeeze,dtcplx2.cfs,'uni',0); cfs1 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs1,'uni',0); cfs2 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs2,'uni',0); dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0)); dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0)); dtenr1 = dtenr1(:); dtenr2 = dtenr2(:); table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})
ans =
7x3 table
Level dtenr1 dtenr2
_____ ______ ______
'D1' 4.936 4.936
'D2' 6.6691 6.6691
'D3' 12.682 12.611
'D4' 8.3891 8.4808
'D5' 5.8868 5.8791
'D6' 3.053 3.0415
'A6' 58.384 58.382
Двойное древовидное преобразование производит сопоставимый дисперсионный анализ шкалой для исходного сигнала и его циркулярной переключенной версии.
Стандартная реализация DWT в 2D использовании отделимая фильтрация столбцов и строки изображения. LH, HL и вейвлеты HH для Добечиса меньше всего - асимметричный вейвлет фазы с 4 исчезающими моментами (sym4) показывают в следующей фигуре.
figure; J = 5; L = 3*2^(J+1); N = L/2^J; Y = zeros(L,3*L); dt = dddtree2('dwt',Y,J,'sym4'); dt.cfs{J}(N/3,N/2,1) = 1; dt.cfs{J}(N/2,N/2+N,2) = 1; dt.cfs{J}(N/2,N/2+2*N,3) = 1; dwtImage = idddtree2(dt); imagesc(dwtImage); axis xy; axis off; title({'Critically Sampled DWT';'2D separable wavelets (sym4) -- LH, HL, HH'});
Обратите внимание на то, что LH и вейвлеты HL имеют ясные горизонтальные и вертикальные ориентации соответственно. Однако вейвлет HH на ультраправых смесях и +45 и-45 направлений степени, производя артефакт шахматной доски. Это смешивание ориентаций происходит из-за использования отделимых фильтров с действительным знаком. HH отделимый фильтр с действительным знаком имеет полосы пропускания во всех четырех высокочастотных углах 2D плоскости частоты.
Двойной древовидный DWT достигает направленной селективности при помощи вейвлетов, которые приблизительно аналитичны, означая, что у них есть поддержка только на одной половине оси частоты. В двойном древовидном DWT существует шесть поддиапазонов и для действительных и для мнимых частей. Эти шесть действительных частей сформированы путем добавления выходных параметров фильтрации столбца, сопровождаемой фильтрацией строки входного изображения в этих двух деревьях. Эти шесть мнимых частей сформированы путем вычитания выходных параметров фильтрации столбца, сопровождаемой фильтрацией строки.
Фильтры применились к столбцам, и строки могут быть от той же пары фильтра,
или
, или от различных пар фильтра
. Следующий код показывает, что ориентация этих 12 вейвлетов, соответствующих действительным и мнимым частям комплекса, ориентировала двойной древовидный DWT.
J = 4; L = 3*2^(J+1); N = L/2^J; Y = zeros(2*L,6*L); wt = dddtree2('cplxdt',Y,J,'dtf3'); wt.cfs{J}(N/2,N/2+0*N,2,2,1) = 1; wt.cfs{J}(N/2,N/2+1*N,3,1,1) = 1; wt.cfs{J}(N/2,N/2+2*N,1,2,1) = 1; wt.cfs{J}(N/2,N/2+3*N,1,1,1) = 1; wt.cfs{J}(N/2,N/2+4*N,3,2,1) = 1; wt.cfs{J}(N/2,N/2+5*N,2,1,1) = 1; wt.cfs{J}(N/2+N,N/2+0*N,2,2,2) = 1; wt.cfs{J}(N/2+N,N/2+1*N,3,1,2) = 1; wt.cfs{J}(N/2+N,N/2+2*N,1,2,2) = 1; wt.cfs{J}(N/2+N,N/2+3*N,1,1,2) = 1; wt.cfs{J}(N/2+N,N/2+4*N,3,2,2) = 1; wt.cfs{J}(N/2+N,N/2+5*N,2,1,2) = 1; waveIm = idddtree2(wt); imagesc(waveIm); axis off; title('Complex Dual-Tree 2D Wavelets');
Верхняя строка предыдущей фигуры показывает, что шесть направленных вейвлетов действительного ориентированного двойного древовидного вейвлета преобразовывают. Вторая строка показывает мнимые части. Вместе действительные и мнимые части формируются, комплекс ориентировался, двойной древовидный вейвлет преобразовывают. Действительные и мнимые части ориентированы в том же направлении. Можно использовать dddtree2 с опцией 'realdt', чтобы получить действительный ориентированный двойной древовидный DWT, который использует только действительные части. Используя действительный ориентированный двойной древовидный вейвлет преобразовывают, можно достигнуть направленной селективности, но вы не получаете полную выгоду от использования аналитических вейвлетов, таких как аппроксимированная инвариантность сдвига.
Аппроксимированная аналитичность и выборочная ориентация комплексных двойных древовидных вейвлетов обеспечивают наилучшее решение по стандартному 2D DWT в представлении ребер в изображениях. Чтобы проиллюстрировать это, мы анализируем тестовые изображения с ребрами, состоящими из строки и особенностей кривой в нескольких направлениях с помощью критически выбранного 2D DWT, и 2D комплекс ориентировался, двойное дерево преобразовывают. Во-первых, анализируйте изображение восьмиугольника, который состоит из особенностей строки.
load woctagon; figure; imagesc(woctagon); colormap gray; title('Original Image'); axis equal; axis off;
Анализируйте изображение вниз к уровню 4 и восстановите приближение изображений на основе коэффициентов детали уровня 4.
dtcplx = dddtree2('cplxdt',woctagon,4,'dtf3'); dtDWT = dddtree2('dwt',woctagon,4,'farras'); dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1})); dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2})); dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3})); dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5})); dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1})); dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2})); dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3})); dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5})); dtImage = idddtree2(dtcplx); dwtImage = idddtree2(dtDWT); subplot(121) imagesc(dtImage); axis equal; axis off; colormap gray; title('Complex Oriented Dual-Tree'); subplot(122) imagesc(dwtImage); axis equal; axis off; colormap gray; title('DWT')
Затем, анализируйте восьмиугольник с гиперболическими ребрами. Ребра в гиперболическом восьмиугольнике являются особенностями кривой.
load woctagonHyperbolic; figure; imagesc(woctagonHyperbolic); colormap gray; title('Octagon with Hyperbolic Edges'); axis equal; axis off;
Снова, анализируйте изображение вниз к уровню 4 и восстановите приближение изображений на основе коэффициентов детали уровня 4 и для стандартного 2D DWT и для комплекса, ориентированного на двойной древовидный DWT.
dtcplx = dddtree2('cplxdt',woctagonHyperbolic,4,'dtf3'); dtDWT = dddtree2('dwt',woctagonHyperbolic,4,'farras'); dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1})); dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2})); dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3})); dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5})); dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1})); dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2})); dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3})); dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5})); dtImage = idddtree2(dtcplx); dwtImage = idddtree2(dtDWT); subplot(121) imagesc(dtImage); axis equal; axis off; colormap gray; title('DT-CWT'); subplot(122) imagesc(dwtImage); axis equal; axis off; colormap gray; title('DWT')
Обратите внимание на то, что звонящие артефакты, очевидные в 2D критически, выбрали DWT, отсутствуют в 2D DT-CWT обоих изображений. DT-CWT более искренне воспроизводит особенности кривой и строка.
Из-за способности изолировать отличные ориентации в отдельных поддиапазонах, двойной древовидный DWT часто может превзойти стандартный отделимый DWT по характеристикам в приложениях как шумоподавление изображений. Чтобы продемонстрировать это, используйте функцию помощника helperCompare2DDenoising. Функция помощника загружает изображение и добавляет нулевой средний белый Гауссов шум с
. Для предоставленной пользователями области значений порогов функция сравнивает шумоподавление с помощью мягкой пороговой обработки для критически выбранного DWT, действительного ориентированного двойного древовидного DWT, и комплекс ориентировал двойной древовидный DWT. Для каждого порогового значения отображены среднеквадратичная (RMS) ошибка и пиковое отношение сигнал-шум (PSNR).
numex = 3;
helperCompare2DDenoising(numex,0:2:100,'PlotMetrics');
И ориентированное действительное и комплекс ориентировались, двойные древовидные DWTs превосходят стандартный DWT по характеристикам по ошибке RMS и PSNR.
Затем, получите изображения denoised для порогового значения 25, который равен стандартному отклонению аддитивного шума.
numex = 3;
helperCompare2DDenoising(numex,25,'PlotImage');
С пороговым значением, равным стандартному отклонению аддитивного шума, ориентировался комплекс, двойное дерево преобразовывают, обеспечивает PSNR на почти 4 дБ выше, чем стандартный 2D DWT.
Звонящие артефакты, наблюдаемые с отделимым DWT в двух измерениях, усилены при расширении анализа вейвлета к более высоким размерностям. DT-CWT позволяет вам поддержать направленную селективность в 3D с минимальным сокращением. В 3D в двойном древовидном преобразовании существует 28 поддиапазонов вейвлета.
Чтобы продемонстрировать направленную селективность 3D двойного древовидного вейвлета преобразовывают, визуализируют пример 3D изоповерхности и 3D двойного дерева и отделимых вейвлетов DWT. Во-первых, визуализируйте действительные и мнимые части отдельно двух двойных древовидных поддиапазонов.
helperVisualize3D('Dual-Tree',28,'separate'); helperVisualize3D('Dual-Tree',25,'separate');
Красный фрагмент графика изоповерхности указывает на положительное отклонение вейвлета от нуля, в то время как синий обозначает отрицательное отклонение. Можно ясно видеть направленную селективность на пробеле действительных и мнимых частей двойных древовидных вейвлетов. Теперь визуализируйте один из двойных древовидных поддиапазонов с действительными и мнимыми графиками, построенными вместе как одна изоповерхность.
helperVisualize3D('Dual-Tree',25,'real-imaginary');
Предыдущий график демонстрирует, что действительные и мнимые части являются переключенными версиями друг друга на пробеле. Это отражает то, что мнимая часть комплексного вейвлета является аппроксимированным Гильбертовым преобразованием действительной части. Затем, визуализируйте изоповерхность действительного ортогонального вейвлета в 3D для сравнения.
helperVisualize3D('DWT',7);
Смешивание ориентаций, наблюдаемых в 2D DWT, является еще более явным в 3D. Так же, как в 2D случае, смешивание ориентаций в 3D приводит к значительному вызову или блокированию артефактов. Чтобы продемонстрировать это, исследуйте 3D DWT и детали вейвлета DT-CWT сферического объема. Сфера 64 64 64.
load sphr; [A,D] = dualtree3(sphr,2,'excludeL1'); A = zeros(size(A)); sphrDTCWT = idualtree3(A,D); montage(reshape(sphrDTCWT,[64 64 1 64]),'DisplayRange',[]) title('DT-CWT Level 2 Details');
Сравните предыдущий график против деталей второго уровня на основе отделимого DWT.
sphrDEC = wavedec3(sphr,2,'sym4','mode','per'); sphrDEC.dec{1} = zeros(size(sphrDEC.dec{1})); for kk = 2:8 sphrDEC.dec{kk} = zeros(size(sphrDEC.dec{kk})); end sphrrecDWT = waverec3(sphrDEC); figure; montage(reshape(sphrrecDWT,[64 64 1 64]),'DisplayRange',[]) title('DWT Level 2 Details');
Увеличьте масштаб изображений и в DT-CWT и в монтаже DWT, и вы будете видеть, как видный блокирующиеся артефакты в деталях DWT сравниваются с DT-CWT.
Подобно 2D случаю направленная селективность 3D DT-CWT часто приводит к улучшениям шумоподавления объема.
Чтобы продемонстрировать это, рассмотрите набор данных MRI, состоящий из 16 срезов. Гауссов шум со стандартным отклонением 10 был добавлен к исходному набору данных. Отобразите шумный набор данных.
load MRI3D; montage(reshape(noisyMRI,[128 128 1 16]),'DisplayRange',[]);
Обратите внимание на то, что исходный ОСШ до шумоподавления составляет приблизительно 11 дБ.
20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))
ans = 11.2997
Denoise набор данных MRI вниз к уровню 4 с помощью и DT-CWT и DWT. Подобные длины фильтра вейвлета используются в обоих случаях. Постройте получившийся ОСШ как функцию порога. Отобразите результаты denoised и для DT-CWT и для DWT, полученного в лучшем ОСШ.
[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI); figure; montage(reshape(imrecDTCWT,[128 128 1 16]),'DisplayRange',[]); title('DT-CWT Denoised Volume'); figure; montage(reshape(imrecDWT,[128 128 1 16]),'DisplayRange',[]); title('DWT Denoised Volume');
Мы показали, что двойной древовидный DWT обладает желательными свойствами близкой инвариантности сдвига и направленной селективности, не достижимой с критически выбранным DWT. Мы продемонстрировали, как эти свойства могут привести к улучшенной производительности в анализе сигнала, представлении ребер в изображениях и объемах, и шумоподавлении объема и изображении.
Chen, H. & Kingsbury, N.G. "Efficient registration of nonrigid 3-D bodies", IEEE Transactions on Image Processing, Vol. 21, No. 1, Jan 2012, pp. 262-272.
Kingsbury, N.G. "Complex Wavelets for Shift Invariant Analysis and Filtering of Signals". Journal of Applied and Computational Harmonic Analysis. Vol 10, Number 3, May 2001, pp. 234-253.
Percival, D.B. and A.T. Walden. "Wavelet Methods for Time Series Analysis", Cambridge University Press, 2000.
Selesnick, I., Baraniuk, R.G., and N.G. Kingsbury. "The Dual-Tree Complex Wavelet Transform." IEEE Signal Processing Magazine. Vol. 22, Number 6, November, 2005, pp. 123-151.
Selesnick, I. "The Double Density DWT". Wavelets in Signal and Image Analysis: From Theory to Practice (A.A Petrosian, F.G. Meyer, eds.), Norwell, MA: Kluwer Academic Publishers, 2001, pp. 39-66.
Selesnick, I. "The Double-Density Dual-Tree Wavelet Transform". IEEE Transactions on Signal Processing. Vol. 52, Number 5, May 2004, pp. 1304-1314.