Пограничные эффекты

Классически DWT задается для последовательностей с длиной некоторой степени 2, и необходимы различные способы расширения выборок других размеров. Методы расширения сигнала включают заполнение нулями, сглаживание заполнения, периодическое расширение и репликацию краевого значения (симметризация).

Основной алгоритм для DWT не ограничивается диадической длиной и основан на простой схеме: свертка и понижающая дискретизация. Как обычно, когда свертка выполняется по сигналам конечной длины, возникают искажения границы.

Расширения сигналов: заполнение нулями, симметризация и сглаживание заполнения

Для борьбы с искажениями границ граница должна рассматриваться по-другому по сравнению с другими частями сигнала.

Для решения этой проблемы доступны различные методы, называемые «вейвлетами на интервале» [1]. Эти интересные конструкции эффективны в теории, но не совсем удовлетворительны с практической точки зрения.

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

Подробная информация об обосновании этих схем содержится в главе 8 книги «Wavelets and Filter Banks» Странга и Нгуена [2].

Доступные режимы расширения сигнала следующие (см. dwtmode):

Заполнение нулями ('zpd'): Этот метод используется в версии DWT, приведенной в предыдущих разделах, и принимает, что сигнал равен нулю вне исходной поддержки.
Недостатком заполнения нулями является то, что на границе искусственно создаются разрывы.
Симметризация ('sym'): Этот метод принимает, что сигналы или изображения могут быть восстановлены вне их исходной поддержки путем симметричной репликации краевых значений.
Это режим по умолчанию вейвлет в тулбоксе.
Симметризация имеет недостаток искусственного создания разрывов первой производной на границе, но этот метод хорошо работает в целом для изображений.
Сглаживайте заполнение порядка 1 ('spd'или 'sp1'): Этот метод принимает, что сигналы или изображения могут быть восстановлены вне их исходной поддержки с помощью простой экстраполяции производной первого порядка: заполнение с помощью линейного расширения для подгонки к первым двум и последним двум значениям.
Сглаживание заполнения работает хорошо в целом для сглаживания сигналов.
Сглаживайте заполнение порядка 0 ('sp0'): Этот метод принимает, что сигналы или изображения могут быть восстановлены вне их исходной поддержки путем простой постоянной экстраполяции. Для расширения сигнала это повторение первого значения слева и последнего значения справа.
Периодическое заполнение (1) ('ppd'): Этот метод принимает, что сигналы или изображения могут быть восстановлены вне их первоначальной поддержки периодическим расширением.
Недостатком периодического заполнения является то, что на границе искусственно создаются разрывы.

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

Периодическое заполнение (2) ('per'): Если длина сигнала нечетная, сигнал сначала расширяется путем добавления дополнительной выборки, равной последнему значению справа. Затем с каждой стороны выполняется минимальное периодическое расширение. Для изображений существует то же правило. Этот режим внутреннего абонента используется для SWT (1-D & 2D).

Этот последний режим производит разложение вейвлет наименьшей длины. Но режим расширения, используемый для IDWT, должен быть тем же самым, чтобы гарантировать идеальную реконструкцию.

Прежде чем рассматривать иллюстративный пример, давайте сравним некоторые свойства теоретического Дискретного Преобразования Вейвлета с фактическим DWT.

Теоретический DWT применяется к сигналам, которые заданы на интервале времени бесконечной длины (Z). Для ортогонального вейвлета это преобразование имеет следующие желательные свойства:

Нормальная консервация
Предположим cA и cD быть приближением и детализацией коэффициентов DWT X сигнала бесконечной длины. Затем l²-норма сохранена:
‖ <reservedrangesplaceholder0> ² =  <reservedrangesplaceholder0> ² +  <reservedrangesplaceholder0> ²
Ортогональность
Позвольте A и D быть восстановленным приближением и деталями. Затем A и D ортогональны и
‖ <reservedrangesplaceholder0> ² =  <reservedrangesplaceholder0> ² +  <reservedrangesplaceholder0> ²
Идеальная реконструкция
X = A + D

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

Обеспечить важнейшее свойство 3 (идеальная реконструкция) для произвольных вариантов

Длина сигнала
Вейвлет
Режим расширения

