%%
% by Ljiljana Milic
% Supplemental material for Chapter VIII:
%% 8. Complementary Filter Pairs
%
% Table of Contents
% 8. Complementary Filter Pairs
%     8.1 FIR Complementary Filter Pairs
%         8.1.1 Delay-Complementary and Magnitude-Complementary FIR Filter Pair
%         Example 8.1
%         8.1.2 Power-Complementary FIR Filter Pairs
%         Example 8.2
%     8.2  IIR complementary filter pairs
%         8.2.1 Class I: Power-Complementary and All-Pass complementary Filter Pairs Implemented as a Parallel Connection of Two All-Pass Subfilters
%         Example 8.3
%         8.2.2 Class II: Power-Complementary Filter Pairs Implemented as a Tapped Cascaded Interconnection of Two Identical All-Pass Subfilters
%         Example 8.4
%%
%% 8.1 FIR Complementary Filter Pairs
%
%% 8.1.1 Delay-Complementary and Magnitude-Complementary FIR Filter Pair
%
% Remainder
%
% Two filters H_{\textrm{LP}} \left(z\right) and H_{\textrm{HP}} \left(z\right) make the delay-complementary filter pair when satisfy
% H_{\textrm{LP}} \left(z\right)+H_{\textrm{HP}} \left(z\right)=z^{\left(N-1\right)/2}
% where N is the filter length, an odd number.
% The magnitude-complementary property is satisfied when H_{\textrm{LP}} \left(z\right) and H_{\textrm{HP}} \left(z\right) are designed to meet the following criteria:
% 1-\delta \le \left|\left(H_{\textrm{LP}} \left(e^{j\omega } \right)\right)\right|\le 1 and \left|\left(H_{\textrm{HP}} \left(e^{j\omega } \right)\right)\right|\le \delta for \omega \in\left\lbrack 0,\omega_p \right\rbrack
% 1-\delta \le \left|\left(H_{\textrm{HP}} \left(e^{j\omega } \right)\right)\right|\le 1 and \left|\left(H_{\textrm{LP}} \left(e^{j\omega } \right)\right)\right|\le \delta for \omega \in\left\lbrack 0,\pi \right\rbrack
% \left|\left(H_{\textrm{LP}} \left(e^{j\omega } \right)\right)\right|=\left|\left(H_{\textrm{HP}} \left(e^{j\omega } \right)\right)\right|=\frac{1}{2}   for \omega =\omega_c
% \delta ={10}^{\frac{{-a}_s }{20}}
% Here, \omega_p and \omega_s are passband/stopband edge frequencies, \omega_c is the crossover frequency and a_s denotes minimal attenuation in the stopband.
%% Example 8.1
% Design magnitude-complementary and delay-complementary halfband FIR filter pair to meet the criteria defined above.
% Filter order N_{\textrm{ord}}=18, the passband/stopband ripple \delta =0\ldotp 04, the crossover frequency \omega_c =0\ldotp 5\pi.
% Solution:

close all, clear all

% Setting the input parameters

Nord = 18;  delta = 0.04; fc = 0.5;   

%  Step 1: Adjusting the ripple value

delp = delta/(2*(1 - delta));          

% Step 2: Design of the intermediate filter Elp(z)

elp = firceqrip(Nord,fc,[delp,delp]);

% Step 3: Modifying the intermediate filter Elp(z)

hlp = (1 - delta)*elp;     

% Step 3: Modifying the intermediate filter Elp(z)

hlp(Nord/2 + 1) = hlp(Nord/2+1)+delta/2;% 

% Zero settings

hhp=zeros(size(1:length(hlp))); 

% Step 4: Complementary highpass filter

hhp(10) = 1;  hhp = hhp - hlp; 

% Display impulse responses

figure 
stem((0:18),hlp), xlabel('n'),ylabel('h_L_P[n]')
title('Lowpass filter impulse response')
axis([0,18,-0.4,0.6])
figure
stem((0:18),hhp), xlabel('n'),ylabel('h_L_P[n]')
title('Highpass filter impulse response')
axis([0,18,-0.4,0.6])

