Создайте банк фильтров 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 с помощью пользовательских биортогональных вейвлет.

Две пары фильтров анализа (разложения) и синтеза (реконструкции) связаны с биортогональным вейвлет. Каждая пара состоит из lowpass и highpass фильтра. Укажите фильтры анализа и синтеза для почти ортогональных биортогональных вейвлеты на основе схемы лаплакской пирамиды Берта и Адельсона (таблица 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 являются lowpass и highpass фильтрами разложения, соответственно. Hr и Gr являются lowpass и highpass фильтрами реконструкции, соответственно.

Создайте анализы и синтез блоков фильтров 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

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

[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. Постройте график суммы квадратов величин частотных характеристик. Поскольку фильтры масштабирования (lowpass) и вейвлет (highpass) не образуют ортогональную квадратурную пару зеркальных фильтров, сумма не равняется 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')