MATLAB Examples

Orthogonal and Biorthogonal Filter Banks

This example shows to construct and use orthogonal and biorthogonal filter banks with the Wavelet Toolbox software. The classic critically sampled two-channel filter bank is shown in the following figure.

Let $\tilde{G}$ and $\tilde{H}$ denote the lowpass and highpass analysis filters and $G$ and $H$ denote the corresponding lowpass and highpass synthesis filters. A two-channel critically sampled filter bank filters the input signal using a lowpass and highpass filter. The subband outputs of the filters are downsampled by two to preserve the overall number of samples. To reconstruct the input, upsample by two and then interpolate the results using the lowpass and highpass synthesis filters. If the filters satisfy certain properties, you can achieve perfect reconstruction of the input. To demonstrate this, filter an ECG signal using Daubechies's extremal phase wavelet with two vanishing moments. The example explains the notion of vanishing moments in a later section.

load wecg;
plot(wecg);
title('ECG Signal')

Obtain the lowpass (scaling) and highpass (wavelet) analysis and synthesis filters.

[gtilde,htilde,g,h] = wfilters('db2');

For this example, set the padding mode for the DWT to periodization. This does not extend the signal.

origmodestatus = dwtmode('status','nodisplay');
dwtmode('per','nodisplay');

Obtain the level-one discrete wavelet transform (DWT) of the ECG signal. This is equivalent to the analysis branch (with downsampling) of the two-channel filter bank in the figure.

[lowpass,highpass] = dwt(wecg,gtilde,htilde);

Upsample and interpolate the lowpass (scaling coefficients) and highpass (wavelet coefficients) subbands with the synthesis filters and demonstrate perfect reconstruction.

xrec = idwt(lowpass,highpass,g,h);
max(abs(wecg-xrec))
subplot(2,1,1)
plot(wecg); title('Original ECG Waveform')
subplot(2,1,2)
plot(xrec); title('Reconstructed ECG Waveform');
ans =

   1.3658e-12

The analysis and synthesis filters for the 'db2' wavelet are just time reverses of each other. You can see this by comparing the following.

scalingFilters = [flip(gtilde); g]
waveletFilters = [flip(htilde); h]
scalingFilters =

    0.4830    0.8365    0.2241   -0.1294
    0.4830    0.8365    0.2241   -0.1294


waveletFilters =

   -0.1294   -0.2241    0.8365   -0.4830
   -0.1294   -0.2241    0.8365   -0.4830

This is the case with all orthogonal wavelet filterbanks. The orthogonal wavelet families supported by the Wavelet Toolbox are 'dbN', 'fkN', 'symN', and 'coifN' where N is a valid filter number.

Instead of providing dwt with the filters in the previous example, you the string 'db2' instead. Using the wavelet family short name and filter number, you do not have to correctly specify the analysis and synthesis filters.

[lowpass,highpass] = dwt(wecg,'db2');
xrec = idwt(lowpass,highpass,'db2');

The filter number in the Daubechies's extremal phase and least asymmetric phase wavelets ('db' and 'sym') refers to the number of vanishing moments. Basically, a wavelet with N vanishing moments removes a polynomial of order N-1 in the wavelet coefficients. To illustrate this, construct a signal which consists of a linear trend with additive noise.

n = (0:511)/512;
x = 2*n+0.2*randn(size(n));
plot(n,x)

A linear trend is a polynomial of degree 1. Therefore, a wavelet with two vanishing moments removes this polynomial. The linear trend is preserved in the scaling coefficients and the wavelet coefficients can be regarded as consisting of only noise. Obtain the level-one DWT of the signal with the 'db2' wavelet (two vanishing moments) and plot the coefficients.

[A,D] = dwt(x,'db2');
subplot(2,1,1)
plot(A); title('Scaling Coefficients');
subplot(2,1,2)
plot(D); title('Wavelet Coefficients');

You can use dwt and idwt to implement a two-channel orthogonal filter bank, but it is often more convenient to implement a multi-level two-channel filter bank using wavedec. The multi-level DWT iterates on the output of the lowpass (scaling) filter. In other words, the input to the second level of the filter bank is the output of the lowpass filter at level 1. A two-level wavelet filter bank is illustrated in the following figure.

