www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/OrthogonalAndBiorthogonalFilterBanksExample.m

    %% 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.
%%
% <<../critfilterbank2.png>>
%%
% 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' extremal phase wavelet with 2 vanishing moments. The
% example explains the notion of vanishing moments in a later section.

% Copyright 2015 The MathWorks, Inc.

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');
%%
% 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]
%%
% 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' 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.
%% 
% <<../twolevel_analysis1.png>>
%%
% At each successive level, the number of scaling and wavelet coeffficients
% 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));
%%
% 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))
%%
% 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');