%%
% by Ljiljana Milic
% Supplemental material for Chapter XII:
%% 12. Examples of Multirate Filter Banks
% Table of Contents
% 12. Examples of Multirate Filter Banks
%     12.1 TWO-CHANNEL FILTER BANKS
%     12.1.1 Quadrature Mirror (QMF) Filter Banks
%         FIR QMF Filter Bank
%         Example 12.1
%     12.1.2 Alias-Free Two-Channel IIR QMF Banks with Magnitude Preserving Property
%         Example 12.2
%     12.1.3 Orthogonal Two-Channel FIR Filter Banks
%
%         Example 12.3
%         Example 12.4
%     12.2 TREE-STRUCTURED FILTER BANKS
%     12.2.1 Filter banks with equal passband widths
%         Example 12.5
%     12.2.2 Octave Filter Banks
%         Example 12.6
%         Example 12.7
%% 12.1 TWO-CHANNEL FILTER BANKS
%
% Remainder
% The block diagram represenring the analysis/synthesis two-channel filter bank with the processing unit between the analysis and synthesis parts is shown in the Figure. The analysis part consists of the lowpass/highpass filter pair H_0 \left(z\right)/ H_1 \left(z\right) and factor-of-2 down-samplers, and the synthesis  part consists of factor-of-2 up-samplers and the lowpass/highpass filter pair  G_0 \left(z\right)/ G_1 \left(z\right).
%
%
% The two-channel analysis/synthesis filter bank.
%
% Since, the v_0 \left\lbrack n\right\rbrackand v_1 \left\lbrack n\right\rbrack are the down-sampled version of x_0 \left\lbrack n\right\rbrackand x_1 \left\lbrack n\right\rbrack, the z-transforms V_0 \left(z\right) and V_1 \left(z\right) are expressible in terms of the z-transformors with the down-sampling factor M = 2:
% V_0 \left(z\right)=\frac{1}{2}\left\lbrack H_0 \left(z^{1/2} \right)X\left(z^{1/2} \right)+H_0 \left({-z}^{1/2} \right)X\left({-z}^{1/2} \right)\right\rbrack
% V_1 \left(z\right)=\frac{1}{2}\left\lbrack H_1 \left(z^{1/2} \right)X\left(z^{1/2} \right)+H_1 \left({-z}^{1/2} \right)X\left({-z}^{1/2} \right)\right\rbrack
% The first terms represent the z-transforms of the desired decimated signal components, whereas the second terms represent the aliasing components that overlap in the basebands of the decimated signals.
% In the synthesis part, signals w_0 \left\lbrack n\right\rbrackand w_1 \left\lbrack n\right\rbrack are up-sampled-by-2, and up-sampled signals u_0 \left\lbrack n\right\rbrackand u_1 \left\lbrack n\right\rbrackare processed by the lowpass filters G_0 \left(z\right) and G_1 \left(z\right). The z-transforms Y_0 \left(z\right) and Y_1 \left(z\right) are given by
% Y_0 \left(z\right)=G_0 \left(z\right)W_0 \left(z^2 \right)
% Y_1 \left(z\right)=G_1 \left(z\right)W_1 \left(z^2 \right)
% With the assumption that the analysis and synthesis parts are connected together, i.e.,w_0 \left\lbrack n\right\rbrack =v_0 \left(n\right) ,w_1 \left\lbrack n\right\rbrack =v_1 \left(n\right), one obtains the relation between z-transforms Y_0 \left(z\right), Y_1 \left(z\right) and the z-transform of the input signal X\left(z\right).
% Y_0 \left(z\right)=\frac{1}{2}\left\lbrack H_0 \left(z\right)X\left(z\right)+H_0 \left(-z\right)X\left(-z\right)\right\rbrack G_0 \left(z\right)
% Y_1 \left(z\right)=\frac{1}{2}\left\lbrack H_1 \left(z\right)X\left(z\right)+H_1 \left(-z\right)X\left(-z\right)\right\rbrack G_1 \left(z\right)
% Finally, the filtered signals y_0 \left\lbrack n\right\rbrackand y_1 \left\lbrack n\right\rbrack are added together to yield the output y\left\lbrack n\right\rbrack. The z-transform Y\left(z\right) is the sum
% Y\left(z\right)=Y_0 \left(z\right)+Y_1 \left(z\right)
% The above relation is usually expressed in the form:
% Y\left(z\right)=T\left(z\right)X\left(z\right)+A\left(z\right)X\left(X\left(-z\right)\right)
% where
% T\left(z\right)=\frac{1}{2}\left\lbrack H_0 \left(z\right)G_0 \left(z\right)+H_1 \left(z\right)G_1 \left(z\right)\right\rbrack
% is the distortion transfer function , and
% A\left(z\right)=\frac{1}{2}\left\lbrack H_0 \left(-z\right)G_0 \left(z\right)+H_1 \left(-z\right)G_1 \left(z\right)\right\rbrack
% is called the aliasing transfer function.
% It is well known from the literature, that the alias free two-channel filter bank can be achieved by the proper combination of  H_0 \left(z\right), H_1 \left(z\right), G_0 \left(z\right) and G_1 \left(z\right). See Chapter XII of the book.
% Perfect reconstruction property
% The T\left(z\right) represents the the distortion of the overall analysis/ synthesis filter bank. In the frequency domain, we consider the frequency characteristic  T\left(e^{j\omega } \right)
% T\left(e^{j\omega } \right)=\left|T\left(e^{j\omega } \right)\right|e^{-\textrm{jarg}\left\lbrack T^{\left(j\omega \right)} \right\rbrack }
% where \left|T\left(e^{j\omega } \right)\right| is the amplitude distortion, and \arg \left\lbrack T^{j\omega } \right\rbrack phase distortion.
% The magnitude-preserving property means that \left|T\left(e^{j\omega } \right)\right| satisfies
% \left|T\left(e^{j\omega } \right)\right|=\textrm{const},\;\textrm{for}\;\textrm{all}\;
% The phase-preserving property is satisfied if  T\left(z\right) is a linear-phase transfer function
% \arg \left\lbrack T^{j\omega } \right\rbrack =K\omega ,\;\;\textrm{for}\;\textrm{all}\;\omega
% The filter bank satisfies the perfect-reconstraction property if the magnitude-preserving and phase-preserving properties are satisfied. The perfect reconstruction means that the signal at the output of the analysis/synthesis filter bank is y\left\lbrack n\right\rbrack is delayed version of the original signal x\left\lbrack n\right\rbrack
% y\left\lbrack n\right\rbrack =x\left\lbrack n-K\right\rbrack
% In many applications, instead the perfect-reconstrauction property it is sufficient to achieve the nearly-perfect reconstruction property property
% y\left\lbrack n\right\rbrack \approx x\left\lbrack n-K\right\rbrack
%
%% 12.1.1 Quadrature Mirror (QMF) Filter Banks
%
% Remainder
% The Quadrature mirror (QMF) filter bank denotes the mirror symetry of the lowpass/highpass halfband filters H_0 \left(z\right) and H_1 \left(z\right)
% H_1 \left(z\right)=H_0 \left(-z\right)
% For cansecing the aliasing transfer function A\left(z\right), the synthesis filters G_0 \left(z\right) and G_1 \left(z\right) are determined by
% G_0 \left(z\right)=2H_0 \left(-z\right),  and G_1 \left(z\right)=-2H_1 \left(z\right)
% The distortion thansfer function T\left(z\right) is expressible in terms of lowpass  filter H_0 \left(z\right)
% T\left(z\right)=H_0^2 \left(z\right)-H_1^2 \left(z\right)=H_0^2 \left(z\right)-H_0^2 \left(-z\right)
% Accordingly, the proper choice of the lowpass filter H_0 \left(z\right) determines the QMF bank performances.
%%
%% FIR QMF Filter Bank
%% Example 12.1
% In this example, we choose the Johnson's FIR filter of the length N = 12 to examine the performances of a linear-phase FIR QMF bank. The next Figure presents the polyphase structure that is used for implementation.
%
%