% Display the pole-zero locations of lowpass filter

figure
zplane(hlp,1)
title('Zero-pole locations of lowpass filter')

% Computing the magnitude responses

f(1)=0;
omega(1)=pi*f(1);
K=(length(hlp)+1)/2;
n=1:K-1;
h_red=hlp(1:K-1);
h_red=fliplr(h_red);
hh_red=hhp(1:K-1);
hh_red=fliplr(hh_red);
for p=1:101
Alp(p)=1/2+2*sum(h_red.*cos(omega(p)*n));
Ahp(p)=1/2+2*sum(hh_red.*cos(omega(p)*n));
f(p+1)=f(p)+0.01;
omega(p+1)=pi*f(p+1);
end

% Displaying magnitude responses of the filter pair

figure
plot(f(1:length(Alp)),Alp,f(1:length(Ahp)),Ahp,'b')
xlabel('Normalized frequency \omega/\pi'),ylabel('Magnitude'),grid
title('Complementary filter pair')
axis([0,1,-0.1,1.1])
disp('END OF EXAMPLE 8.1')

%
%% 8.1.2 Power-Complementary FIR Filter Pairs
%
% The power-complementary property is satisfied when the sum of squared magnitude responses of the lowpass and highpass filters of the pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack is equal to unity for all frequencies,
% {\left|H_{\textrm{LP}} \left(e^{j\omega } \right)\right|}^2 +{\left|H_{\textrm{HP}} \left(e^{j\omega } \right)\right|}^2 =1
% The power-complementary property is achieved when the filter magnitude responses meet the following criteria:
% \sqrt{1-\delta^2 }\le \left|H_{\textrm{LP}} \left(e^{j\omega } \right)\right|\le 1 and \left|H_{\textrm{HP}} \left(e^{j\omega } \right)\right|\le \delta  for  \omega \in \left\lbrack {0,\omega }_p \right\rbrack
% \sqrt{1-\delta^2 }\le \left|H_{\textrm{HP}} \left(e^{j\omega } \right)\right|\le 1 and \left|H_{\textrm{LP}} \left(e^{j\omega } \right)\right|\le \delta  for  \omega \in \left\lbrack \omega_s ,\pi \right\rbrack
% \left|\left(H_{\textrm{LP}} \left(e^{j\omega } \right)\right)\right|=\left|\left(H_{\textrm{HP}} \left(e^{j\omega } \right)\right)\right|=\frac{\sqrt{2}}{2}   for \omega =\omega_c
% \delta ={10}^{\frac{{-a}_s }{20}}
%% Example 8.2
% Design the power-complementary halfband FIR filter pair to meet the specifications given bellow:
% Filter order N_{\textrm{ord}\;} =15, the passband/stopband ripple \delta =0\ldotp 1, the crossover frequency \omega_c =0\ldotp 5\pi.
% Solution:

close all, clear all

% Setting the input papameters

Nord = 15; Nord = 2*Nord;  delta = 0.1;   

% Adjusting the ripple value, Substep 1.1

dfir =(delta^2)/(2*(1 - (delta^2)));    

% Designing the intermediate filter Elp(z), Substep 1.2

elp = firceqrip(Nord,0.5,[dfir,dfir]); 

% Modifying the intermediate filter Elp(z), Substep 1.3

eh = (1-(delta^2))* elp;                                                   

% Modifying the intermediate filter Elp(z), Substep 1.3

eh(Nord/2+1) = eh(Nord/2 + 1) + (delta^2)/2;%          

% Generating the minimum-phase lowpass filter Hlp(z), Step 2

[hlp,ssp,iter] = minphase(eh(Nord/2 + 1:Nord+1));           

%  Indexing the filter coefficients

n = 0:length(hlp) - 1;                                                         

