Около инвариантности сдвига DTCWT

DWT страдает от отклонения сдвига, означая, что маленькие сдвиги во входном сигнале или изображении могут вызвать существенные изменения в распределении энергии сигнала/изображения через шкалы в коэффициентах DWT. DTCWT является приблизительно инвариантом сдвига.

Чтобы продемонстрировать это на тестовом сигнале, создайте два переключенных импульса дискретного времени 128 выборок в длине. Один сигнал имеет модульный импульс в демонстрационных 60, в то время как другой сигнал имеет модульный импульс в демонстрационных 64. Оба сигнала ясно имеют модульную энергию ($$l^2$$ норма.

kronDelta1 = zeros(128,1);
kronDelta1(60) = 1;
kronDelta2 = zeros(128,1);
kronDelta2(64) = 1;

Установите дополнительный режим DWT на периодический. Получите DWT и DTCWT двух сигналов вниз к уровню 3 с вейвлетом и масштабирующимися фильтрами длины 14. Извлеките коэффициенты детали уровня 3 для сравнения.

origmode = dwtmode('status','nodisplay');
dwtmode('per','nodisp')
J = 3;
[dwt1C,dwt1L] = wavedec(kronDelta1,J,'sym7');
[dwt2C,dwt2L] = wavedec(kronDelta2,J,'sym7');
dwt1Cfs = detcoef(dwt1C,dwt1L,3);
dwt2Cfs = detcoef(dwt2C,dwt2L,3);

[dt1A,dt1D] = dualtree(kronDelta1,'Level',J,'FilterLength',14);
[dt2A,dt2D] = dualtree(kronDelta2,'Level',J,'FilterLength',14);
dt1Cfs = dt1D{3};
dt2Cfs = dt2D{3};

Постройте абсолютные значения DWT и коэффициентов DTCWT для двух сигналов на уровне 3 и вычислите энергию (в квадрате $$l^2$$ нормы) коэффициентов. Постройте коэффициенты по той же шкале. Четыре выборки переключаются на нижний регистр, сигнал вызвал существенное изменение в энергии уровня 3 коэффициенты DWT. Энергия на уровне 3 коэффициенты DTCWT изменилась только на приблизительно 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)
ylim([0 0.4])
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'DWT';['Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])

figure
subplot(1,2,1)
stem(abs(dt1Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree CWT';['Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title({'Dual-tree CWT';['Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)]},...
    'fontsize',10)
ylim([0 0.4])

Чтобы продемонстрировать утилиту аппроксимированной инвариантности сдвига в действительных данных, мы анализируем электрокардиограмму (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 изолируются наиболее заметно в тех шкалах для данной частоты дискретизации.

figure
J = 6; 
[df,rf] = dtfilters('farras');
[dtDWT1,L1] = wavedec(wecg,J,df(:,1),df(:,2));
details = zeros(2048,3);
details(2:4:end,2) = detcoef(dtDWT1,L1,2);
details(4:8:end,3) = detcoef(dtDWT1,L1,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.

[dtcplxA,dtcplxD] = dualtree(wecg,'Level',J,'FilterLength',14);
details = zeros(2048,3);
details(2:4:end,2) = dtcplxD{2};
details(4:8:end,3) = dtcplxD{3};
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,L2] = wavedec(wecgShift,J,df(:,1),df(:,2));

detCfs1 = detcoef(dtDWT1,L1,1:J,'cells');
apxCfs1 = appcoef(dtDWT1,L1,rf(:,1),rf(:,2),J);
cfs1 = horzcat(detCfs1,{apxCfs1});
detCfs2 = detcoef(dtDWT2,L2,1:J,'cells');
apxCfs2 = appcoef(dtDWT2,L2,rf(:,1),rf(:,2),J);
cfs2 = horzcat(detCfs2,{apxCfs2});

sigenrgy = norm(wecg,2)^2;
enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0));
enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'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%-е изменения в энергии между исходным и переключенным сигналом. Затем мы повторяем, что этот анализ с помощью комплексного двойного древовидного дискретного вейвлета преобразовывает.

[dtcplx2A,dtcplx2D] = dualtree(wecgShift,'Level',J,'FilterLength',14);
dtCfs1 = vertcat(dtcplxD,{dtcplxA});
dtCfs2 = vertcat(dtcplx2D,{dtcplx2A});

dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtCfs1,'uni',0));
dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtCfs2,'uni',0));
dtenr1 = dtenr1(:);
dtenr2 = dtenr2(:);
table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})

