Drawing2.jpg

by Ljiljana Milic
Supplemental material for Chapter VIII:

8. Complementary Filter Pairs

Table of Contents

8.1 FIR Complementary Filter Pairs

8.1.1 Delay-Complementary and Magnitude-Complementary FIR Filter Pair

Remainder
Two filters and make the delay-complementary filter pair when satisfy
where N is the filter length, an odd number.
The magnitude-complementary property is satisfied when and are designed to meet the following criteria:
and for
and for
for
Here, and are passband/stopband edge frequencies, is the crossover frequency and 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 =18, the passband/stopband ripple , the crossover frequency .
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')
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 is equal to unity for all frequencies,
The power-complementary property is achieved when the filter magnitude responses meet the following criteria:
and for
and for
for

Example 8.2

Design the power-complementary halfband FIR filter pair to meet the specifications given bellow:
Filter order , the passband/stopband ripple , the crossover frequency .
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')
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 is an odd order minimum-phase IIR transfer function whose magnitude response meets the following criteria
and for
and for
for
The implementation structure of the filter pair is based on the parallel connection of two all-pass subfilters and ,
Figure8_13.jpg

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
[A0,f] = freqz([beta0,0,1],[1,0,beta0],512,2);
Frequency response of
[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')
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 exhibits power-complementary property as defined for example 8.2.
The solution for Class II filter pair is achieved with the aid of the following transfer functions:
a) FIR prototype filter pair designed to provide passbans/stopband ripple and power-complementary property of the overall pair .
b) IIR prototype filter pair designed to provide edge frequencies of the pair.

Example 8.4

Design the Class II power-complementary IIR filter pair to meet the following specifications:
Minimal attenuation in the stopband . The edge frequencies of the filter pair are: and .
Implementation structure of the Class II power-complementary IIR filter pair:
Figure8_16.jpg
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)])
lattice coefficient k1 = 2.382963175
disp([sprintf(' lattice coefficient k3 = %15.9f', k3)])
lattice coefficient k3 = -0.541364466
disp([sprintf(' lattice coefficient k5 = %15.9f', k5)])
lattice coefficient k5 = 0.108528870
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
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.')
Specifications are met with two very low order filters.
disp('END OF EXAMPLE 8.4')
END OF EXAMPLE 8.4
disp(' END OF CHAPTER VIII')
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