%  Complementary highpass filter, Step 3

hhp = fliplr(hlp).*((-1).^n);                                                  

% Displaying the results

figure 
stem((0:15),hlp), xlabel('n'),ylabel('h_L_P[n]')
title('Lowpass filter impulse response')
axis([0,15,-0.5,0.6])
figure
stem((0:15),hhp),xlabel('n'),ylabel('h_H_P[n]')
title('Highpass filter impulse response')
axis([0,15,-0.5,0.6])

figure
zplane(hlp,1)
title('Lowpass filter, pole-zaro locations')
figure
zplane(hhp,1)
title('Highpass filter, pole-zaro locations')

figure 
[Hlp,f]=freqz(hlp,1,512,2);
[Hhp,f]=freqz(hhp,1,512,2);
plot(f,abs(Hlp),f,abs(Hhp),'b')
xlabel('Normalized frequency \omega/\pi'),ylabel('Magnitude'),grid
title('Magnitude response of the complementary filter pair')
axis([0,1,-0.1,1.1])
g(1)=axes('Position',[0.20 0.70 0.25 0.1]);
plot(f(1:222),(abs(Hlp(1:222)))), axis([0,0.45,0.99,1.002]),grid
g(2)=axes('Position',[0.62 0.70 0.25 0.1]);
plot(f(285:512),(abs(Hhp(285:512)))), axis([0.55,1,0.99,1.002]),grid

disp('END OF EXAMPLE 8.2')

%
%% 8.2  IIR complementary filter pairs
%
%% 8.2.1 Class I: Power-Complementary and All-Pass complementary Filter Pairs Implemented as a Parallel Connection of Two All-Pass Subfilters
%
% Remainder
%
% The lowpass (highpass) filter transfer function G_{\textrm{LP}} \left(z\right)\;\left(G_{\textrm{HP}} \left(z\right)\right) is an odd order minimum-phase IIR transfer function whose magnitude response meets the following criteria
% \sqrt{1-\delta^2 }\le \left|G_{\textrm{LP}} \left(e^{j\omega } \right)\right|\le 1 and \left|G_{\textrm{HP}} \left(e^{j\omega } \right)\right|\le \delta  for  \omega \in \left\lbrack {0,\omega }_p \right\rbrack
% \sqrt{1-\delta^2 }\le \left|G_{\textrm{HP}} \left(e^{j\omega } \right)\right|\le 1 and \left|G_{\textrm{LP}} \left(e^{j\omega } \right)\right|\le \delta  for  \omega \in \left\lbrack \omega_s ,\pi \right\rbrack
% \left|\left(G_{\textrm{LP}} \left(e^{j\omega } \right)\right)\right|=\left|\left(G_{\textrm{HP}} \left(e^{j\omega } \right)\right)\right|=\frac{\sqrt{2}}{2}   for \omega =\omega_c
% \delta ={10}^{\frac{{-a}_s }{20}}
% The implementation structure of the filter pair \left\lbrack G_{\textrm{LP}} \left(z\right),G_{\textrm{HP}} \left(z\right)\right\rbrack is based on the parallel connection of two all-pass subfilters A_0 \left(z\right) and A_1 \left(z\right),
%
% G_{\textrm{LP}} \left(z\right)=\frac{1}{2}\left(A_0 \left(z\right)\right)+\left(A_1 \left(z\right)\right)
% G_{\textrm{HP}} \left(z\right)=\frac{1}{2}\left(A_0 \left(z\right)\right)-\left(A_1 \left(z\right)\right)
%% Example 8.3

close all, clear all

% Setting the input parameters for lowpass filter

N = 5;   fp = 0.4;  

%     Lowpass halfband filter design

[b,a,z,p,k] = halfbandiir(N,fp);                                             

%  Sorting the filter poles

pabs = sort(abs(p)); 

% Coefficients of allpass branches

beta0 = (abs(pabs(2)))^2;    beta1=(abs(pabs(4)))^2;      

