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

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

Этот пример добавляет Бейлкина (18) пара фильтра QMF к тулбоксу и показывает, как впоследствии использовать фильтр в дискретном анализе вейвлета. Пример затем демонстрирует, как проверить необходимые и достаточные условия для пары 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);

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

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')

Наконец, проверьте, что новый фильтр удовлетворяет условия для ортогональной пары 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');

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

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

Итоговый критерий утверждает, что сумма значений в квадрате преобразований Фурье масштабирования и фильтров вейвлета на каждой частоте равна 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;

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

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

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

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

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

Чтобы быть распознанным 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)

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

[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')

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

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

synthesisLowerBound = 0.9800

synthesisUpperBound = 1.0509

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

Удалите фильтр Берта-Адельсона из Тулбокса.

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

Ссылки

[1] Daubechies, я. Десять лекций по вейвлетам. CBMS-NSF региональный ряд конференции в прикладной математике. Филадельфия, PA: общество промышленной и прикладной математики, 1992.

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