%IIR Lab ECE609 C55x
%Infinity norm scaling to prevent state variable overflow
%Illustrated using Elliptic Filter Example
%lcpotter, Apr2007

close all;clear all
Fs = 48000;  % Sampling Frequency
N      = 6;      % Order
Fpass1 = 9600;   % First Passband Edge Frequency
Fpass2 = 12000;  % Second Passband Edge Frequency
Apass  = 1;      % Passband Ripple (dB)
Astop  = 50;     % Stopband Attenuation (dB)

% Calculate the zpk values using the ELLIP function.
[B,A] = ellip(N/2, Apass, Astop, [Fpass1 Fpass2]/(Fs/2), 'stop');

% To avoid round-off errors, do not use the transfer function.  Instead,
% convert transfer function to second-order sections.
[sos,g] = tf2sos(B,A,'UP','none');
sos(1,1:3) = g*sos(1,1:3);%incorp gain into first monic numerator
nsec=size(sos,1); %number of sections
Nf=1024;%number of frequency samples for inf norm (i.e., max abs)
H(1,:)=freqz(sos(1,1:3),sos(1,4:6),Nf);%first stage
G(1,:) = freqz(1,sos(1,4:6),Nf);%gain to first state variable
s(1) = max(abs(G(1,:)));
pre = 1/s(1);%prescale factor
for k=2:nsec
    H(k,:) = freqz(sos(k,1:3),sos(k,4:6),Nf);
    G(k,:) = freqz(1,sos(k,4:6),Nf).' .* H(k-1,:);
    s(k) = max(abs(G(k,:)));
    scale(k-1) = s(k-1)/s(k);%scale factor for kth section
end
tmp = ceil(max(abs(s(nsec)*sos(nsec,1:3))));%at last numerator
if(tmp > 1), asl = nextpow2(tmp); else asl=0;end;% pull out left shifts
scale(nsec) = s(nsec)/( 2^asl );
%modify coefficient list
coeffs = sos(:,[1 2 3 5 6]);%a0=1, so omit
coeffs(:,4) = coeffs(:,4)/2;%a1 coefficients divide by 2
for k=1:nsec
    coeffs(k,1:3) = coeffs(k,1:3)*scale(k);%scale numerators
end
%Finally, display results
coeffs
pre
asl