свойства 1 и 2 могут быть потеряны. Эти свойства верны для расширенного сигнала длины, обычно большей, чем длина исходного сигнала. Так что всегда сохраняется только идеальное свойство реконструкции. Тем не менее, если DWT выполняется с использованием периодического режима расширения ('per') и если длина сигнала делится на 2^J, где J - максимальное разложение уровня, свойства 1, 2 и 3 остаются истинными.

Практические факторы

Открыть скрипт

Итерационный шаг DWT состоит из фильтрации с последующей понижающей дискретизацией:

Примените lowpass, затем понижайте значение на два, чтобы получить коэффициенты приближения.
Примените фильтр highpass, затем понижайте значение на два, чтобы получить коэффициенты детализации.

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

В реальном мире мы имеем дело с сигналами конечной длины. С желанием применения теоретического алгоритма DWT к практическому необходимо рассмотреть вопрос о граничных условиях: как следует расширять сигнал?

Перед исследованием различных сценариев сохраните текущий режим расширения контура.

origmodestatus = dwtmode('status','nodisplay');

Периодический, Степень 2

Рассмотрим следующий пример. Загрузите noisdopp данные. Сигнал имеет 1024 выборки, которая является степенью 2. Использование dwtmode для установки периодического режима внутреннего абонента. Затем используйте wavedec для получения DWT уровня 3 сигнала с помощью ортогональной db4 вейвлет.

load noisdopp;
x = noisdopp;
lev = 3;
wav = 'db4';
dwtmode('per','nodisp')
[c,bk] = wavedec(x,lev,wav);
bk

bk =

         128         128         256         512        1024

Вектор бухгалтерии bk содержит количество коэффициентов по уровням. На каждом этапе количество коэффициентов детализации уменьшается точно в 2 раза. В конце присутствуют $1024/2^3 = 128$ приближения коэффициента.

Сравните $l^2$ -нормы.

fprintf('l2-norm difference: %.5g\n',sum(x.^2)-sum(c.^2))

l2-norm difference: 9.0658e-09

Получите восстановленные приближения и детали путем установки на 0 соответствующих сегментов вектора коэффициентов c и взятия обратного DWT.

cx = c;
cx(bk(1)+1:end) = 0;
reconApp = waverec(cx,bk,wav);
cx = c;
cx(1:bk(1)) = 0;
reconDet = waverec(cx,bk,wav);

Проверяйте ортогональность.

fprintf('Orthogonality difference %.4g\n',...
    sum(x.^2)-(sum(reconApp.^2)+sum(reconDet.^2)))

Orthogonality difference 1.816e-08

Проверяйте на идеальную реконструкцию.

fprintf('Perfect reconstruction difference: %.5g\n',...
    max(abs(x-(reconApp+reconDet))));

Perfect reconstruction difference: 1.674e-11

Сохраняются три теоретических свойства DWT.

Периодический, Не Степень 2

Теперь получите трехуровневый DWT сигнала с 1026 выборками. Используйте тот же вейвлет, и режим расширения выше. Количество коэффициентов на этапе n неравномерно делит длину сигнала.

x = [0 0 noisdopp];
[c,bk] = wavedec(x,lev,wav);
bk

bk =

         129         129         257         513        1026

Проверяйте на $l^2$ -нормовую консервацию, ортогональность и идеальную реконструкцию.

cx = c;
cx(bk(1)+1:end) = 0;
reconApp = waverec(cx,bk,wav);
cx = c;
cx(1:bk(1)) = 0;
reconDet = waverec(cx,bk,wav);

fprintf('l2-norm difference: %.5g\n',sum(x.^2)-sum(c.^2))
fprintf('Orthogonality difference %.4g\n',...
    sum(x.^2)-(sum(reconApp.^2)+sum(reconDet.^2)))
fprintf('Perfect reconstruction difference: %.5g\n',...
    max(abs(x-(reconApp+reconDet))));

l2-norm difference: -1.4028
Orthogonality difference -0.3319
Perfect reconstruction difference: 1.6858e-11

Совершенная реконструкция удовлетворена, но $l^2$ -норма и ортогональность не сохранились.

Не Периодический, Степень 2

Получите трехуровневый DWT сигнала с 1024. Используйте тот же вейвлет, что и выше, но на этот раз измените режим расширения, чтобы сгладить расширение порядка 1. Количество коэффициентов на этапе n неравномерно делит длину сигнала.