clear all, close all

% First, we design the 12th order linear-phase halfband filter H_0 \left(z\right) using the first half of the Johnston's coefficients

B1=[-0.006443977,0.02745539,-0.00758164,-0.0913825,0.09808522,0.4807962];

% Generating the lowpass filter H_0 \left(z\right)

h0 = [B1,fliplr(B1)];     
figure (1)
subplot(2,1,1)
stem(0:length(h0)-1,h0);title('Impulse response of low pass filter')
axis([0,11,-0.6,0.6])
xlabel('n')
% Generate the complementary highpass filter

% Generating the highpass filter H_1 \left(z\right)

for k = 1:length(h0)      
   h1(k) = ((-1)^k)*h0(k);
end
figure (1)
subplot(2,1,2)
stem(0:length(h1)-1,h1);title('Impulse response of high pass filter')
axis([0,11,-0.6,0.6])
xlabel('n')

% Compute the gain responses of the analysis filters

[H1z,w]=freqz(h0,1,256);
H0=abs(H1z); g1=20*log10(H0);
[H2z,w]=freqz(h1,1,256); 
H1=abs(H2z); g2=20*log10(H1);
figure (2)

% Plot the gain responses of the analysis filters

subplot(4,1,1:3)
plot(w/pi,g1,'-',w/pi,g2,'-b'),grid
title('Gain responses of H_0(z) and H_1(z)')
ylabel('Gain, dB')
axis([0,1,-60,2])

% Compute the magnitude response of the distortion transfer function

ver=H0.^2+H1.^2;
figure (2)
subplot(4,1,4)
plot(w/pi,10*log10(ver))
xlabel('Normalized frequency \omega/\pi'), ylabel('Gain, dB')
title('Distortion transfer function')

