Добавление квадратурных зеркальных и биортогональных фильтров вейвлет

Открыть Live Script

В этом примере показано, как добавить пару ортогонального квадратурного зеркального фильтра (QMF) и биортогональный вейвлет в четверик к Wavelet Toolbox™. В то время как Wavelet Toolbox™ уже содержит многие наиболее широко используемые ортогональные и биортогональные семейства вейвлет, включая экстремальную фазу Daubechies, наименее асимметричную фазу Daubechies, койфлет, фильтры Фейера-Коровкина и биортогональные сплайновые вейвлеты, вы можете легко добавить свои собственные

Этот пример добавляет пару фильтров QMF Бейлкина (18) в тулбокс и показывает, как впоследствии использовать фильтр в дискретном вейвлет. Затем пример демонстрирует, как проверить необходимые и достаточные условия для того, чтобы пара QMF составляла масштабный и вейвлет. После добавления пары QMF в пример добавляется почти ортогональный биортогональный вейвлет, основанный на лаплакийской пирамидной схеме Берта и Адельсона (таблица 8.4 на стр. 283 в [1]).

Добавление QMF

Во-первых, у вас должен быть какой-то способ получить коэффициенты. В этом случае вот коэффициенты для lowpass (масштабирование) фильтра Бейлкина (18). Вам нужен только допустимый масштабирующий фильтр, wfilters создает соответствующий вейвлет для вас.

beyl = [9.93057653743539270E-02
    4.24215360812961410E-01
    6.99825214056600590E-01 
    4.49718251149468670E-01
    -1.10927598348234300E-01
    -2.64497231446384820E-01
    2.69003088036903200E-02
    1.55538731877093800E-01
    -1.75207462665296490E-02
    -8.85436306229248350E-02
    1.96798660443221200E-02
    4.29163872741922730E-02
    -1.74604086960288290E-02
    -1.43658079688526110E-02
    1.00404118446319900E-02
    1.48423478247234610E-03
    -2.73603162625860610E-03
    6.40485328521245350E-04];

Сохраните фильтр Бейлкина (18) и добавьте новый фильтр в тулбокс.

save beyl beyl

Использование wavemngr для добавления вейвлета фильтра к тулбоксу. Задайте имя семейства вейвлет и краткое имя, используемое для доступа к фильтру. Определите тип вейвлета как 1. Вейвлеты типа 1 являются ортогональными вейвлетами в тулбоксе. Поскольку вы добавляете только один вейвлет в это семейство, задайте вход переменной NUMS в wavemngr быть пустой строкой.

familyName      = 'beylkin';
familyShortName = 'beyl';
familyWaveType  = 1;
familyNums      = '';
fileWaveName    = 'beyl.mat';

Добавьте вейвлет с помощью wavemngr.

wavemngr('add',familyName,familyShortName,familyWaveType, ...
    familyNums,fileWaveName)

Убедитесь, что вейвлет добавлен в тулбокс.

wavemngr('read')

ans = 19x35 char array
    '==================================='
    'Haar              ->->haar           '
    'Daubechies        ->->db             '
    'Symlets           ->->sym            '
    'Coiflets          ->->coif           '
    'BiorSplines       ->->bior           '
    'ReverseBior       ->->rbio           '
    'Meyer             ->->meyr           '
    'DMeyer            ->->dmey           '
    'Gaussian          ->->gaus           '
    'Mexican_hat       ->->mexh           '
    'Morlet            ->->morl           '
    'Complex Gaussian  ->->cgau           '
    'Shannon           ->->shan           '
    'Frequency B-Spline->->fbsp           '
    'Complex Morlet    ->->cmor           '
    'Fejer-Korovkin    ->->fk             '
    'beylkin           ->->beyl           '
    '==================================='

Теперь можно использовать вейвлет для анализа сигналов или изображений. Для примера загрузите сигнал ECG и получите MODWT сигнала до четвертого уровня с помощью фильтра Бейлкина (18).

load wecg
wtecg = modwt(wecg,'beyl',4);

Загрузите прямоугольное изображение, получите 2-D DWT с помощью фильтра Бейлкина (18). Отобразите диагональные коэффициенты детализации уровня 1.

load xbox
[C,S] = wavedec2(xbox,1,'beyl');
[H,V,D] = detcoef2('all',C,S,1);
subplot(2,1,1)
imagesc(xbox)
axis off
title('Original Image')
subplot(2,1,2)
imagesc(D)
axis off
title('Level-One Diagonal Coefficients')

Figure contains 2 axes. Axes 1 with title Original Image contains an object of type image. Axes 2 with title Level-One Diagonal Coefficients contains an object of type image.

