Загрузка 23-канальных данных ЭЭГ Espiga3 [3]. Данные дискретизируются на частоте 200 Гц.

load Espiga3

Вычислите 1-D DWPT данных с помощью sym3 вейвлет вниз до уровня 4. Получение узлов пакетного вейвлета терминала, вектора учета и центральных частот приблизительных полос пропускания.

[wpt,bk,~,f] = dwpt(Espiga3,'sym3','Level',4);

Продукция wpt является $${\mathrm{}}^{\mathrm{}}$$массивом ячеек 1 на 24. Каждый элемент wpt является матрицей. Выберите любой терминальный узел и подтвердите размер матрицы 23-by-M, где M - последний элемент вектора бухгалтерского учета. bk.

nd = 13;
size(wpt{nd})

ans = 1×2

    23    66

bk(end)

ans = 66

Извлеките коэффициенты конечного уровня пятого канала.

p5 = cell2mat(cellfun(@(x) x(5,:).',wpt,'UniformOutput',false));
size(p5)

ans = 1×2

    66    16

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

plot(200*f,'x')
title('Center Frequencies')
ylabel('Hz')

В этом примере показано, как взять выражение пары биорогональных фильтров и построить фильтры нижних и верхних частот для создания пары совершенной реконструкции (PR) в Toolbox™ вейвлета.

Фильтр LeGall 5/3 - это вейвлет, используемый в JPEG2000 для сжатия изображения без потерь. Фильтры нижних частот (масштабирование) для вейвлета LeGall 5/3 имеют пять и три ненулевых коэффициента соответственно. Выражения для этих двух фильтров:

Создайте эти фильтры.

H0 = 1/8*[-1 2 6 2 -1];
H1 = 1/2*[1 2 1];

Многие из дискретных вейвлет-преобразований и вейвлет-пакетов в инструментарии вейвлета полагаются на то, что фильтры имеют четную длину и одинаковую длину для того, чтобы получить идеальный набор восстановительных фильтров, связанный с этими преобразованиями. Эти преобразования также требуют конкретной нормализации коэффициентов в фильтрах для алгоритмов, чтобы создать банк PR-фильтров. Используйте biorfilt функция прототипа нижних частот для создания набора импульсных фильтров PR.

[LoD,HiD,LoR,HiR] = biorfilt(H0,H1);

Сумма фильтров низкочастотного анализа и синтеза теперь равна $$\sqrt{2}$$.

sum(LoD)

ans = 1.4142

sum(LoR)

ans = 1.4142

Вейвлет-фильтры суммируют, при необходимости, до нуля. L2-norms низкочастотного анализа и фильтров синтеза верхних частот равны. То же самое относится к фильтрам синтеза нижних частот и анализа верхних частот.

Теперь вы можете использовать эти фильтры в дискретных преобразованиях вейвлет и вейвлет-пакетов и получить банк фильтров эйвлет-пакетов PR. Чтобы продемонстрировать это, загрузите и постройте график сигнала ЭКГ.

load wecg
plot(wecg)
axis tight
grid on

Получение дискретного вейвлет-пакетного преобразования ЭКГ-сигнала с использованием набора фильтров LeGall 5/3.

[wpt,L] = dwpt(wecg,LoD,HiD);

Теперь используйте фильтры реконструкции (синтеза) для восстановления сигнала и демонстрации идеальной реконструкции.

xrec = idwpt(wpt,L,LoR,HiR);
plot([wecg xrec])
axis tight, grid on;

norm(wecg-xrec,'Inf')

ans = 3.3307e-15

Этот банк фильтров можно также использовать в 1-D и 2-D дискретных вейвлет-преобразованиях. Чтение и печать изображения.

im = imread('woodsculp256.jpg');
image(im); axis off;

Получение 2-D вейвлет-преобразования с помощью фильтров анализа LeGall 5/3.

[C,S] = wavedec2(im,3,LoD,HiD);

Реконструируйте изображение с помощью фильтров синтеза.

imrec = waverec2(C,S,LoR,HiR);
image(uint8(imrec)); axis off;

Фильтр LeGall 5/3 эквивалентен встроенному 'bior2.2' вейвлет в Wavelet Toolbox. Используйте 'bior2.2' и сравните с фильтрами LeGall 5/3.

[LD,HD,LR,HR] = wfilters('bior2.2');
subplot(2,2,1)
hl = stem([LD' LoD']);
hl(1).MarkerFaceColor = [0 0 1];
hl(1).Marker = 'o';
hl(2).MarkerFaceColor = [1 0 0];
hl(2).Marker = '^';
grid on
title('Lowpass Analysis')
subplot(2,2,2)
hl = stem([HD' HiD']);
hl(1).MarkerFaceColor = [0 0 1];
hl(1).Marker = 'o';
hl(2).MarkerFaceColor = [1 0 0];
hl(2).Marker = '^';
grid on
title('Highpass Analysis')
subplot(2,2,3)
hl = stem([LR' LoR']);
hl(1).MarkerFaceColor = [0 0 1];
hl(1).Marker = 'o';
hl(2).MarkerFaceColor = [1 0 0];
hl(2).Marker = '^';
grid on
title('Lowpass Synthesis')
subplot(2,2,4)
hl = stem([HR' HiR']);
hl(1).MarkerFaceColor = [0 0 1];
hl(1).Marker = 'o';
hl(2).MarkerFaceColor = [1 0 0];
hl(2).Marker = '^';
grid on
title('Highpass Synthesis')