Банки ортогональных и биортогональных фильтров

Этот пример показывает, чтобы создать и использовать ортогональные и биортогональные банки фильтров с программным обеспечением Wavelet Toolbox. Классический критически выбранный двухканальный блок фильтров показан на следующем рисунке.

Пусть$\tilde{G}$ и$\tilde{H}$ обозначают lowpass и highpass фильтры анализа и$G$$H$ обозначают соответствующие lowpass и highpass фильтры синтеза. Двухканальный критически дискретизированный банк фильтров фильтрует входной сигнал с помощью lowpass и highpass фильтра. Выходы поддиапазона фильтров уменьшаются на две, чтобы сохранить общее количество выборок. Чтобы восстановить вход, увеличьте значение на два и затем интерполируйте результаты с помощью фильтров синтеза lowpass и highpass. Если фильтры удовлетворяют определенным свойствам, можно добиться идеальной реконструкции входов. Чтобы продемонстрировать это, фильтруйте сигнал ЭКГ, используя экстремальную фазу вейвлет Daubechies с двумя моментами исчезновения. Пример объясняет понятие моментов исчезновения в более позднем разделе.

load wecg;
plot(wecg);
title('ECG Signal')

Получите lowpass (масштабирование) и highpass (вейвлет) фильтры анализа и синтеза.

[gtilde,htilde,g,h] = wfilters('db2');

В данном примере установите режим заполнения для DWT в режим периодизации. Это не расширяет сигнал.

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

Получите дискретный вейвлет преобразование (DWT) уровня 1 сигнала ECG. Это эквивалентно ветви анализа (с понижающей дискретизацией) двухканальной группы фильтров на рисунке.

[lowpass,highpass] = dwt(wecg,gtilde,htilde);

Преобразуйте и интерполируйте поддиапазоны lowpass (коэффициенты масштабирования) и highpass (коэффициенты вейвлета) с помощью фильтров синтеза и демонстрируйте совершенную реконструкцию.

xrec = idwt(lowpass,highpass,g,h);
max(abs(wecg-xrec))
subplot(2,1,1)
plot(wecg); title('Original ECG Waveform')
subplot(2,1,2)
plot(xrec); title('Reconstructed ECG Waveform');
ans =

   1.3658e-12

Фильтры анализа и синтеза для вейвлета 'db2' являются всего лишь временными противоположными. Это можно увидеть, сравнив следующее.

scalingFilters = [flip(gtilde); g]
waveletFilters = [flip(htilde); h]
scalingFilters =

    0.4830    0.8365    0.2241   -0.1294
    0.4830    0.8365    0.2241   -0.1294


waveletFilters =

   -0.1294   -0.2241    0.8365   -0.4830
   -0.1294   -0.2241    0.8365   -0.4830

Это относится ко всем группам ортогональных вейвлет. Ортогональными семействами вейвлет, поддерживаемыми Wavelet Toolbox, являются 'dbN', 'fkN', 'symN' и 'coifN', где N является допустимым номером фильтра.

Вместо предоставления dwt с фильтрами в предыдущем примере, вы строка 'db2' вместо этого. Используя краткое имя семейства вейвлет и номер фильтра, вы не должны правильно задавать фильтры анализа и синтеза.

[lowpass,highpass] = dwt(wecg,'db2');
xrec = idwt(lowpass,highpass,'db2');

Число фильтров в экстремальной фазе Daubechies и наименее асимметричных фазовых вейвлетах ('db' и 'sym') относится к количеству исчезающих моментов. В основном, вейвлет с N моментами исчезновения удаляет полином порядка N-1 в вейвлет-коэффициентах. Чтобы проиллюстрировать это, создайте сигнал, который состоит из линейного тренда с аддитивным шумом.

n = (0:511)/512;
x = 2*n+0.2*randn(size(n));
plot(n,x)

Линейный тренд является полиномом степени 1. Поэтому вейвлет с двумя моментами исчезновения удаляет этот полином. Линейный тренд сохраняется в коэффициентах масштабирования, и коэффициенты вейвлета могут быть расценены как состоящие только из шума. Получите DWT уровня 1 сигнала с вейвлет 'db2' (два момента исчезновения) и постройте график коэффициентов.

[A,D] = dwt(x,'db2');
subplot(2,1,1)
plot(A); title('Scaling Coefficients');
subplot(2,1,2)
plot(D); title('Wavelet Coefficients');

Можно использовать dwt и idwt для реализации двухканальной ортогональной группы фильтров, но часто удобнее реализовать многоуровневую двухканальную группу фильтров с помощью wavedec. Многоуровневый DWT изменяется на выходе lowpass (масштабирование) фильтра. Другими словами, вход на второй уровень группы фильтров является выходом lowpass на уровне 1. Двухуровневая группа вейвлет проиллюстрирована на следующем рисунке.

