Алгоритм быстрого Вейвлета преобразования (FWT)

В 1988 году Маллат выпустил алгоритм быстрого вейвлета разложения и реконструкции [1]. Алгоритм Mallat для дискретного вейвлет (DWT) является, по сути, классической схемой в сообществе обработки сигналов, известной как двухканальный поддиапазонный кодер, использующий сопряженные квадратурные фильтры или квадратурные зеркальные фильтры (QMF).

Алгоритм разложения начинается с сигнала s, затем вычисляет координаты A 1 и _D 1, а затем координаты A 2 и D 2 и так далее.
Алгоритм реконструкции, называемый обратным дискретным вейвлет (IDWT), начинается с координат _AJ и _DJ, затем вычисляет координаты _AJ _-1, а затем с помощью координат _AJ -1 и _DJ -1 вычисляет координаты _AJ -2 и так далее.

В этом разделе рассматриваются следующие темы:

Фильтры, используемые для вычисления DWT и IDWT

Для ортогонального вейвлета, в среде мультиразрешения, мы начинаем с функции масштабирования, и функции вейвлета. Одним из фундаментальных отношений является отношение двойной шкалы (уравнение расширения или уравнение уточнения):

$\frac{1}{2} ϕ (\frac{x}{2}) = \sum_{n \in Z} w_{n} ϕ (x - n)$

Все фильтры, используемые в DWT и IDWT, тесно связаны с последовательностью

_{(<reservedrangesplaceholder1>)} n∊Z

Очевидно, что, если и поддерживается компактно, последовательность (_wn) является конечной и может рассматриваться как фильтр. Фильтрующая W, которая называется масштабирующим фильтром (nonnormalized), является

Конечная импульсная характеристика (КИХ)
Длины 2 N
Из суммы 1
Нормы $\frac{1}{\sqrt{2}}$
Нормы 1
Фильтр lowpass

Для примера, для db3 масштабирующий фильтр,

load db3 
db3
    db3 =
    0.2352   0.5706   0.3252  -0.0955  -0.0604   0.0249

sum(db3)
    ans =
         1.0000

norm(db3)
    ans =
         0.7071

Из фильтра W задаем четыре КИХ-фильтра, длиной 2 N и нормой 1, организованные следующим образом .

Фильтры	Низкочастотные	High-Pass
Разложение	Lo_D	Hi_D
Реконструкция	Lo_R	Hi_R

Четыре фильтра вычисляются с помощью следующей схемы.

где qmf таковы, что Hi_R и Lo_R являются квадратурными зеркальными фильтрами (т.е. Hi_R (k) = (-1 )^kLo_R (2 N + 1 - k)) для k = 1, 2,..., 2 N.

Обратите внимание, что wrev разворачивает коэффициенты фильтра. Так Hi_D и Lo_D также являются квадратурными зеркальными фильтрами. Расчет из этих фильтров выполняется с помощью orthfilt. Далее мы иллюстрируем эти свойства с db6 вейвлет.

Загрузите фильтр экстремального масштабирования фазы Daubechies и постройте график коэффициентов.

load db6;
subplot(421); stem(db6,'markerfacecolor',[0 0 1]);
title('Original scaling filter');

Использовать orthfilt для возврата фильтров анализа (разложения) и синтеза (реконструкции).

Получите дискретные преобразования Фурье (DFT) фильтров lowpass и highpass анализа. Постройте график модуля ДПФ.

LoDFT = fft(Lo_D,64);
HiDFT = fft(Hi_D,64);
freq = -pi+(2*pi)/64:(2*pi)/64:pi;
subplot(427); plot(freq,fftshift(abs(LoDFT)));
set(gca,'xlim',[-pi,pi]); xlabel('Radians/sample');
title('DFT Modulus - Lowpass Filter')
subplot(428); plot(freq,fftshift(abs(HiDFT)));
set(gca,'xlim',[-pi,pi]); xlabel('Radians/sample');
title('Highpass Filter');

Алгоритмы

Учитывая s сигнала длины N, DWT состоит не более чем из N каскадов _log2. Начиная с s, первый шаг производит два набора коэффициентов: коэффициенты аппроксимации cA 1 и коэффициенты детализации _cD 1. Эти векторы получаются путем свертки s с Lo_D фильтра нижних частот для аппроксимации и с Hi_D фильтра высоких частот для детализации с последующей диадической децимацией.

Точнее, первый шаг -