ans=7×3 table
    Level     dtenr1    dtenr2
    ______    ______    ______

    {'D1'}    5.3533    5.3619
    {'D2'}    6.2672    6.2763
    {'D3'}    12.155     12.19
    {'D4'}    8.2831     8.325
    {'D5'}      5.81    5.8577
    {'D6'}    3.1768    3.0526
    {'A6'}    58.403    58.384

Двойное древовидное преобразование производит сопоставимый дисперсионный анализ шкалой для исходного сигнала и его циркулярной переключенной версии.

Направленная селективность в обработке изображений

Стандартная реализация DWT в 2D использовании отделимая фильтрация столбцов и строки изображения. Используйте функцию помощника helperPlotCritSampDWT построить LH, HL и вейвлеты HH для Добечиса меньше всего - асимметричный вейвлет фазы с 4 исчезающими моментами, sym4.

 helperPlotCritSampDWT

Обратите внимание на то, что LH и вейвлеты HL имеют ясные горизонтальные и вертикальные ориентации соответственно. Однако вейвлет HH на ультраправых смесях и +45 и-45 направлений степени, производя артефакт шахматной доски. Это смешивание ориентаций происходит из-за использования отделимых фильтров с действительным знаком. HH отделимый фильтр с действительным знаком имеет полосы пропускания во всех четырех высокочастотных углах 2D плоскости частоты.

Двойной древовидный DWT достигает направленной селективности при помощи вейвлетов, которые приблизительно аналитичны, означая, что у них есть поддержка только на одной половине оси частоты. В двойном древовидном DWT существует шесть поддиапазонов и для действительных и для мнимых частей. Эти шесть действительных частей сформированы путем добавления выходных параметров фильтрации столбца, сопровождаемой фильтрацией строки входного изображения в этих двух деревьях. Эти шесть мнимых частей сформированы путем вычитания выходных параметров фильтрации столбца, сопровождаемой фильтрацией строки.

Фильтры применились к столбцам, и строки могут быть от той же пары фильтра, $$\{h_0,h_1\}$$ или $$\{g_0,g_1\}$$, или от различных пар фильтра, $$\{h_0,g_1\},\{g_0,h_1\}$$. Используйте функцию помощника helperPlotWaveletDTCWT построить ориентацию этих 12 вейвлетов, соответствующих действительным и мнимым частям DTCWT. Верхняя строка предыдущей фигуры показывает действительные части этих шести вейвлетов, и вторая строка показывает мнимые части.

helperPlotWaveletDTCWT

Представление ребра в двух измерениях

Аппроксимированная аналитичность и выборочная ориентация комплексных двойных древовидных вейвлетов обеспечивают наилучшее решение по стандартному 2D DWT в представлении ребер в изображениях. Чтобы проиллюстрировать это, мы анализируем тестовые изображения с ребрами, состоящими из линии и сингулярности кривой в нескольких направлениях с помощью критически произведенного 2D DWT, и 2D комплекс ориентировался, двойное дерево преобразовывают. Во-первых, анализируйте изображение восьмиугольника, который состоит из сингулярности линии.

load woctagon
imagesc(woctagon)
colormap gray
title('Original Image')
axis equal
axis off

Используйте функцию помощника helprPlotOctagon чтобы анализировать изображение вниз к уровню 4 и восстановить приближение изображений на основе уровня 4 детализируют коэффициенты.