% Frequency response of A_0 \left(z\right)

[A0,f] = freqz([beta0,0,1],[1,0,beta0],512,2);    

% Frequency response of A_1 \left(z\right)

[A1,f] = freqz([0,beta1,0,1],[1,0,beta1],512,2);  

%   Lowpass filter frequency response

GLp = (A0 + A1)/2;                                                                 

% Highpass filter frequency response

GHp = (A0 - A1)/2;                                                               

% Displaying the magnitude response of the complementary filter pair

figure 
plot(f,(abs(GLp)),f,abs(GHp),'b')
xlabel ('Normalized frequency \omega/\pi'),ylabel('Magnitude'), grid,
axis([0,1,0,1.05])
title('Complementary IIR filter pair [G_0(z),G_1(z)]')
g(1)=axes('Position',[0.24 0.70 0.25 0.1]);
plot(f(1:205),(abs(GLp(1:205)))), axis([0,0.4,0.9998,1.00002])
g(2)=axes('Position',[0.62 0.70 0.25 0.1]);
plot(f(307:512),(abs(GHp(307:512)))), axis([0.6,1,0.9998,1.00002])
disp('END OF EXAMPLE 8.3')

%
%% 8.2.2 Class II: Power-Complementary Filter Pairs Implemented as a Tapped Cascaded Interconnection of Two Identical All-Pass Subfilters
%
%  Remainder
%
% The filter pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack exhibits power-complementary property as defined for example 8.2.
% The solution for Class II filter pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack is achieved with the aid of the following transfer functions:
% a) FIR prototype filter pair \left\lbrack F_{\textrm{LP}} \left(z\right),F_{\textrm{HP}} \left(z\right)\right\rbrack designed to provide passbans/stopband ripple and power-complementary property of the overall pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack.
% b) IIR prototype filter pair \left\lbrack G_{\textrm{LP}} \left(z\right),G_{\textrm{HP}} \left(z\right)\right\rbrack designed to provide edge frequencies of the pair\left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack.
%
%% Example 8.4
% Design the Class II power-complementary IIR filter pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack to meet the following specifications:
% Minimal attenuation in the stopband a_s =60\;\textrm{dB}. The edge frequencies of the filter pair are: \omega_s =0\ldotp 6\pi and \omega_p =\pi -\omega_s =0\ldotp 4\pi.
% Implementation structure of the Class II power-complementary IIR filter pair:
%
% Solution:
% Functions firceqrip, minphase were used.

close all, clear all

% Design of the FIR prototype filter pair

Nord = 5; Nord = 2*Nord;
delta = 0.001;
dfir =(delta^2)/(2*(1-(delta^2)));
hl = firceqrip(Nord,0.5,[dfir,dfir]); 
eh = (1-(delta^2))*hl;
eh(Nord/2+1) = eh(Nord/2+1)+(delta^2)/2; 
[alp,ssp,iter] = minphase(eh(Nord/2+1:Nord+1)); 
n = 0:length(alp)-1;
ahp = fliplr(alp).*((-1).^n);
len = Nord/2;
figure 
stem((0:len),alp), xlabel('n'),ylabel('a_L[n]')
title('Impulse response of lowpass FIR filter')
axis([0,len,-0.5,0.6])
figure
stem((0:len),ahp),xlabel('n'),ylabel('a_H[n]')
title('Impulse response of highpass FIR filter')
axis([0,len,-0.5,0.6])
figure 
[Hlp,f]=freqz(alp,1,512,2);
[Hhp,f]=freqz(ahp,1,512,2);
plot(f,20*log10(abs(Hlp)),f,20*log10(abs(Hhp)),'b')
xlabel('\omega/\pi'),xlabel('\theta/\pi'),ylabel('Gain [dB]'),grid
title('Gain responses of prototype FIR filter pair')
axis([0,1,-80,2])


% Design of the IIR prototype filter pair