dwtmode('sp1','nodisp')
[c,bk] = wavedec(x,lev,wav);
bk

bk =

         134         134         261         516        1026

Проверяйте на $l^2$ -нормовую консервацию, ортогональность и идеальную реконструкцию.

cx = c;
cx(bk(1)+1:end) = 0;
reconApp = waverec(cx,bk,wav);
cx = c;
cx(1:bk(1)) = 0;
reconDet = waverec(cx,bk,wav);

fprintf('l2-norm difference: %.5g\n',sum(x.^2)-sum(c.^2))
fprintf('Orthogonality difference %.4g\n',...
    sum(x.^2)-(sum(reconApp.^2)+sum(reconDet.^2)))
fprintf('Perfect reconstruction difference: %.5g\n',...
    max(abs(x-(reconApp+reconDet))));

l2-norm difference: -113.58
Orthogonality difference -2.678
Perfect reconstruction difference: 1.6372e-11

Опять же, удовлетворена только идеальная реконструкция.

Восстановите исходный режим расширения.

dwtmode(origmodestatus,'nodisplay');

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

Система координат является набором функций $\{\phi_k\}$ , которые удовлетворяют следующему условию: существуют такие константы $0 < A \leq B$ , что для любой функции $f$ удерживается неравенство систем координат: $A{\|f\|}^2 \leq \Sigma_k {|\langle f,\phi_k\rangle|}^2 \leq B {\|f\|}^2.$

Функции в системе координат обычно не являются линейно независимыми. Это означает, что функция $f$ не имеет уникального расширения в. $\phi_k$ Daubechies [3] показывает, что если является $\{\tilde{\phi}_k\}$ двойственной системой координат и для $f = \Sigma_{k\in K} c_k\phi_k$ некоторых, а $c = \left(c_k\right)\in l^2(K)$ если не все равны $c_k$ , то. $\langle f,\tilde{\phi}_k\rangle$ $\Sigma_{k\in K}{|c_k|}^2 \geq \Sigma_{k\in K}{|\langle f,\tilde{\phi}_k\rangle|}^2$

Если, $A = B$ система координат называется плотной системой координат. Если $A = B = 1$ и $\|\phi_k\|^2 = 1$ для всех, $\phi_k$ система координат является ортонормированным базисом. Если, $A \neq B$ то энергия не обязательно сохраняется, и общее количество коэффициентов может превысить длину сигнала. Если режим расширения периодический, вейвлет ортогональен, и длина сигнала делится на, $2^J$ где является $J$ максимальным уровнем разложения вейвлета, все три теоретических свойства DWT удовлетворяются.

Произвольные расширения

Открыть Live Script

Интересно заметить, что, если используется произвольное расширение и разложение выполняется с помощью схемы свертки-понижающей дискретизации, совершенная реконструкция восстанавливается с помощью idwt или idwt2.

Создайте сигнал и получите фильтры, сопоставленные со db9 вейвлет.

x = sin(0.3*[1:451]);
w = 'db9';
[LoD,HiD,LoR,HiR] = wfilters(w);

Добавьте и подготовьте length(LoD) случайные числа к сигналу. Постройте график исходных и расширенных сигналов.

lx = length(x);
lf = length(LoD);
ex = [randn(1,lf) x randn(1,lf)];
ymin = min(ex);
ymax = max(ex);
subplot(2,1,1)
plot(lf+1:lf+lx,x)
axis([1 lx+2*lf ymin ymax]);
title('Original Signal')
subplot(2,1,2)
plot(ex)
title('Extended Signal')
axis([1 lx+2*lf ymin ymax])

Figure contains 2 axes. Axes 1 with title Original Signal contains an object of type line. Axes 2 with title Extended Signal contains an object of type line.

Используйте lowpass и highpass вейвлета фильтры разложения, чтобы получить одноуровневый вейвлет разложение расширенного сигнала.

la = floor((lx+lf-1)/2);
ar = wkeep(dyaddown(conv(ex,LoD)),la);
dr = wkeep(dyaddown(conv(ex,HiD)),la);

Подтвердите идеальную реконструкцию сигнала.

xr = idwt(ar,dr,w,lx);
err0 = max(abs(x-xr))

err0 = 5.4700e-11

Сравнение различий в расширениях

