Инвариация ближнего сдвига 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$$ норм) коэффициентов. Постройте график коэффициентов в той же шкале. Сдвиг четырех выборок в сигнале вызвал значительное изменение энергии коэффициентов DWT уровня 3. Энергия в коэффициентах DTCWT уровня 3 изменилась только приблизительно на 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 составляет 1/180 секунд. Данные взяты из Percival & Walden [3], стр. 25 (данные первоначально предоставлены Уильямом Константином и Пером Рейнхоллом, Университет Вашингтона). Для удобства мы принимаем данные, чтобы начать с 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')

Как критически выбранные, так и двухдревовидные вейвлет локализуют важную функцию формы волны ЭКГ в сходные шкалы

Важным применением вейвлетов в 1-D сигналах является получение анализа отклонения по шкале. Это означает, что этот дисперсионный анализ не должен быть чувствительным к круговым сдвигам в входном сигнале. К сожалению, это не так с критически выбранным DWT. Чтобы продемонстрировать это, мы круто сдвигаем сигнал ЭКГ на 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 в 2-D использует разделяемую фильтрацию столбцов и строк изображения. Используйте функцию helper helperPlotCritSampDWT построение графика LH, HL и HH вейвлетов для наименее асимметричнейшая фаза вейвлета Daubechies с 4 исчезающими моментами, sym4.

 helperPlotCritSampDWT

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

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

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

helperPlotWaveletDTCWT

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

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

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

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

helperPlotOctagon(woctagon)

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

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

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

helperPlotOctagon(woctagonHyperbolic)

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

Шумоподавление изображений

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

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

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

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

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

При пороге значении, равном стандартному отклонению аддитивного шума, DTCWT обеспечивает PSNR почти на 4 дБ выше, чем стандартный 2-D 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)

Смешение ориентаций, наблюдаемое в 2-D DWT, ещё более выражено в 3-D. Так же, как и в 2-D случае, смешение ориентаций в 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.

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

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

Чтобы продемонстрировать это, рассмотрим набор данных МРТ, состоящий из 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

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

Ссылки