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