То В этом примере показано, как двойной древовидный комплексный дискретный вейвлет преобразовывает (DT-CWT), обеспечивает преимущества перед критически произведенным DWT для сигнала, изображения и обработки объема. Двойной древовидный DWT реализован, когда два разделяют двухканальные наборы фильтров. Получать преимущества описало в этом примере, вы не можете произвольно выбрать масштабирование и фильтры вейвлета, используемые в этих двух деревьях. Lowpass (масштабирование) и highpass (вейвлет) фильтры одного дерева, , должен сгенерировать масштабирующуюся функцию и вейвлет, которые являются аппроксимированными преобразованиями Гильберта масштабирующейся функции и вейвлета, сгенерированного lowpass и highpass фильтрами другого дерева, . Поэтому функции масштабирования с комплексным знаком и вейвлеты, сформированные из этих двух деревьев, приблизительно аналитичны.
В результате двойной древовидный DWT показывает меньше отклонения сдвига и больше направленной селективности, чем критически произведенный DWT с только a фактор сокращения для - размерные данные. Сокращение в двойном древовидном DWT значительно меньше сокращения в неподкошенном (стационарном) DWT.
Этот пример иллюстрирует, что аппроксимированная инвариантность сдвига двойного древовидного DWT, выборочная ориентация вейвлетов анализа двойного дерева в 2D и 3-D, и использование двойного древовидного комплексного вейвлета преобразовывают в шумоподавление объема и изображение.
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 [3], 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(3,1,1) stem(t,details(:,2),'Marker','none','ShowBaseline','off') title('Level 2') ylabel('mV') subplot(3,1,2) stem(t,details(:,3),'Marker','none','ShowBaseline','off') title('Level 3') ylabel('mV') subplot(3,1,3) 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(3,1,1) stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off') title('Level 2') ylabel('mV') subplot(3,1,2) stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off') title('Level 3') ylabel('mV') subplot(3,1,3) 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=7×3 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=7×3 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';'2-D 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 2-D 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(1,2,1) imagesc(dtImage) axis equal axis off colormap gray title('Complex Oriented Dual-Tree') subplot(1,2,2) 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(1,2,1) imagesc(dtImage) axis equal axis off colormap gray title('DT-CWT') subplot(1,2,2) 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 позволяет вам обеспечить направленную селективность в 3-D с минимальным сокращением. В 3-D в двойном древовидном преобразовании существует 28 поддиапазонов вейвлета.
Чтобы продемонстрировать направленную селективность 3-D двойного древовидного вейвлета преобразовывают, визуализируют пример 3-D изоповерхности и 3-D двойного дерева и отделимых вейвлетов DWT. Во-первых, визуализируйте действительные и мнимые части отдельно двух двойных древовидных поддиапазонов.
helperVisualize3D('Dual-Tree',28,'separate');
helperVisualize3D('Dual-Tree',25,'separate');
Красный фрагмент графика изоповерхности указывает на положительное отклонение вейвлета от нуля, в то время как синий обозначает отрицательное отклонение. Можно ясно видеть направленную селективность на пробеле действительных и мнимых частей двойных древовидных вейвлетов. Теперь визуализируйте один из двойных древовидных поддиапазонов с действительными и мнимыми графиками, построенными вместе как одна изоповерхность.
helperVisualize3D('Dual-Tree',25,'real-imaginary');
Предыдущий график демонстрирует, что действительные и мнимые части являются переключенными версиями друг друга на пробеле. Это отражает то, что мнимая часть комплексного вейвлета является аппроксимированным преобразованием Гильберта действительной части. Затем визуализируйте изоповерхность действительного ортогонального вейвлета в 3-D для сравнения.
helperVisualize3D('DWT',7);
Смешивание ориентаций, наблюдаемых в 2D DWT, является еще более явным в 3-D. Так же, как в 2D случае, смешивание ориентаций в 3-D приводит к значительному вызову или блокированию артефактов. Чтобы продемонстрировать это, исследуйте 3-D 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 случаю направленная селективность 3-D 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. Мы продемонстрировали, как эти свойства могут привести к улучшенной производительности в анализе сигнала, представлении ребер в изображениях и объемах, и шумоподавлении объема и изображении.