% Generate the rectangular test signal x

x = [zeros(size(1:100)),ones(size(101:250)),zeros(size(251:511))]; % Test signal

% Polyphase decomposition

e0 = h0(1:2:length(h0));         % Polyphase component E_0(z)
e1 = h0(2:2:length(h0));         % Polyphase component E_1(z)

% Signal analysis

x0 = x(1:2:length(x))/2;   % Down-sampling, branch E_0
x1 = [0,x(2:2:length(x)-1)]/2;  % Down-sampling, branch E_1
v00 = filter(e0,1,x0);         % Filtering with E_0(z)
v11 = filter(e1,1,x1);         % Filtering with E_1(z)
v0 = v00 + v11; % Output of the lowpass filter
v1 = v00 - v11; % Output of the highpass filter
figure (3)
subplot(7,1,1:2)
plot(1:511,x), ylabel('x[n]')
title('Signal decomposition and reconstruction')
axis([0,512,-0.1,1.1])
subplot(7,1,3:4)
plot(0:255,v0),ylabel('v_0[n]')
axis([0,256,-0.1,1.1])
subplot(7,1,5)
plot(0:255,v1),,ylabel('v_1[n]')
axis([0,256,-0.1,0.1])

% Signal synthesis

w00 = v0 + v1; % Input to polyphase component E_1(z), synthesis part
w11 = v0 - v1; % Input to polyphase component E_0(z), synthesis part
u00 = 2*filter(e1,1,w00);   % Filtering with E_1(z)
u11 = 2*filter(e0,1,w11);   % Filtering with E_0(z)	
y0 = zeros(size(1:512));
y0(1:2:512) = u00; % Upsampling, branch E_0
y1 = zeros(size(1:512));
y1(1:2:512) = u11; % Upsampling, branch E_1
y0 = [0,y0(1:511)]; % Inserting the delay z^(-1)
y = 2*(y0 + y1); % Reconstructed signal
figure (3)
subplot(7,1,6:7)
plot(0:511,y), xlabel('Time index n'), ylabel('y[n]')
axis([0,512,-0.1,1.1])
disp('END OF EXAMPLE 12.1')

%
%% 12.1.2 Alias-Free Two-Channel IIR QMF Banks with Magnitude Preserving Property
%
%% Example 12.2
% In this example, we choose the 5th order Butterworth filter to examine the performances of the IIR QMF bank. The next Figure presents the polyphase structure that is used for implementation.
%
%

clear all, close all

% Butterworth filter design

N = 5;  % Selecting the Butterworth filter order
[z,p,k]=butter(N,0.5);     % Butterworth filter design

% Computing and displaying the magnitude responses of the filter pair H_0 \left(z\right), H_1 \left(z\right)

figure (1)
subplot(4,1,1:3)
B=k*poly(z); A=poly(p);     % Lowpass filter H0(z)
[H0,f]=freqz(B,A,256,2);    % Lowpass filter frequency response
plot(f,abs(H0))
hold on
BB=k*poly(-z);   AA = A;           % Highpass filter H1(z)
[H1,f]=freqz(BB,AA,256,2);     % Highpass filter frequency response
plot(f,abs(H1),'b'),ylabel('Magnitude'), grid
title('Magnitude responses of H_0(z) and H_1(z)')
axis([0,1,0,1.1])

% Compute the magnitude response of the distortion transfer function

ver=abs(H0).^2+abs(H1).^2;
figure (1)
subplot(4,1,4)
%plot(f,10*log10(ver))
plot(f,ver)
xlabel('\omega/\pi'), % ylabel('Gain, dB')
title('Magnitude distortion')

% Computing and displaying the group delay of the overall analysis/synthesis QMF bank


[gd,f] = grpdelay(conv(BB,B),conv(AA,A),512,2);
figure (2)
subplot(2,1,1)
plot(f,gd), grid
title('Group delay of the analysis/synthesis two-channel filter bank')
xlabel('\omega/\pi'),ylabel('Group delay, samples')
axis([0,1,0,15])

% Composing the allpass branches A_0 \left(z\right)and A_1 \left(z\right)

c=abs(p).^2;               % Squared moduli of the filter poles
g=sort(c(1:2:length(c)));  
bet0=g(2:2:length(g));    % Selecting beta coefficients for A0(z)
bet1=g(3:2:length(g));    % Selecting beta coefficients for A1(z)
P0=poly(-bet0); % Allpass branch A_0(z)=Q0(z)/P0(z);
Q0=fliplr(P0);  % Allpass branch A_0(z)=Q0(z)/P0(z);
P1=poly(-bet1); % Allpass branch A_1(z)=Q1(z)/P1(z);
Q1=fliplr(P1); % Allpass branch A_1(z)=Q1(z)/P1(z);

% Generating the test signal x