Наконец, проверьте, что новый фильтр удовлетворяет условиям для ортогональной пары QMF. Получите фильтры масштабирования (lowpass) и вейвлет (highpass).

[Lo,Hi] = wfilters('beyl');

Суммируйте коэффициенты lowpass, чтобы убедиться, что сумма равна $\sqrt{2}$ . Суммируйте коэффициенты вейвлета фильтра и проверьте, что сумма равна 0.

sum(Lo)

ans = 1.4142

sum(Hi)

ans = -1.9873e-16

Проверьте, что автокорреляция масштабирующего и вейвлет при всех ненулевых лагах равна 0. У вас должны быть Signal Processing Toolbox™ для использования xcorr.

[Clow,lags] = xcorr(Lo,Lo,10);
Chigh = xcorr(Hi,Hi,10);
subplot(2,1,1)
stem(lags,Clow,'markerfacecolor',[0 0 1])
grid on;
title('Autocorrelation of Scaling Filter');
subplot(2,1,2)
stem(lags,Chigh,'markerfacecolor',[0 0 1])
grid on;
title('Autocorrelation of Wavelet Filter');

Figure contains 2 axes. Axes 1 with title Autocorrelation of Scaling Filter contains an object of type stem. Axes 2 with title Autocorrelation of Wavelet Filter contains an object of type stem.

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

[xcr,lags] = xcorr(Lo,Hi,10);
figure
stem(lags,xcr,'markerfacecolor',[0 0 1]);
title('Cross-correlation of QMF filters')

Figure contains an axes. The axes with title Cross-correlation of QMF filters contains an object of type stem.

Конечный критерий утверждает, что сумма квадратов величин преобразований Фурье масштабирования и вейвлет фильтров на каждой частоте равна 2. Другими словами, позвольте $G (f)$ быть преобразованием Фурье масштабирующего фильтра и $H (f)$ быть преобразованием Фурье вейвлет. Следующее удерживает для всех $f$ : $| H (f) |^{2} + | G (f) |^{2} = 2$ . Версия этого равенства для ДПФ: $| G_{2^{m} k m o d N} |^{2} + | H_{2^{m} k m o d N} |^{2} = 2$ для любого $m$ . Проверьте фильтр Бейлкина (18) с $m = 0$ .

N = numel(Lo);
LoDFT = fft(Lo);
HiDFT = fft(Hi);
k = 0:N-1;
m = 0;
sumDFTmags = abs(LoDFT(1+mod(2^m*k,N))).^2+abs(HiDFT(1+mod(2^m*k,N))).^2

sumDFTmags = 18×1

    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
      ⋮

Все значения равны 2, как и ожидалось. Чтобы понять, почему эти фильтры называются квадратурными зеркальными фильтрами, визуализируйте квадратурные частотные характеристики масштабирующего и вейвлет.

nfft = 512;
F = 0:1/nfft:1/2;
LoDFT = fft(Lo,nfft);
HiDFT = fft(Hi,nfft);
figure
plot(F,abs(LoDFT(1:nfft/2+1)).^2,'DisplayName','Scaling Filter');
hold on
plot(F,abs(HiDFT(1:nfft/2+1)).^2,'r','DisplayName','Wavelet Filter');
xlabel('Frequency'); ylabel('Squared Magnitude')
title('Beylkin(18) QMF Filter Pair')
grid on
plot([1/4 1/4], [0 2],'k','HandleVisibility','off');
legend show;

Figure contains an axes. The axes with title Beylkin(18) QMF Filter Pair contains 2 objects of type line. These objects represent Scaling Filter, Wavelet Filter.

Обратите внимание, что величины отклики являются симметричными или зеркальными изображениями друг друга вокруг квадратурной частоты 1/4.

Следующий код удаляет вейвлет Бейлкина (18).

wavemngr('del',familyShortName);
delete('beyl.mat')

Добавление биортогонального вейвлета

Добавление биортогонального вейвлета к тулбоксу аналогично добавлению QMF. Вы предоставляете допустимую пару lowpass (масштабирующих) фильтров, используемых в анализе и синтезе. The wfilters функция сгенерирует высокочастотные фильтры.

Подлежит распознаванию wfilters, фильтр масштабирования анализа должен быть назначен переменной Df, и синтезирующий масштабирующий фильтр должен быть назначен переменной Rf. Биортогональные масштабирующие фильтры не должны быть даже равной длины. Созданные выходы пары биортогональных фильтров будут иметь даже равные длины. Вот пары масштабирующих функций почти ортогонального биортогонального вейвлет, основанного на лаплакийской пирамидной схеме Берта и Адельсона.

Df = [-1 5 12 5 -1]/20*sqrt(2);
Rf = [-3 -15 73 170 73 -15 -3]/280*sqrt(2);

