Exponenta Banner

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

Скрипт Open Live Script

В этом примере показано, как добавить ортогональный квадратурный фильтр зеркала (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, чтобы проверить, что сумма равняется 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. Версия ДПФ этого равенства: |G2mkmodN|2+|H2mkmodN|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.

Документация Wavelet Toolbox
Поддержка