Теперь давайте проиллюстрировать различия между первыми тремя методами как для 1-D, так и для 2-D сигналов.

Заполнение нулями

Используя приложение Wavelet Analysis, мы исследуем эффекты заполнения нули.

Из MATLAB^® приглашение, тип
```
dwtmode('zpd')
```
Из подсказки MATLAB введите waveletAnalyzer.
Появится Wavelet Analyzer.
Выберите меню Вейвлета 1-D элемента. Появляется инструмент дискретного вейвлет для 1-D данных о сигнале.
В меню File выберите параметр Example Analysis и выберите Basic Signals > with db2 at level 5 > Two nearly disprinuities.
Выберите Режим отображения: Показать и прокрутить.
Коэффициенты детализации четко показывают эффекты конца сигнала.
Симметричное расширение
Из подсказки MATLAB введите
```
dwtmode('sym')
```
Выберите меню Вейвлета 1-D элемента.
Появляется инструмент дискретного вейвлет для 1-D данных о сигнале.
В меню File выберите параметр Example Analysis и выберите Basic Signals > with db2 at level 5 > Two nearly disprinuities .
Выберите Режим отображения: Показать и прокрутить.
Коэффициенты детализации четко показывают эффекты конца сигнала.
Сглаживайте заполнение
Из подсказки MATLAB введите
```
dwtmode('spd')
```
Выберите меню Вейвлета 1-D элемента.
Появляется инструмент дискретного вейвлет для 1-D данных о сигнале.
В меню File выберите параметр Example Analysis и выберите Basic Signals > with db2 at level 5 > Two nearly disprinuities.
Выберите Режим отображения: Показать и прокрутить.
Коэффициенты детализации показывают, что эффекты конца сигнала не присутствуют, и разрывы хорошо детектируются.

Расширения изображений

Открыть Live Script

Давайте рассмотрим пример изображения. Сохраните текущий режим расширения. Загрузите и отобразите geometry изображение.

origmodestatus = dwtmode('status','nodisplay');
load geometry
nbcol = size(map,1);
colormap(pink(nbcol))
image(wcodemat(X,nbcol))

Figure contains an axes. The axes contains an object of type image.

Заполнение нулями

Установите режим расширения в нулевое заполнение и выполните разложение изображения до уровня 3 с помощью sym4 вейвлет. Затем восстановите приближение уровня 3.

lev = 3;
wname = 'sym4';
dwtmode('zpd','nodisp')
[c,s] = wavedec2(X,lev,wname);
a = wrcoef2('a',c,s,wname,lev);
image(wcodemat(a,nbcol))

Figure contains an axes. The axes contains an object of type image.

Симметричное расширение

Установите симметричный режим расширения и выполните разложение изображения до уровня 3 с помощью sym4 вейвлет. Затем восстановите приближение уровня 3.

dwtmode('sym','nodisp')
[c,s] = wavedec2(X,lev,wname);
a = wrcoef2('a',c,s,wname,lev);
image(wcodemat(a,nbcol))

Figure contains an axes. The axes contains an object of type image.

Сглаживайте заполнение

Установите режим расширения для сглаживания заполнения и выполните разложение изображения до уровня 3 с помощью sym4 вейвлет. Затем восстановите приближение уровня 3.

dwtmode('spd','nodisp')
[c,s] = wavedec2(X,lev,wname);
a = wrcoef2('a',c,s,wname,lev);
image(wcodemat(a,nbcol))

Figure contains an axes. The axes contains an object of type image.

Восстановите исходный режим расширения.

dwtmode(origmodestatus,'nodisplay')

Ссылки

[1] Cohen, A., I. Daubechies, B. Jawerth, and P. Vial. «Мультирезолюционный анализ, вейвлеты и быстрые алгоритмы на интервале». Comptes Rendus Academ. Science Paris Sér. А, том 316, стр. 417-421, 1993.

[2] Странг, Г. и Т. Нгуен. Вейвлеты и банки фильтров. Wellesley, MA: Wellesley-Cambridge Press, 1996.

[3] Daubechies, I. Десять лекций по вейвлетам, серия региональных конференций CBMS-NSF по прикладной математике. Филадельфия, Пенсильвания: СИАМ Эд, 1992.

См. также

dwtfilterbank | wavedec

Документация