helperPlotOctagon(woctagon)

Затем анализируйте восьмиугольник с гиперболическими ребрами. Ребра в гиперболическом восьмиугольнике являются сингулярностью кривой.

load woctagonHyperbolic
figure
imagesc(woctagonHyperbolic)
colormap gray
title('Octagon with Hyperbolic Edges')
axis equal
axis off

Снова, используйте функцию помощника helprPlotOctagon анализировать изображение вниз к уровню 4 и восстановить приближение изображений на основе коэффициентов детали уровня 4 и для стандартного 2D DWT и для комплекса ориентировали двойной древовидный DWT.

helperPlotOctagon(woctagonHyperbolic)

Обратите внимание на то, что звонящие артефакты, очевидные в 2D критически, произвели DWT, отсутствуют в 2D DTCWT обоих изображений. DTCWT более искренне воспроизводит сингулярность кривой и линия.

Отобразите шумоподавление

Из-за способности изолировать отличные ориентации в отдельных поддиапазонах, DTCWT часто может превзойти стандартный отделимый DWT по характеристикам в приложениях как шумоподавление изображений. Чтобы продемонстрировать это, используйте функцию помощника helperCompare2DDenoisingDTCWT. Функция помощника загружает изображение и добавляет нулевой средний белый Гауссов шум с $$\sigma =25$$. Для предоставленной пользователями области значений порогов функция сравнивает шумоподавление с помощью мягкой пороговой обработки для критически произведенного DWT и DTCWT. Для каждого порогового значения отображены среднеквадратичная (RMS) ошибка и пиковое отношение сигнал-шум (PSNR).

numex = 3;
helperCompare2DDenoisingDTCWT(numex,0:2:100,'PlotMetrics');

DTCWT превосходит стандартный DWT по характеристикам по ошибке RMS и PSNR.

Затем получите изображения denoised для порогового значения 25, который равен стандартному отклонению аддитивного шума.

numex = 3;
helperCompare2DDenoisingDTCWT(numex,25,'PlotImage');

С пороговым значением, равным стандартному отклонению аддитивного шума, DTCWT обеспечивает PSNR на почти 4 дБ выше, чем стандартный 2D DWT.

Направленная селективность в 3-D

Звонящие артефакты, наблюдаемые с отделимым DWT в двух измерениях, усилены при расширении анализа вейвлета к более высоким размерностям. DTCWT позволяет вам обеспечить направленную селективность в 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 и детали вейвлета DTCWT сферического объема. Сфера 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('DTCWT 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')

Увеличьте масштаб изображений и в DTCWT и в монтаже DWT, и вы будете видеть, как видный блокирующиеся артефакты в деталях DWT сравниваются с DTCWT.

Шумоподавление объема

Подобно 2D случаю направленная селективность 3-D DTCWT часто приводит к улучшениям шумоподавления объема.

Чтобы продемонстрировать это, рассмотрите набор данных 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 с помощью и DTCWT и DWT. Подобные длины фильтра вейвлета используются в обоих случаях. Постройте получившийся ОСШ в зависимости от порога. Отобразите denoised результаты и для DTCWT и для DWT, полученного в лучшем ОСШ.

[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI);

figure
montage(reshape(imrecDTCWT,[128 128 1 16]),'DisplayRange',[])
title('DTCWT Denoised Volume')

figure
montage(reshape(imrecDWT,[128 128 1 16]),'DisplayRange',[])
title('DWT Denoised Volume')

Восстановите исходный дополнительный режим.

dwtmode(origmode,'nodisplay')

Сводные данные

Мы показали, что двойной древовидный CWT обладает желательными свойствами близкой инвариантности сдвига и направленной селективности, не достижимой с критически произведенным DWT. Мы продемонстрировали, как эти свойства могут привести к улучшенной производительности в анализе сигнала, представлении ребер в изображениях и объемах, и шумоподавлении объема и изображении.

Ссылки