Сохраните фильтры в .mat файл.

save burt Df Rf

Использование wavemngr для добавления биортогональных фильтров вейвлета к тулбоксу. Задайте имя семейства вейвлет и краткое имя, используемое для доступа к фильтру. Поскольку вейвлеты являются биортогональными, установите тип вейвлета равным 2. Поскольку вы добавляете только один вейвлет в это семейство, задайте вход переменной NUMS в wavemngr быть пустой строкой.

familyName      = 'burtAdelson';
familyShortName = 'burt';
familyWaveType  = 2;
familyNums      = '';
fileWaveName    = 'burt.mat';
wavemngr('add',familyName,familyShortName,familyWaveType,...
    familyNums,fileWaveName)

Проверьте, что биортогональный вейвлет добавлен в тулбокс.

wavemngr('read')

ans = 19x35 char array
    '==================================='
    'Haar              ->->haar           '
    'Daubechies        ->->db             '
    'Symlets           ->->sym            '
    'Coiflets          ->->coif           '
    'BiorSplines       ->->bior           '
    'ReverseBior       ->->rbio           '
    'Meyer             ->->meyr           '
    'DMeyer            ->->dmey           '
    'Gaussian          ->->gaus           '
    'Mexican_hat       ->->mexh           '
    'Morlet            ->->morl           '
    'Complex Gaussian  ->->cgau           '
    'Shannon           ->->shan           '
    'Frequency B-Spline->->fbsp           '
    'Complex Morlet    ->->cmor           '
    'Fejer-Korovkin    ->->fk             '
    'burtAdelson       ->->burt           '
    '==================================='

Теперь можно использовать вейвлет в тулбоксе. Создайте банк фильтров DWT для анализа с помощью burt вейвлет. Подтвердите, что группа фильтров DWT является биортогональной. Постройте график величины частотных характеристик вейвлета полосно-пропускающих фильтров и функции масштабирования с наибольшим разрешением.

fb = dwtfilterbank('Wavelet','burt');
isBiorthogonal(fb)

ans = logical
   1

freqz(fb)

Figure contains an axes. The axes with title DWT Filter Bank burt contains 8 objects of type line. These objects represent D 1, D 2, D 3, D 4, D 5, D 6, D 7, A 7.

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

[fb_phi,t] = scalingfunctions(fb);
[fb_psi,~] = wavelets(fb);
subplot(2,1,1)
plot(t,fb_phi(end,:))
axis tight
grid on
title('Analysis - Scaling')
subplot(2,1,2)
plot(t,fb_psi(end,:))
axis tight
grid on
title('Analysis - Wavelet')

Figure contains 2 axes. Axes 1 with title Analysis - Scaling contains an object of type line. Axes 2 with title Analysis - Wavelet contains an object of type line.

Создайте синтез набора фильтров DWT с помощью burt вейвлет. Вычислите границы фрейма.

fb2 = dwtfilterbank('Wavelet','burt','FilterType','Synthesis','Level',4);
[synthesisLowerBound,synthesisUpperBound] = framebounds(fb2)

synthesisLowerBound = 0.9800

synthesisUpperBound = 1.0509

Получите lowpass и highpass фильтры анализа и синтеза, сопоставленные с burt. Обратите внимание, что выходные фильтры имеют одинаковую четную длину. Подтвердите сумму коэффициентов lowpass к sqrt(2) и коэффициенты высокочастотного фильтра равны 0.

[LoD,HiD,LoR,HiR] = wfilters('burt');
[LoD' HiD' LoR' HiR']

ans = 8×4

         0    0.0152   -0.0152         0
         0   -0.0758   -0.0758         0
   -0.0707   -0.3687    0.3687   -0.0707
    0.3536    0.8586    0.8586   -0.3536
    0.8485   -0.3687    0.3687    0.8485
    0.3536   -0.0758   -0.0758   -0.3536
   -0.0707    0.0152   -0.0152   -0.0707
         0         0         0         0

sum([(LoD'/sqrt(2)) HiD' (LoR'/sqrt(2)) HiR'])

ans = 1×4

    1.0000   -0.0000    1.0000         0

Удалите фильтр Burt-Adelson из тулбокса.

wavemngr('del',familyShortName);
delete('burt.mat')

Ссылки

[1] Daubechies, I. Ten Lectures on Wavelets. Серия региональных конференций CBMS-NSF по прикладной математике. Филадельфия, Пенсильвания: Общество промышленной и прикладной математики, 1992.

Документация

Добавление квадратурных зеркальных и биортогональных фильтров вейвлет

Добавление QMF

Добавление биортогонального вейвлета

Ссылки

Документация по Wavelet Toolbox

Поддержка