Adding a Quadrature Mirror Filter
This example shows how to add an orthogonal quadrature mirror filter (QMF) pair to the Wavelet Toolbox™. While Wavelet Toolbox™ already contains many of the most widely used orthogonal QMF families, including the Daubechies' extremal-phase, the Daubechies' least-asymmetric phase, the coiflet, and the Fejer-Korovkin filters, you can easily add your own QMF filters and use the filter in any of the discrete wavelet or wavelet packet algorithms. This example adds the Beylkin(18) QMF filter to the toolbox and shows how to subsequently use the filter in discrete wavelet analysis. Finally, the example demonstrates how to verify the necessary and sufficient conditions for the QMF pair to constitute a scaling and wavelet filter.
First, you must have some way of obtaining the coefficients. In this case, here are the coefficients for the lowpass (scaling) Beylkin(18) filter. You only need a valid scaling filter, wfilters creates the corresponding wavelet filter for you.
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];
Save the Beylkin(18) filter and add the new filter to the toolbox.
save beyl beyl
Use wavemngr to add the wavelet filter to the toolbox. Define the wavelet family name and the short name used to access the filter. Define the wavelet type to be 1. Type 1 wavelets are orthogonal wavelets in the toolbox. Because you are adding only one wavelet in this family, define the NUMS variable input to wavemngr to be an empty string.
familyName = 'beylkin'; familyShortName = 'beyl'; familyWaveType = 1; familyNums = ''; fileWaveName = 'beyl.mat';
Add the wavelet using wavemngr.
wavemngr('add',familyName,familyShortName,familyWaveType, ... familyNums,fileWaveName)
Verify that the wavelet has been added to the toolbox.
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 ' '==================================='
You can now use the wavelet to analyze signals or images. For example, load an ECG signal and obtain the MODWT of the signal down to level four using the Beylkin(18) filter.
load wecg wtecg = modwt(wecg,'beyl',4);
Load a box image, obtain the 2-D DWT using the Beylkin(18) filter. Show the level-one diagonal detail coefficients.
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')
Finally, verify that the new filter satisfies the conditions for an orthogonal QMF pair. Obtain the scaling (lowpass) and wavelet (highpass) filters.
[Lo,Hi] = wfilters('beyl');
Sum the lowpass filter coefficients to verify that the sum equals . Sum the wavelet filter coefficients and verify that the sum is 0.
ans = 1.4142 ans = -1.9873e-16
Verify that the autocorrelation of the scaling and wavelet filters at all even nonzero lags is 0. You have to have the Signal Processing Toolbox™ to use 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');
Note that the autocorrelation values in both plots is zero for nonzero even lags. Show that the cross-correlation of the scaling and wavelet filter is zero at all even lags.
[xcr,lags] = xcorr(Lo,Hi,10); figure stem(lags,xcr,'markerfacecolor',[0 0 1]); title('Cross-correlation of QMF filters')
The final criterion states the sum of squared magnitudes of the Fourier transforms of scaling and wavelet filters at each frequency is equal to 2. In other words, let be the Fourier transform of the scaling filter and be the Fourier transform of the wavelet filter. The following holds for all : . The DFT version of this equality is: for any . Check this for the Beylkin(18) filter with .
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 = 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
All the values are equal to 2 as expected. To understand why these filters are called quadrature mirror filters, visualize the squared-magnitude frequency responses of the scaling and wavelet filters.
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;
Note the magnitude responses are symmetric, or mirror images, of each other around the quadrature frequency of 1/4.
The following code removes the Beylkin(18) wavelet filter.