Этот пример показывает построению и использованию ортогональные и биоортогональные наборы фильтров с программным обеспечением Wavelet Toolbox. Классика критически произвела двухканальный набор фильтров, показан в следующем рисунке.
Позвольте и обозначьте lowpass и highpass аналитические фильтры и и обозначьте соответствующий lowpass и highpass фильтры синтеза. Двухканальный критически произведенный набор фильтров фильтрует входной сигнал с помощью фильтра highpass и lowpass. Поддиапазон выходные параметры фильтров прорежен два, чтобы сохранить общее количество выборок. Чтобы восстановить вход, сверхдискретизируйте два и затем интерполируйте результаты с помощью lowpass и highpass фильтров синтеза. Если фильтры удовлетворяют определенным свойствам, можно достигнуть совершенной реконструкции входа. Чтобы продемонстрировать это, отфильтруйте сигнал ECG использование экстремального вейвлета фазы Добечиса с двумя исчезающими моментами. Пример объясняет понятие исчезающих моментов в более позднем разделе.
load wecg; plot(wecg); title('ECG Signal')
Получите lowpass (масштабирование) и highpass (вейвлет) фильтры синтеза и анализ.
[gtilde,htilde,g,h] = wfilters('db2');
В данном примере установите дополнительный режим для DWT к periodization. Это не расширяет сигнал.
origmodestatus = dwtmode('status','nodisplay'); dwtmode('per','nodisplay');
Получите уровень один дискретный вейвлет преобразовывает (DWT) сигнала 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');
Номер фильтра в экстремальной фазе Добечиса и наименьшем количестве асимметричных вейвлетов фазы ('дб' и 'sym') относится к номеру исчезающих моментов. В основном вейвлет с N исчезающие моменты удаляет полином порядка n-1 в коэффициентах вейвлета. Чтобы проиллюстрировать это, создайте сигнал, который состоит из линейного тренда с аддитивным шумом.
n = (0:511)/512; x = 2*n+0.2*randn(size(n)); plot(n,x)
Линейный тренд является полиномом степени 1. Поэтому вейвлет с двумя исчезающими моментами удаляет этот полином. Линейный тренд сохраняется в масштабных коэффициентах, и коэффициенты вейвлета могут рассматриваться как состоящий только из шума. Получите уровень один DWT сигнала с '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 сигнала ECG использование '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 с 2 048 выборками это должно произойти когда.
[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');
Если вы исследуете аналитические фильтры (ЛОД, HiD) и фильтры синтеза (LoR, HiR), вы видите, что они очень отличаются. Эти наборы фильтров все еще обеспечивают совершенную реконструкцию.
[A,D] = dwt(wecg,LoD,HiD); xrec = idwt(A,D,LoR,HiR); max(abs(wecg-xrec))
ans = 6.6613e-16
Биоортогональные фильтры полезны, когда линейная фаза является требованием для вашего набора фильтров. Ортогональные фильтры не могут иметь линейной фазы за исключением фильтра вейвлета Хаара. Если у вас есть 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');