beta = 0.545577589; % For explanation see pages 266-267 in the book.
Aa0 = [1,0,beta];
Ab0 = fliplr(Aa0);
Ab1 = [0,1];
Aa1 = 1;
[H0,f] = freqz(Ab0,Aa0,512,2);
[H1,f] = freqz(Ab1,Aa1,512,2);
Hiil=(H0+H1);
Hiih=(H0-H1);
figure 
subplot(5,1,1:2)
plot(f,20*log10(abs(Hiil/2)),f,20*log10(abs(Hiih/2)),'b'),grid
xlabel('\omega/\pi'),ylabel('Gain [dB]')
title('Gain responses of prototype IIR filter pair')
axis([0,1,-30,2])
a = alp;

% Lattice coefficients of the 5th order minimum phase halfband FIR filter

a5=a/a(1);
b5=a5;
b5(1:2:6)=-a5(1:2:6);
b5=fliplr(b5);
k5=a5(6);
a3=(a5+k5*b5)/(1+k5^2);
a3=a3(1:4);
b3=(b5-k5*a5)/(1+k5^2);
b3=b3(3:6);
k3=a3(4);
a1=(a3+k3*b3)/(1+k3^2);
a1=a1(1:2);
b1=(b3-k3*a3)/(1+k3^2);
b1=b1(3:4);
k1=a1(2);
disp('=================================================')
disp([sprintf(' lattice coefficient k1 = %15.9f', k1)])
disp([sprintf(' lattice coefficient k3 = %15.9f', k3)])
disp([sprintf(' lattice coefficient k5 = %15.9f', k5)])
disp('=================================================')

% Coefficients of  allpass branches

Q=[beta,0,1];
P=fliplr(Q);
QQ=conv(Q,Q);
PP=conv(P,P);

% Computing and plotting the gain response of the overall filter pair \left\lbrack H_{\textrm{LP}} \left(z\right),H_{\textrm{HP}} \left(z\right)\right\rbrack

u=zeros(size(1:1024));
u(1)=1;
v1=[0,1];
v2=[0,0,1];
x1=filter(Q,P,u)-k1*filter(v1,1,u);
y1=k1*filter(Q,P,u)+filter(v1,1,u);
x3=filter(QQ,PP,x1)-k3*filter(v2,1,y1);
y3=k3*filter(QQ,PP,x1)+filter(v2,1,y1);
x5=filter(QQ,PP,x3)-k5*filter(v2,1,y3);
y5=k5*filter(QQ,PP,x3)+filter(v2,1,y3);
yLP = x5;
yHP = y5;
[H0,w]=freqz(yLP,1,4096);
dop=20*log10(max(abs(H0)));
figure 
plot(w/(pi),20*log10(abs(H0))-dop)
hold on
[H1,w]=freqz(yHP,1,4096);
plot(w/(pi),20*log10(abs(H1))-dop,'b'), grid
xlabel('\omega/\pi'),ylabel('Gain [dB]')
title('Complementary filter pair [H_0(z),H_1(z)], gain response')
axis([0,1,-80,2])
disp('Specifications are met with two very low order filters.')
disp('END OF EXAMPLE 8.4')
disp('                     END OF CHAPTER VIII')

%

function [y,ssp,iter]=minphase(h);
% Reproduced with permission of IEEE from: IEEE Trans. Circuits and Systems-I: Fundamental Theory
% and Applications-I: 50(3), 365-375.

y=[1 zeros(1,length(h)-1)];
ssp=realmax;
ss=ssp/2;
iter=0;
d=0;
while (ss<ssp)
    y=y+d';
    ssp=ss;
    iter=iter+1;
    Ar=toeplitz([y(1),zeros(1,length(h)-1)],y);
    Al=fliplr(toeplitz([y(length(h)),zeros(1,length(h)-1)],fliplr(y)));
    A=Al+Ar;
    b=h'-Al*y';
    d=A\b;
    ss=norm(d);
end;
end