Длина каждого фильтра равна 2 L. Результат свертки сигнала N длиной 2 L фильтра - N + 2 L -1. Поэтому F и G сигналов имеют длину N + 2 L - 1. После понижающей дискретизации на 2 векторы коэффициентов cA ₁ и cD 1 имеют длину

$⌊ \frac{N - 1}{2} + L ⌋ .$

Следующий шаг разделяет коэффициенты аппроксимации cA 1 на две части с помощью одной и той же схемы, заменяя s на cA 1 и получая _cA 2 и cD 2 и так далее.

Таким образом, вейвлет разложение s сигнала, анализируемого на уровне j, имеет следующую структуру: [_cAj, _cDj,..., cD ₁].

Эта структура содержит для J = 3 терминальные узлы следующего дерева.

С другой стороны, начинаясь с <reservedrangesplaceholder5> <reservedrangesplaceholder4> и <reservedrangesplaceholder3> <reservedrangesplaceholder2>, IDWT восстанавливает <reservedrangesplaceholder1> <reservedrangesplaceholder0>-1, инвертируя шаг разложения, вставляя нули и скручивая результаты с фильтрами реконструкции.
Для изображений аналогичный алгоритм возможен для двумерных вейвлетов и функций масштабирования, полученных из 1-D вейвлетов тензорным продуктом.
Этот вид 2-D DWT приводит к разложению коэффициентов аппроксимации на уровне j в четырех компонентах: приближение на уровне j + 1 и детали в трех ориентациях (горизонтальной, вертикальной и диагональной).
Следующие графики описывают основные шаги разложения и реконструкции для изображений.

Так, для J = 2 2-D вейвлет имеет следующую форму.

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

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

Почему такой алгоритм существует?

Предыдущий абзац описывает алгоритмы, разработанные для сигналов или изображений конечной длины. Чтобы понять обоснование, мы должны рассмотреть сигналы бесконечной длины. Методы расширения заданного сигнала конечной длины описаны в Border Effects.

Давайте обозначим h = Lo_R и g = Hi_R и сосредоточимся на 1-D случае.

Сначала мы обосновываем, как перейти от уровня j к уровню j + 1, для вектора приближения. Это основной шаг алгоритма разложения для расчета приближений. Детали вычисляются таким же образом с помощью g фильтра вместо h фильтра.

Предположим (_Ak^(j)) k ∊ Z быть координатами _{вектора} Aj:

$A_{j} = \sum_{k} A_{k}^{(j)} ϕ_{j, k}$

и _Ak^(j+1) координаты вектора _Aj ₊ 1:

$A_{j + 1} = \sum_{k} A_{k}^{(j + 1)} ϕ_{j + 1, k}$

_Ak^(j+1) вычисляется по формуле

$A_{k}^{(j + 1)} = \sum_{n} h_{n - 2 k} A_{n}^{(j)}$

Эта формула напоминает формулу свертки.

Расчет очень простое.

Давайте определим

$\tilde{h} (k) = h (- k), и F_{k}^{(j + 1)} = \sum_{n} {\tilde{h}}_{k - n} A n^{(j)} .$

Последовательность F^(j+1) - отфильтрованный выход последовательности A^(j) фильтром $\tilde{h}$ .

Получаем

_Ak^(j+1) = <reservedrangesplaceholder1> 2 <reservedrangesplaceholder0>^(j+1)

Мы должны взять четные значения индекса F. Это понижающая дискретизация.

Последовательность A^(j+1) - версия последовательности с понижающей дискретизацией F^(j+1).

Инициализация осуществляется с помощью _Ak⁽⁰⁾ = s (k), где s (k) - значение сигналов в момент k.

Существует несколько причин для этого удивительного результата, все из которых связаны с ситуацией мультиразрешения и с некоторыми свойствами функций и _j,k _j,k .

Давайте опишем некоторые из них.

Семейство $(ϕ_{0, k}, k \in Z)$ формируется из ортонормальных функций. Как следствие для любого j, семейство $(ϕ_{j, k}, k \in Z)$ ортонормальный.
Семейство с двойным индексом

