Инвариантность ближнего сдвига DTCWT

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

Чтобы продемонстрировать это на тестовом сигнале, создайте два сдвинутых дискретно-временных импульса длиной 128 выборок. Один сигнал имеет единичный импульс в выборке 60, в то время как другой сигнал имеет единичный импульс в выборке 64. Оба сигнала явно имеют единичную энергию ($${}^{\mathrm{норма\; l2}}$$).

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 и вычислите энергию (возведенные в квадрат $${}^{\mathrm{}}$$нормы l2) коэффициентов. Постройте график коэффициентов на одной шкале. Сдвиг четырех выборок в сигнале вызвал значительное изменение энергии коэффициентов 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])

Чтобы продемонстрировать полезность аппроксимационной инвариантности сдвига в реальных данных, анализируем сигнал электрокардиограммы (ЭКГ). Интервал дискретизации для сигнала ЭКГ составляет 1/180 секунд. Данные взяты из Percival & Walden [3], стр. 125 (данные первоначально предоставлены Уильямом Константином и Пером Рейнхаллом, Вашингтонский университет). Для удобства мы берем данные, чтобы начать с t = 0.

load wecg
dt = 1/180;
t = 0:dt:(length(wecg)*dt)-dt;
figure
plot(t,wecg)
xlabel('Seconds')
ylabel('Millivolts')

Большие положительные пики с интервалом приблизительно 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 использует разделяемую фильтрацию столбцов и строк изображения. Использовать функцию помощника helperPlotCritSampDWT для построения графика LH, HL и HH вейвлетов для наименее асимметричного фазового вейвлета Daubechies с 4 моментами исчезновения, sym4.

 helperPlotCritSampDWT

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

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

Фильтры, применяемые к столбцам и строкам, могут быть из одной пары фильтров, $$\{{}_{\mathrm{h0}},{}_{\mathrm{}}h1$$} или $${\{}_{\mathrm{}}g0{,}_{\mathrm{}}$$g1}, или из разных пар фильтров, $${}_{\mathrm{\{}}{\mathrm{h0}}_{\mathrm{,}}g1{\}}_{\mathrm{,}}{\{}_{\mathrm{}}g0$$, h1}. Использовать функцию помощникаhelperPlotWaveletDTCWT для построения графика ориентации 12 вейвлетов, соответствующих действительной и мнимой частям DTCWT. Верхняя строка предыдущего рисунка показывает реальные части шести вейвлетов, а вторая строка показывает воображаемые части.

helperPlotWaveletDTCWT

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

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

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 как для стандартного 2-D DWT, так и для комплексного ориентированного DWT с двойным деревом.

helperPlotOctagon(woctagonHyperbolic)

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

Обличение изображения

Благодаря способности выделять различные ориентации в отдельных поддиапазонах, DTCWT часто способен превосходить стандартный разделяемый DWT в таких применениях, как деноизирование изображения. Чтобы продемонстрировать это, используйте функцию помощника helperCompare2DDenoisingDTCWT. Вспомогательная функция загружает изображение и добавляет нулевое среднее значение белого гауссова шума при $$\lambda \mathrm{=}\mathrm{}$$25. Для предоставленного пользователем диапазона пороговых значений функция сравнивает отрицание с использованием мягкого порогового значения для критически отобранного DWT и DTCWT. Для каждого порогового значения отображаются среднеквадратическая ошибка (среднеквадратичная) и пиковое отношение сигнал/шум (PSNR).

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

DTCWT превосходит стандартный DWT в RMS error и PSNR.

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

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

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

Избирательность по направлению в 3-D

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

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

Чтобы продемонстрировать это, рассмотрим набор данных МРТ, состоящий из 16 срезов. Гауссов шум со стандартным отклонением 10 был добавлен к исходному набору данных. Отображение набора шумных данных.

load MRI3D
montage(reshape(noisyMRI,[128 128 1 16]),'DisplayRange',[])

Следует отметить, что исходное SNR до деноизирования составляет приблизительно 11 дБ.

20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))

ans = 11.2997

Деноизируйте набор данных МРТ до уровня 4, используя как DTCWT, так и DWT. В обоих случаях используются одинаковые длины вейвлет-фильтров. Постройте график результирующего SNR как функции порогового значения. Отображение деноизированных результатов для DTCWT и DWT, полученных при лучшем SNR.

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

Ссылки