At each successive level, the number of scaling and wavelet coefficients is downsampled by two so the total number of coefficients are preserved. Obtain the level three DWT of the ECG signal using the 'sym4' orthogonal filter bank.

[C,L] = wavedec(wecg,3,'sym4');

The number of coefficients by level is contained in the vector, L. The first element of L is equal to 256, which represents the number of scaling coefficients at level 3 (the final level). The second element of L is the number of wavelet coefficients at level 3. Subsequent elements give the number of wavelet coefficients at higher levels until you reach the final element of L. The final element of L is equal to the number of samples in the original signal. The scaling and wavelet coefficients are stored in the vector C in the same order. To extract the scaling or wavelet coefficients, use appcoef or detcoef. Extract all the wavelet coefficients in a cell array and final-level scaling coefficients.

wavcoefs = detcoef(C,L,'dcells');
a3 = appcoef(C,L,'sym4');

You can plot the wavelet and scaling coefficients at their approximate positions.

cfsmatrix = zeros(numel(wecg),4);
cfsmatrix(1:2:end,1) = wavcoefs{1};
cfsmatrix(1:4:end,2) = wavcoefs{2};
cfsmatrix(1:8:end,3) = wavcoefs{3};
cfsmatrix(1:8:end,4) = a3;
subplot(5,1,1)
plot(wecg); title('Original Signal');
axis tight;
for kk = 2:4
    subplot(5,1,kk)
    stem(cfsmatrix(:,kk-1),'marker','none','ShowBaseLine','off');
    ylabel(['D' num2str(kk-1)]);
    axis tight;
end
subplot(5,1,5);
stem(cfsmatrix(:,end),'marker','none','ShowBaseLine','off');
ylabel('A3'); xlabel('Sample');
axis tight;

Because the critically sampled wavelet filter bank downsamples the data at each level, the analysis must stop when you have only one coefficient left. In the case of the ECG signal with 2048 samples, this must occur when $L = \log_2{2048}$.

[C,L] = wavedec(wecg,log2(numel(wecg)),'sym4');
fprintf('The number of coefficients at the final level is %d. \n',L(1));
The number of coefficients at the final level is 1. 

If you wish to implement an orthogonal wavelet filter bank without downsampling, you can use modwt.

ecgmodwt = modwt(wecg,'sym4',3);
ecgmra = modwtmra(ecgmodwt,'sym4');
subplot(5,1,1);
plot(wecg); title('Original Signal');

title('MODWT-Based Multiresolution Analysis');
for kk = 2:4
    subplot(5,1,kk)
    plot(ecgmra(kk-1,:));
    ylabel(['D' num2str(kk-1)]);
end
subplot(5,1,5);
plot(ecgmra(end,:));
ylabel('A3'); xlabel('Sample');

In a biorthogonal filter bank, the synthesis filters are not simply time-reversed versions of the analysis filters. The family of biorthogonal spline wavelet filters are an example of such filter banks.

[LoD,HiD,LoR,HiR] = wfilters('bior3.5');

If you examine the analysis filters (LoD,HiD) and the synthesis filters (LoR,HiR), you see that they are very different. These filter banks still provide perfect reconstruction.

[A,D] = dwt(wecg,LoD,HiD);
xrec = idwt(A,D,LoR,HiR);
max(abs(wecg-xrec))
ans =

   4.4409e-16

Biorthogonal filters are useful when linear phase is a requirement for your filter bank. Orthogonal filters cannot have linear phase with the exception of the Haar wavelet filter. If you have the Signal Processing Toolbox software, you can look at the phase responses for an orthogonal and biorthogonal pair of wavelet filters.

[Lodb6,Hidb6] = wfilters('db6');
[PHIdb6,W] = phasez(Hidb6,1,512);
PHIbior35 = phasez(HiD,1,512);
figure;
subplot(2,1,1)
plot(W./(2*pi),PHIdb6); title('Phase Response for db6 Wavelet');
grid on;
xlabel('Cycles/Sample'); ylabel('Radians');
subplot(2,1,2)
plot(W./(2*pi),PHIbior35); title('Phase Response for bior3.5 Wavelet');
grid on;
xlabel('Cycles/Sample'); ylabel('Radians');

Set the dwtmode back to the original setting.

dwtmode(origmodestatus,'nodisplay');