x = [zeros(size(1:100)),ones(size(101:250)),zeros(size(251:511))]; 

% Signal decomposition

xx=x/2;  % Scaling the input signal by 1/2
x0=xx(1:2:length(xx));  % Down-sampling, branch A0
x1=[0,xx(2:2:length(xx))]; % Down-sampling, branch A1
v00=filter(Q0,P0,x0); % Filtering with A_0(z)
v11=filter(Q1,P1,x1); % Filtering with A_1(z)
v0=v00+v11;  % Lowpass output
v1=v00-v11;  % Highpass output

% Signal reconstruction

w00=v0+v1; % Input to A_0(z)
w11=v0-v1; % Input to B_0(z)
u00=filter(Q1,P1,w00); % Filtering with A_1(z)
u11=filter(Q0,P0,w11); % Filtering with A_0(z)
y0=zeros(size(1:512));
y0(1:2:512)=u00;      % Up-sampling, branch A_1
y1=zeros(size(1:512));
y1(1:2:512)=u11;  % Up-sampling, branch B_1
y0=[0,y0(1:511)];  %  Inserting delay z^(-1)
y=y0+y1;  % Reconstructed signal

% Displaying decomposition and reconstruction of the test signal x

figure (3)
subplot(7,1,1:2)
plot(1:511,x), ylabel('x[n]')
title('Signal decomposition and reconstruction')
axis([0,512,-0.3,1.3])
subplot(7,1,3:4)
plot(0:255,v0),ylabel('v_0[n]')
axis([0,256,-0.3,1.3])
subplot(7,1,5)
plot(0:255,v1),,ylabel('v_1[n]')
axis([0,256,-0.3,0.3])
subplot(7,1,6:7)
plot(0:511,y), xlabel('Time index n'), ylabel('y[n]')
axis([0,512,-0.3,1.3])
disp('END OF EXAMPLE 12.2')

