Filter Realizations - Introduction to Digital Signal Processing and Filter Design - page 336

[b,a]=cheby2(N,As,Wn);
end
if ftype==4
disp('Elliptic Lowpass Filter')
[N,Wn]=ellipord(Wp,Ws,Ap,As);
M=mod(N,2);
if M==0
N=N+1
end
[b,a]=ellip(N,Ap,As,Wn);
end
[h0,w]=freqz(b,a,256);
H0=abs(h0);
plot(w/pi,H0);grid
axis([0.0 1.0 0.0 1.0])
title('MAGNITUDE OF SPECIFIED LP FILTER')
ylabel('Magnitude')
xlabel('Normalized frequency')
% TO FIND Q(z)
k=sum(a)/sum(b);
b=b*k;
fliped a= fliplr(a);
%R(z)= Q2(z)=P2(z)-z^-N D(z^-1)D(z)
R=conv(b,b)-conv(a,fliped a);
% Calculate Q
Q(1)=R(1)^(0.5);
Q(2)=R(2)/(2*Q(1));
for n=2:N
term=0;
for k=1:n-1
term=Q(k+1)*Q(n-k+1)+term;
end
Q(n+1)=(R(n+1)-term)/(2*Q(1));
end
%Zeros of P+Q is calculated
j=1;
k=0;
P plus Q=b+Q;
zeros=roots(P plus Q);
for i=1:N
if abs(zeros(i))<1
zero in(j)=zeros(i);
j=j+1;
else
k=k+1;