$(ψ_{j, k}, j \in Z, k \in Z)$
ортонормальный.
Для любого j, $(ϕ_{j, k}, k \in Z)$ являются ортогональными $(ψ_{j^{'}, k}, j^{'} \leq j, k \in Z)$ .
Между двумя последовательными шкалами мы имеем фундаментальное отношение, называемое отношением двойной шкалы.

Двойное соотношение шкалы для $ϕ$

$ϕ_{1, 0} = \sum_{k \in Z} h_{k} ϕ_{0, k}$

$ϕ_{j + 1, 0} = \sum_{k \in Z} h_{k} ϕ_{j, k}$

Это отношение представляет h фильтр алгоритма $(h_{n} = \sqrt{2 ω_{n}})$ . Для получения дополнительной информации см. раздел «Фильтры, используемые для вычисления DWT и IDWT».
Проверяем, что:
1. Координата и _j, 1,0 _наj,k , hk и не зависит от j.
2. Координата _{φ <reservedrangesplaceholder2> +1} , n на φ <reservedrangesplaceholder0> равна $〈 ϕ_{j + 1, n}, ϕ_{j, k} 〉 = h_{k - 2 n}$ .
Эти отношения обеспечивают ингредиенты для алгоритма.
До сих пор мы использовали фильтр h. Высокочастотный фильтрующий g используется в двойном соотношении шкал, связывающем функции Между двумя последовательными шкалами у нас есть следующее двухкомпонентное фундаментальное соотношение.

Двухкомпонентное соотношение между $ψ$ и $ϕ$

$ψ_{1, 0} = \sum_{k \in Z} g_{k} ϕ_{0, k}$

$ψ_{j + 1, 0} = \sum_{k \in Z} g_{k} ϕ_{j, k}$
После шага разложения мы обосновываем теперь алгоритм реконструкции, построив его. Давайте упростим обозначение. Начиная с A 1 и _D 1, давайте изучим A 0 ₌ A 1 ₊ Dj 1. Процедура такая же, чтобы вычислить _A ₌ Aj ₊ ₁ + Dj + 1.
Давайте определим _{α <reservedrangesplaceholder1>} , δ <reservedrangesplaceholder0>, $α_{k}^{0}$ около
$A_{1} = \sum_{n} a_{n} ϕ_{1, n} D_{1} = \sum_{n} δ_{n} ψ_{1, n} A_{0} = \sum_{k} a_{k}^{0} ϕ_{0, k}$
Давайте оценим $α_{k}^{0}$ координаты как
$\begin{array}{l} a_{k}^{0} = 〈 A_{0}, ϕ_{0, k} 〉 = 〈 A_{1} + D_{1}, ϕ_{0, k} 〉 = 〈 A_{1}, ϕ_{0, k} 〉 + 〈 D_{1}, ϕ_{0, k} 〉 \\ = \sum_{n} a_{n} 〈 ϕ_{1, n}, ϕ_{0, k} 〉 + \sum_{n} δ_{n} 〈 ψ_{1, n}, ϕ_{0, k} 〉 \\ = \sum_{n} a_{n} h_{k - 2 n} + \sum_{n} δ_{n} g_{k - 2 n} \end{array}$

Двойное соотношение шкалы для $ϕ$
$ϕ_{1, 0} = \sum_{k \in Z} h_{k} ϕ_{0, k}$	$ϕ_{j + 1, 0} = \sum_{k \in Z} h_{k} ϕ_{j, k}$

Двухкомпонентное соотношение между $ψ$ и $ϕ$
$ψ_{1, 0} = \sum_{k \in Z} g_{k} ϕ_{0, k}$	$ψ_{j + 1, 0} = \sum_{k \in Z} g_{k} ϕ_{j, k}$

Мы сосредоточим наше исследование на первой сумме $\sum_{n} a_{n} h_{k - 2 n};$ вторая сумма $\sum_{n} δ_{n} g_{k - 2 n}$ обрабатывается аналогичным образом.

Вычисления легко организовать, если заметим, что (взятие k = 0 в предыдущих формулах, делает вещи проще)

$\begin{array}{l} \sum_{n} a_{n} h_{- 2 n} = \dots + a_{- 1} h_{2} + a_{0} h_{0} + a_{1} h_{- 2} + a_{2} h_{- 4} + \dots \\ = \dots + a_{- 1} h_{2} + 0 h_{1} + a_{0} h_{0} + 0 h_{- 1} + a_{1} h_{- 2} + 0 h_{- 3} + a_{2} h_{- 4} + \dots \end{array}$

Если мы преобразуем $(α_{n})$ последовательность в новую последовательность $({\tilde{α}}_{n})$ определяется как

_{..., α–1, 0, α0, 0, α1, 0, α2, 0, ...} то есть именно

${\tilde{a}}_{2 n} = a_{n}, {\tilde{a}}_{2 n} + 1 = 0$

Тогда

$\sum_{n} a_{n} h_{- 2 n} = \sum_{n} {\tilde{a}}_{n} h_{- n}$

и продлением

$\sum_{n} a_{n} h_{k - 2 n} = \sum_{n} {\tilde{a}}_{n} h_{k - n}$

С тех пор

$a_{k}^{0} = \sum_{n} {\tilde{a}}_{n} h_{k - n} + \sum_{n} {\tilde{δ}}_{n} g_{k - n}$

этапы реконструкции:

Замените α и δ последовательности сверхдискретизированными версиями α ˜ и $\tilde{δ}$ вставка нулей.
Фильтрация по h и g соответственно.
Суммируйте полученные последовательности.

1-D вейвлет

Основные 1-D объекты

	Объекты	Описание
Сигнал в исходном времени	s _Ak, 0 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0> _Dk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	Исходный сигнал Приближение на уровне k Детализация на уровне k
Коэффициенты во времени, связанном с масштабом	_cAk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0> _cDk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0> [_cAj, _cDj..., <reservedrangesplaceholder0> 1]	Приближения на уровне k Коэффициенты детализации на уровне k Вейвлет на уровне j, j ≥ 1

Объекты

Описание

Сигнал в исходном времени

_Ak, 0 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>

_Dk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>

Исходный сигнал

Приближение на уровне k

Детализация на уровне k

Коэффициенты во времени, связанном с масштабом

_cAk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>

_cDk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>

[_cAj, _cDj..., <reservedrangesplaceholder0> 1]

Приближения на уровне k

Коэффициенты детализации на уровне k

Вейвлет на уровне j, j ≥ 1

Возможности разложения и анализа

Цель	Вход	Выход	Файл
Одноуровневое разложение	s	<reservedrangesplaceholder1> 1, <reservedrangesplaceholder0> 1	`dwt`
Одноуровневое разложение	_cAj	_{<reservedrangesplaceholder1> +1}, <reservedrangesplaceholder0> +1	`dwt`
Разложение	s	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	`wavedec`

Возможности синтеза-реконструкции

Цель	Вход	Выход	Файл
Одноуровневая реконструкция	<reservedrangesplaceholder1> 1, <reservedrangesplaceholder0> 1	s или A 0	`idwt`
Одноуровневая реконструкция	_{<reservedrangesplaceholder1> +1}, <reservedrangesplaceholder0> +1	_cAj	`idwt`
Полная реконструкция	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	s или A 0	`waverec`
Выборочная реконструкция	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	_Al, _Dm	`wrcoef`

Утилиты структуры разложения

Цель	Вход	Выход	Файл
Экстракция коэффициентов детализации	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	_cDk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	`detcoef`
Экстракция коэффициентов приближения	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	_cAk, 0 <reservedrangesplaceholder1> <reservedrangesplaceholder0>	`appcoef`
Перекомпозиция структуры разложения	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	[_cAk, _cDk..., <reservedrangesplaceholder2> 1] 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	`upwlev`

Чтобы проиллюстрировать режим командной строки для 1-D возможностей, смотрите 1-D Анализ с использованием командной строки..

2-D вейвлет

Основные 2-D объекты

	Объекты	Описание
Изображение в исходном разрешении	s	Оригинальное изображение
	<reservedrangesplaceholder0> 0	Приближение на уровне 0
	_Ak, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	Приближение на уровне k
	_Dk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	Детали на уровне k
Коэффициенты в масштабном разрешении	_cAk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	Приближения на уровне k
	_cDk, 1 ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>	Коэффициенты детализации на уровне k
	[_cAj, _cDj..., <reservedrangesplaceholder0> 1]	Вейвлет на уровне j

_Dk означает $[D_{k}^{(h)}, D_{k}^{(v)}, D_{k}^{(d)}]$ , горизонтальные, вертикальные и диагональные детали на уровне k.

То же самое относится и к _cDk, которая обозначает $[c D_{k}^{(h)}, c D_{k}^{(v)}, c D_{k}^{(d)}]$ .

Файлы 2-D те же, что и для 1-D случая, но с 2 в конце команды.

Для примера, idwt становится idwt2. Для получения дополнительной информации смотрите 1-D Возможности вейвлета.

Чтобы проиллюстрировать режим командной строки для возможностей 2-D, смотрите Вейвлет Изображения Анализ и сжатие..

Ссылки

[1] Mallat, S. G. «A Theory for Multirresolution Signal Decomposition: The Wavelet Representation», IEEE Transactions on Pattern Analysis and Machine Intelligence. Том 11, выпуск 7, июль 1989 года, стр. 674-693.

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