%
%% 12.1.3 Orthogonal Two-Channel FIR Filter Banks
%
% Remainder
% The analysis/synthesis filters are nonlinear-phase halfband FIR filters whose transfer functions H_0 \left(z\right), H_1 \left(z\right), G_0 \left(z\right), G_1 \left(z\right) are related in such a manner that the overall analysis/synthesis bank is an alias-free filter bank satisfying the perfect-reconstruction property.
% The transfer functions of the lopass filters H_0 \left(z\right) and  G_0 \left(z\right), [ \textrm{highpass}\;\textrm{filters}\;H_1 \left(z\right) and G_1 \left(z\right)] satisfy the mirror-image symetry
% G_0 \left(z\right)=2z^{-\left(N-1\right)} H_0 \left(z^{-1} \right),    G_1 \left(z\right)=2z^{-\left(N-1\right)} H_1 \left(z^{-1} \right)
% where N is the filter length, an oddnumber.
% The transfer function of the lowpass/highpass analysis filter [H_0 \left(z\right), H_1 \left(z\right)] of the orthogonal filter bank are related by
% H_1 \left(z\right)=-z^{-\left(N-1\right)} H_0 \left(-z^{-1} \right)
% The analysis/synthesis filters are related as follows
% G_0 \left(z\right)=2H_1 \left(-z\right)=2z^{-\left(N-1\right)} H_0 \left(z^{-1} \right)
% G_1 \left(z\right)=2H_0 \left(-z\right)=2z^{-\left(N-1\right)} H_1 \left(z^{-1} \right)
%%
%% Example 12.3
% In this example, we design the orthogonal two-channel filter bank with the analysis filter pair [\left.H_0 \left(z\right),H_1 \left(z\right)\right\rbrack and the synthesis pair [\left.G_0 \left(z\right),G_1 \left(z\right)\right\rbrack]. The filter length is N = 18, and the passband edge frequency of lowpass halfband filters is \omega_p =0\ldotp 4\pi.
%

clear all, close all

% The MATLAB function firpr2chfb returns the impulse responses h_0 \left\lbrack n\right\rbrack ,h_1 \left\lbrack n\right\rbrack ,g_0 \left\lbrack n\right\rbrack ,g_1 \left\lbrack n\right\rbrack \;\;for the input parameters

N = 18; fp = 0.4;

% Computing the impulse responses h_0 \left\lbrack n\right\rbrack ,h_1 \left\lbrack n\right\rbrack ,g_0 \left\lbrack n\right\rbrack ,g_1 \left\lbrack n\right\rbrack \;\;

[h0,h1,g0,g1] = firpr2chfb(N-1,0.4); 

% Display the impulse responses h_0 \left\lbrack n\right\rbrack ,h_1 \left\lbrack n\right\rbrack ,g_0 \left\lbrack n\right\rbrack ,g_1 \left\lbrack n\right\rbrack \;\;

n=0:N-1;
figure (1)
subplot(4,1,1)
stem(n,h0)
title('Impulse responses h_0[n], h_1[n], g_0[n], g_1[n]')
ylabel('h_0[n]')
subplot(4,1,2)
stem(n,h1)
ylabel('h_1[n]')
subplot(4,1,3)
stem(n,g0)
ylabel('g_0[n]')
subplot(4,1,4)
stem(n,g1)
ylabel('g_1[n]')
xlabel('Time index n')

% Computie and plot the gain responses of the analysis filter pair [\left.H_0 \left(z\right),H_1 \left(z\right)\right\rbrack

[H0,f]=freqz(h0,1,1024,2);
[H1,f]=freqz(h1,1,1024,2);
figure (2)
plot(f,20*log10(abs(H0)),f,20*log10(abs(H1)),'b','LineWidth',1); grid
title('Two-channel filter bank')
xlabel('\omega/\pi'), ylabel('Gain, dB')
axis([0,1,-60,2])

% Verification of the perfect-recponstruction property of the analysis/synthesis filter bank

t = (conv(h0,g0)+conv(h1,g1))/2;
figure(3)
stem(0:2*(N-1),t)
xlabel('Time index n'), ylabel('t[n]')
title('Impulse response of the analysis/synthesis filter bank' )
axis([0,2*(N-1),-0.2,1.1])
disp('END OF EXAMPLE 12.3')

%
%% Example 12.4
%
% In this example, we use the function wfilters from the Wavelet Toolbox to construct four filters of the analysis/synthesis bank. We show the design and performances of the two-channel orthogonal filter bank based on Daubechies' db5 wavelet.
% For the analysis/synthesis filters, we use notations from the Wavelet Toolbox:
% Lo_D, Hi_D for the analysis filters  H_0 \left(z\right),H_1 \left(z\right)
% Lo_R, Hi_R for the synthesis filtersG_0 \left(z\right),G_1 \left(z\right)

clear all, close all

% Set the wavelet name

wname = 'db5'; 
% Compute the four filters associated with wavelet name given 
% by the input string wname. 

% Design the four filters associated with the selected wavelet name

[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname); 

% Plot the impulse responses of decomposition and reconstruction filters

N=length(Lo_D) % Filter length
figure (1)
subplot(221); stem(0:N-1,Lo_D); 
title('Decomposition lowpass filter'); 
subplot(222); stem(0:N-1,Hi_D); 
title('Decomposition highpass filter'); 
subplot(223); stem(0:N-1,Lo_R); xlabel('n')
title('Reconstruction lowpass filter'); 
subplot(224); stem(0:N-1,Hi_R); xlabel('n')
title('Reconstruction highpass filter') 

% Display the pole/zero locations of decomposition and reconstruction filters

figure (2)
subplot(2,2,1)
zplane(Lo_D,1)
title('Decomposition lowpass filter'); 
subplot(2,2,2)
zplane(Hi_D,1)
title('Decomposition highpass filterThe'); 
subplot(2,2,3)
zplane(Lo_R,1)
title('Reconstruction lowpass filter'); 
subplot(2,2,4)
zplane(Hi_R,1)
title('Reconstruction highpass filter'); 
disp('The pole/zero plots illustrste the following properties:')
disp('The 5th order zero at the point z=-1 for the lowpassfilters')
disp('The 5th order zero at at the point z=+1 for the highpass filters.')
disp('The analysis part is composed of:')
disp('      -the maximum-phase lowpass filter and')
disp('      -the minimum-phase highpass filter')
disp('The synthesis part is composed of:')
disp('      -the minimum-phase lowpass filter and')
disp('      -the maximum-phase highpass filter.')


% Computie and plot the magnitude responses of the analysis filter pair [\left.H_0 \left(z\right),H_1 \left(z\right)\right\rbrack

[HLo_D,f]=freqz(Lo_D,1,512,2);
[HHi_D,f]=freqz(Hi_D,1,512,2);
figure (3)
plot(f,abs(HLo_D))
title('Tho-channel filter bank')
xlabel('Normalized frequency \omega/\pi'), ylabel ('Magnitude')
hold on
plot(f,abs(HHi_D),'b'),grid


% Impulse response of the analysis/synthesis filter bank

N=length(Lo_D);
n=0:N-1;
figure (4)
stem(0:2*(N-1),conv(Lo_D,Lo_R)+conv(Hi_D,Hi_R))
title('Impulse response of the analysis/synthesis filter bank')
xlabel('Time index n'), ylabel('t[n]')
axis([0,18,-0.3,2.2])
disp('END OF EXAMPLE 12.4');

%% 12.2 TREE-STRUCTURED FILTER BANKS
%
% Remainder
% An approach for constructing a multichannel filter bank is based on a tree structure, which uses the two-channel filter banks as building blocks. With the use of this approach, the multichannel filter banks with uniform and nonuniform separation between the channels can be generated. If the two-channel filter bank satisfy the perfect-reconstruction (nearly-perfect reconstruction) property, the overall tree-structured filter bank also satisfies the perfect-reconstruction (nearly-perfect reconstruction) property.
% For the sake of simplicity, we use in this section the symbolic representations of the analysis and synthesis two-channel filter banks as shown in the Figure. The single-input/two-output device A^{\left(k\right)} \left(z\right) symbolises a two-channel analysis filter bank, and two-input/single-output device S^{\left(k\right)} \left(z\right) symbolises a two-channel synthesis bank.
%
%
% Symbolic representation of two-channel bank: (a) Analysis bank. (b) Synthesis bank.
%
%% 12.2.1 Filter banks with equal passband widths
%
% Remainder
% A two-level four-channel analysis/synthesis filter bank is indicated in the next Figure, which illustrates how the tree structured filter bank is developed. The output y[n] represents the reconstructed signal. The reconstructed signal is delayed replica of the input signal x[n] under the following conditions:
% The processing unit is omited from the structure, i.e., outputs of the analysis bank are connected directly to the inputs of the synthesis bank.
% The analysis and synthesis two-channel filter banks A^{\left(k\right)} \left(z\right) and S^{\left(k\right)} \left(z\right) satisfy the perfect-reconstruction property.
%
%
% Example of the two-level four-channel filter bank
%
% The single-stage equivalent of the four channel filter bank is shown in the next Figure. Recol that this single stage equivalent which coresponds to the four channel uniform filter bank is developed from the Third and Sixth Identities, see Examples 2.5 and 2.6.
%
% Single-stage equivalent of the two-level four-channel filter bank
%%
%
%% Example 12.5
% In this example, we use the two-channel orthogonal filter bank of Example 12.3 to construct the uniform four-channel filter bank using the tree structure. Here, the MATLAB function firpr2chfb computes the coefficients for the analysis/synthesis filters of the two-channel filter banks.

clear all, close all

% Setting the filter lengths

N=18

% Design of the two-channel filter bank for perfect reconstruction

[h0,h1,f0,f1]=firpr2chfb(N-1,0.4);

% Computing and plotting magnitude responses of the analysis filter pair

[H0,f]=freqz(h0,1,1024,2);
[H1,f]=freqz(h1,1,1024,2);
figure (1)
plot(f,abs(H0),f,abs(H1),'b');
title('Two-channel bank')
xlabel('Normaliyed frequency')
ylabel('Magnitude')
axis([0,1,0,1.1])

% Computing the impulse responses of the  four-channel analysis filter bank

hh0 = conv(h0,upsample(h0,2));  % Impulse response, 2nd level 
hh1 = conv(h0,upsample(h1,2));  % Impulse response, 2nd level 
hh2 = conv(h1,upsample(h0,2));  % Impulse response, 2nd level 
hh3 = conv(h1,upsample(h1,2));  % Impulse response, 2nd level 

% Analysis filters frequency responses

% Computing the frequency responses of the four-channel analysis filter bank

[HH0,f] = freqz(hh0,1,1024,2);  
[HH1,f] = freqz(hh1,1,1024,2);
[HH2,f] = freqz(hh2,1,1024,2);
[HH3,f] = freqz(hh3,1,1024,2);


% Computing and plotting the gain responses of the four-channel analysis filter bank

figure (2)
subplot(4,1,1:3)
plot(f,20*log10(abs(HH0)),f,20*log10(abs(HH1)),f,20*log10(abs(HH2)),'g',...
    f,20*log10(abs(HH3)))
ylabel('Gain, dB')
axis([0,1,-50,2])

% Computing the magnitude distortion function \left|T\left(e^{j\omega } \right)\right|

T= abs(HH0).^2+abs(HH1).^2+abs(HH2).^2+abs(HH3).^2; % Verification of the power-complementary property
figure (2)
subplot(4,1,4)
plot(f,T), xlabel('Normalized frequency \omega/\pi'),ylabel('|T(e^j^\omega)|')
axis([0,1,0.9,1.1])

% Comment:
% The bottom subfigure shows that  \left|T\left(e^{j\omega } \right)\right| is a constant and therefore, confirms the magnitude-preserving property of the overall analysis/synthesis four-channel filter bank.
% The two-channel filter banks are linear-phase FIR filters , the overall filter bank is also of linear phase.
% Accordingli, the four-channel filter bank satisfies the perfect-reconstruction property.

disp('END OF EXAMPLE 12.5')

%
%
%% 12.2.2 Octave Filter Banks
%
% Remainder
% An octave filter bank is a multilevel filter bank generated of the two-channel filter banks. Development process of an analysis bank is indicated in the Figure bellow, which displayes an example of the four-level five-channel analysis filter bank. The next Figure illustrates a typical magnitude response of an octave filter bank.
%
% Generation of the orthogonal four-level five-channel analysis filter bank.
%
% Magnitude responses of the example  five-channel orthogonal filter bank.The octave bank is constructed by using the orthogonal two-channel filter bank of Example 12.3.
%
% The structure of the synthesis octave bank for the case of four-level five-channel filter bank is indicated in the next Figure. The overall bank is composed of the two-channel synthesis filter banks.
%
% Generation of the orthogonal four-level five-channel synthesis filter bank.
%
% In two examples that follow, Example 12.6 and Exanple 12.7, we demonstrate how the process of signal decomposition and reconstruction is developed through the four-level analysis/synthesis octave bank.For the investigation purposes, the synthesis bank is connected to the analysis bank through delays that are appropriately entered between them, as shown in the Figure.
%
% Analysis/synthesis four-level octave bank.
%
%% Example 12.6
% In this example, we use the orthogonal two-channel filter bank of Example 12.3 to construct the four-level five-channel analysis/synthesis octave bank. The filter bank is used to perform the decomposition and reconstruction of the rectangular shape signal.l

clear all, close all

% Design the analysis/synthesis ortthogonal two-channel filter bank using the MATLAB function firpr2chfb for designing the analysis and synthesis filters H_0 \left(z\right), H_1 \left(z\right), G_0 \left(z\right)and G_1 \left(z\right) of the length N = 18.

N =18;  % Setting the filter length

% Computing the coefficients of the analysis and synthesis filters.

[h0,h1,g0,g1] = firpr2chfb(N-1,0.4);

% Generating and displaying the test signal x

x=[zeros(size(1:20)),ones(size(21:50)),zeros(size(51:511))];
figure (1)
subplot(3,1,3)
plot(x), title('Test signal'),xlabel('Time index n')
axis([0,512,-0.2,1.2])

% Analysis part: signal decomposition
% Level 1;

x0 = filter(h0,1,x);
x1 = filter(h1,1,x);
v0 = downsample(x0,2);
v1 = downsample(x1,2);

% Level 2

x2 = v0;
x02 = filter(h0,1,x2);
x12 = filter(h1,1,x2);
v02 = downsample(x02,2);
v12 = downsample(x12,2);
v2 = v12;

% Level 3

x3 = v02;
x03 = filter(h0,1,x3);
x13 = filter(h1,1,x3);
v03 = downsample(x03,2);
v13 = downsample(x13,2);
v3 = v13;

% Level 4

x4 = v03;
x04 = filter(h0,1,x4);
x14 = filter(h1,1,x4);
v04 = downsample(x04,2);
v14 = downsample(x14,2);
v4 = v14;
v5 = v04;

% Displaying the results of the signal decomposition

figure (2)
subplot(5,1,1),plot(v1)
ylabel('v_1[n]'),axis([0,256,-0.2,0.2])
title('Signal decomposition')
subplot(5,1,2)
plot(v2)
ylabel('v_2[n]'),axis([0,256,-0.2,0.2])
subplot(5,1,3)
plot(v3)
ylabel('v_3[n]'),axis([0,256,-0.2,0.2])
subplot(5,1,4)
plot(v4)
ylabel('v_4[n]'),axis([0,256,-0.2,0.2])
subplot(5,1,5)
plot(v5)
ylabel('v_5[n]'),xlabel('Time index n')
axis([0,256,-0.2,1.2])

% Synthesis part: process of signal reconstruction
% Level 4

w04=upsample(v5,2);
w14=upsample(v4,2);
y3=filter(g0,1,w04) + filter(g1,1,w14);

% Level 3

w13 = [zeros(size(1:N-1)),v3(1:length(v3)-(N-1))];
w13 = upsample(w13,2);
yu3 = upsample(y3,2);
y2 = filter(g0,1,yu3) + filter(g1,1,w13);

% Level 2

w12 = [zeros(size(1:3*(N-1))),v2(1:length(v2)-3*(N-1))];
w12 = upsample(w12,2);
yu2 = upsample(y2,2);
y1 = filter(g0,1,yu2) + filter(g1,1,w12);

% Level 1

w11 = [zeros(size(1:7*(N-1))),v1(1:length(v1)-7*(N-1))];
w11 = upsample(w11,2);
yu1 = upsample(y1,2);

% Reconstructed signal

y = filter(g0,1,yu1) + filter(g1,1,w11); % Reconstructed signal

% Displaying the process of the signal reconstruction

figure (3)
subplot(4,1,1)
plot(y3)
ylabel('y_3(n)'),axis([0,512,-0.2,1.2])
title('Signal reconstruction')
subplot(4,1,2)
plot(y2)
ylabel('y_2[n]'),axis([0,512,-0.2,1.2])
subplot(4,1,3)
plot(y1)
ylabel('y_1[n]'),axis([0,512,-0.2,1.2])
subplot(4,1,4)
plot(y)
ylabel('y[n]'),xlabel('Time index n')
axis([0,512,-0.2,1.2])
length(y);
disp('=================================================================')
disp('Reconstructed signal y is a delayed replica of the test signal x.')
disp('=================================================================')
disp('END OF EXAMPLE 12.6')

%
%% Example 12.7
% In this example, we design  orthogonal two-channel filter bank using the Daubeshies db9 wavelet, and compose the four- level five-channel analysis/synthesis octave bank. The MATLAB function wfilters(wname), with wname = 'db9', computes the coefficients of the analysis/synthesis filters H_0 \left(z\right), H_1 \left(z\right), G_0 \left(z\right)and G_1 \left(z\right). We use this filter bank to perform the decomposition and reconstruction of the rectangular shape signal.

clear all, close all

% Generating four filters for the orthogonal filter bank
% Set wavelet name:

wname = 'db9'; 

% Design the four filters associated with wavelet name given by the input string wname.

[h0,h1,g0,g1] = wfilters(wname); 
N = 18; % The length of db9 filter

% Compute and plot magnitude responses of the analysis filters H_0 \left(z\right) and H_1 \left(z\right)

[H0,f]=freqz(h0,1,512,2);
[H1,f]=freqz(h1,1,512,2);
figure (1)
plot(f,(abs(H0))/sqrt(2),'b',f,(abs(H1))/sqrt(2),'b'), grid
xlabel('Normalized frequency \omega/\pi'), ylabel('Magnitude')
axis([0,1,0,1.1])



% Generating and displaying test signal x

x=[zeros(size(1:20)),ones(size(21:50)),zeros(size(51:511))];
figure (2)
subplot(3,1,3)
plot(x)
title('Test signal')
xlabel('Time index n'),axis([0,512,-0.2,1.2])


% Analysis part: signal decomposition
% Level 1

x0 = filter(h0,1,x);
x1 = filter(h1,1,x);
v0 = downsample(x0,2);
v1 = downsample(x1,2);

% Level 2

x2 = v0;
x02 = filter(h0,1,x2);
x12 = filter(h1,1,x2);
v02 = downsample(x02,2);
v12 = downsample(x12,2);
v2 = v12;

% Level 3

x3 = v02;
x03 = filter(h0,1,x3);
x13 = filter(h1,1,x3);
v03 = downsample(x03,2);
v13 = downsample(x13,2);
v3 = v13;

% Level 4

x4 = v03;
x04 = filter(h0,1,x4);
x14 = filter(h1,1,x4);
v04 = downsample(x04,2);
v14 = downsample(x14,2);
v4 = v14;
v5 = v04;


% Displaying the results of the signal decomposition

figure (3)
subplot(5,1,1)
plot(v1)
ylabel('v_1[n]'),axis([0,256,-0.5,0.5])
title('Signal decomposition')
subplot(5,1,2)
plot(v2)
ylabel('v_2[n]'),axis([0,256,-0.5,0.5])
subplot(5,1,3)
plot(v3)
ylabel('v_3[n]'),axis([0,256,-0.5,1])
subplot(5,1,4)
plot(v4)
ylabel('v_4[n]'),axis([0,256,-0.5,0.5])
subplot(5,1,5)
plot(v5)
ylabel('v_5[n]'),xlabel('Time index n')
axis([0,256,-0.5,4])

% Synthesis part: Signal reconstruction
% Level 4

w04=upsample(v5,2);
w14=upsample(v4,2);
y3=filter(g0,1,w04) + filter(g1,1,w14);

% Level 3

w13 = [zeros(size(1:N-1)),v3(1:length(v3)-(N-1))];
w13 = upsample(w13,2);
yu3 = upsample(y3,2);
y2 = filter(g0,1,yu3) + filter(g1,1,w13);

% Level 2

w12 = [zeros(size(1:3*(N-1))),v2(1:length(v2)-3*(N-1))];
w12 = upsample(w12,2);
yu2 = upsample(y2,2);
y1 = filter(g0,1,yu2) + filter(g1,1,w12);

% Level 1

w11 = [zeros(size(1:7*(N-1))),v1(1:length(v1)-7*(N-1))];
w11 = upsample(w11,2);
yu1 = upsample(y1,2);

% Reconstructed signal y

y = filter(g0,1,yu1) + filter(g1,1,w11); % Reconstructed signal

% Displaying the process of signal reconstruction

figure (4)
subplot(4,1,1)
plot(y3)
ylabel('y_3(n)'), axis([0,512,-0.4,3.4])
title('Signal reconstruction')
subplot(4,1,2)
plot(y2)
ylabel('y_2[n]'),axis([0,512,-0.4,2.4])
subplot(4,1,3)
plot(y1)
ylabel('y_1[n]'),axis([0,512,-0.4,2.4])
subplot(4,1,4)
plot(y)
ylabel('y[n]'),xlabel('Time index n')
axis([0,512,-0.2,1.2])
disp('=================================================================')
disp('Reconstructed signal y is a delayed replica of the test signal x.')
disp('=================================================================')
disp('END OF EXAMPLE 12.7')
disp('                       END OF CHAPTER 12')