Создайте банк фильтров DWT с использованием значений по умолчанию.

fb = dwtfilterbank

fb = 
  dwtfilterbank with properties:

                Wavelet: 'sym4'
           SignalLength: 1024
                  Level: 7
      SamplingFrequency: 1
             FilterType: 'Analysis'
    CustomWaveletFilter: []
    CustomScalingFilter: []

Постройте график частотно-амплитудных характеристик вейвлетов и масштабной функции грубого масштаба. Откройте график в отдельном окне рисунка. Легенда графика в окне является интерактивной. Чтобы скрыть определенную частотную характеристику, щелкните ее имя.

freqz(fb)

Получают и строят график центрированных по времени вейвлетов, соответствующих полосовым фильтрам вейвлетов.

[psi,t] = wavelets(fb);
plot(t,psi')
grid on
title('Time-Centered Wavelets')
xlabel('Time')
ylabel('Magnitude')

В этом примере показано, как создать банк фильтров DWT с помощью пользовательских биоргональных вейвлет-фильтров.

Две пары фильтров анализа (разложения) и синтеза (реконструкции) связаны с биоргональным вейвлетом. Каждая пара состоит из фильтра нижних и верхних частот. Укажите фильтры анализа и синтеза для почти ортогональных биоргональных вейвлетов на основе схемы лапласианской пирамиды Берта и Адельсона (таблица 8.4 на стр. 283 в [1]). Поскольку набор инструментов требует, чтобы все фильтры, связанные с биоргональным вейвлетом или ортогональным вейвлетом, имели одинаковую четную длину, фильтры дополняются 0.

Hd = [0 -1 5 12 5 -1 0 0]/20*sqrt(2);
Gd = [0 3 -15 -73 170 -73 -15 3]/280*sqrt(2);
Hr = [0 -3 -15 73 170 73 -15 -3]/280*sqrt(2);
Gr = [0 -1 -5 12 -5 -1 0 0]/20*sqrt(2);

Hd и Gd - фильтры разложения нижних и верхних частот соответственно. Hr и Gr являются фильтрами восстановления нижних и верхних частот соответственно.

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

fbAna = dwtfilterbank('Wavelet','Custom',...
    'CustomScalingFilter',[Hd' Hr'],'CustomWaveletFilter',[Gd' Gr']);
fbSyn = dwtfilterbank('Wavelet','Custom',...
    'CustomScalingFilter',[Hd' Hr'],'CustomWaveletFilter',[Gd' Gr'],...
    'FilterType','Synthesis');
fprintf('fbAna: isOrthogonal = %d\tisBiorthogonal = %d\n',...
    isOrthogonal(fbAna),isBiorthogonal(fbAna));

fbAna: isOrthogonal = 0	isBiorthogonal = 1

fprintf('fbSyn: isOrthogonal = %d\tisBiorthogonal = %d\n',...
    isOrthogonal(fbSyn),isBiorthogonal(fbSyn ));

fbSyn: isOrthogonal = 0	isBiorthogonal = 1

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

[fbAna_phi,t] = scalingfunctions(fbAna);
[fbAna_psi,~] = wavelets(fbAna);
[fbSyn_phi,~] = scalingfunctions(fbSyn);
[fbSyn_psi,~] = wavelets(fbSyn);
subplot(2,2,1)
plot(t,fbAna_phi(end,:))
grid on
title('Analysis - Scaling')
subplot(2,2,2)
plot(t,fbAna_psi(end,:))
grid on
title('Analysis - Wavelet')
subplot(2,2,3)
plot(t,fbSyn_phi(end,:))
grid on
title('Synthesis - Scaling')
subplot(2,2,4)
plot(t,fbSyn_psi(end,:))
grid on
title('Synthesis - Wavelet')

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

[a1,a2] = framebounds(fbAna)

a1 = 0.9505

a2 = 1.0211

[b1,b2] = framebounds(fbSyn)

b1 = 0.9800

b2 = 1.0528

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

[psidft,f,phidft] = freqz(fbAna);
flen = length(f);
figure
plot(f(flen/2+1:end),abs(phidft(1,flen/2+1:end)))
hold on
plot(f(flen/2+1:end),abs(psidft(1,flen/2+1:end)))
grid on
legend('Scaling','Wavelet')
title('Frequency Responses')
xlabel('Normalized Frequency')
ylabel('Magnitude')

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

figure
plot(f(flen/2+1:end),abs(phidft(1,flen/2+1:end)).^2 + abs(psidft(1,flen/2+1:end)).^2)
grid on
title('Sum of Squared Frequency Responses')
xlabel('Normalized Frequency')
ylabel('Sum of Magnitudes')