На каждом последующем уровне количество коэффициентов масштабирования и вейвлет уменьшается на два, поэтому общее количество коэффициентов сохраняется. Получите DWT третьего уровня сигнала ЭКГ с помощью банка ортогональных фильтров 'sym4'.

[C,L] = wavedec(wecg,3,'sym4');

Количество коэффициентов по уровню содержится в векторе L. Первый элемент L равен 256, что представляет количество масштабирующих коэффициентов на уровне 3 (конечный уровень). Второй элемент L является количеством коэффициентов вейвлета на уровне 3. Последующие элементы дают количество коэффициентов вейвлета на более высоких уровнях, пока вы не достигнете конечного элемента L. Конечный элемент L равен количеству выборок в исходном сигнале. Масштабирование и вейвлет сохраняются в векторе C в том же порядке. Чтобы извлечь масштабирующие или вейвлет, используйте appcoef или detcoef. Извлеките все коэффициенты вейвлета в массиве ячеек и коэффициенты масштабирования последнего уровня.

wavcoefs = detcoef(C,L,'dcells');
a3 = appcoef(C,L,'sym4');

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

cfsmatrix = zeros(numel(wecg),4);
cfsmatrix(1:2:end,1) = wavcoefs{1};
cfsmatrix(1:4:end,2) = wavcoefs{2};
cfsmatrix(1:8:end,3) = wavcoefs{3};
cfsmatrix(1:8:end,4) = a3;
subplot(5,1,1)
plot(wecg); title('Original Signal');
axis tight;
for kk = 2:4
    subplot(5,1,kk)
    stem(cfsmatrix(:,kk-1),'marker','none','ShowBaseLine','off');
    ylabel(['D' num2str(kk-1)]);
    axis tight;
end
subplot(5,1,5);
stem(cfsmatrix(:,end),'marker','none','ShowBaseLine','off');
ylabel('A3'); xlabel('Sample');
axis tight;

Поскольку критически выбранная группа вейвлет понижает дискретизацию данных на каждом уровне, анализ должен прекратиться, когда у вас остается только один коэффициент. В случае сигнала ECG с 2048 выборками, это должно произойти, когда.$L = \log_2{2048}$

[C,L] = wavedec(wecg,log2(numel(wecg)),'sym4');
fprintf('The number of coefficients at the final level is %d. \n',L(1));
The number of coefficients at the final level is 1. 

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

ecgmodwt = modwt(wecg,'sym4',3);
ecgmra = modwtmra(ecgmodwt,'sym4');
subplot(5,1,1);
plot(wecg); title('Original Signal');

title('MODWT-Based Multiresolution Analysis');
for kk = 2:4
    subplot(5,1,kk)
    plot(ecgmra(kk-1,:));
    ylabel(['D' num2str(kk-1)]);
end
subplot(5,1,5);
plot(ecgmra(end,:));
ylabel('A3'); xlabel('Sample');

В группе биортогональных фильтров фильтры синтеза являются не просто реверсированными по времени версиями фильтров анализа. Семейство биортогональных сплайн вейвлет является примером таких фильтрующих блоков.

[LoD,HiD,LoR,HiR] = wfilters('bior3.5');

Если вы исследуете фильтры анализа (LoD, HiD) и фильтры синтеза (LoR, HiR), вы видите, что они сильно отличаются. Эти фильтрующие банки все еще обеспечивают идеальную реконструкцию.

[A,D] = dwt(wecg,LoD,HiD);
xrec = idwt(A,D,LoR,HiR);
max(abs(wecg-xrec))
ans =

   6.6613e-16

Биортогональные фильтры применяются, когда линейная фаза является требованием для вашего банка фильтров. Ортогональные фильтры не могут иметь линейную фазу за исключением вейвлет Haar. Если у вас есть Signal Processing Toolbox™, можно просмотреть фазовые отклики для ортогональной и биортогональной пары вейвлета фильтров.

[Lodb6,Hidb6] = wfilters('db6');
[PHIdb6,W] = phasez(Hidb6,1,512);
PHIbior35 = phasez(HiD,1,512);
figure;
subplot(2,1,1)
plot(W./(2*pi),PHIdb6); title('Phase Response for db6 Wavelet');
grid on;
xlabel('Cycles/Sample'); ylabel('Radians');
subplot(2,1,2)
plot(W./(2*pi),PHIbior35); title('Phase Response for bior3.5 Wavelet');
grid on;
xlabel('Cycles/Sample'); ylabel('Radians');

Верните dwtmode к исходной настройке.

dwtmode(origmodestatus,'nodisplay');