piir

///// known function case
z=%z;
ts = 1;                // sample time
fz=(z+0.4)*(z-0.2)/(z-0.6)/(z+0.1);    //target function
fs0=0.01*1/ts/2;       // start freq.
fsn=0.49*1/ts/2;       // stop freq.
ff = fs0+(fsn-fs0)*(0:20)/20;
zf = exp(%i*ff'*2*%pi*ts);
qf = horner(fz,zf);
at =[ff',abs(qf),imag(log(qf))/%pi*180];  //  prepare data in the file

if MSDOS then unix('del foo');
else unix('rm -f foo'); end
write('foo',at)

px = piir(ts,2,'foo',0.01);  // The order of px is 2. Phase error is 0.01

// px is almost same as fz

//// arbitrary gain IIR filter design

z=%z;
ts = 1;                // sample time
at = [                 // gain spec
      0.01,0.01,0;     // gain is 0.01 at 0.01Hz
      0.1,1    ,0;     //         0.1  at 0.1
      0.2,2    ,0;
      0.3,1.0  ,0;
      0.4,0.01 ,0;
      0.49,0.01,0];
if MSDOS then unix('del foo');
else unix('rm -f foo'); end
write('foo',at)

q1=log_fir(ts,80,'foo',.03);

px = piir(ts,14,'apgain.txt',0.1);  // The order of px is 14. Phase error is 0.1

// px is the filter having specified gain

///// approximation of exp(1/z)
z=%z;
ts = 1;                // sample time
fz=1/z;    //target function
fs0=0.01*1/ts/2;       // start freq.
fsn=0.49*1/ts/2;       // stop freq.
ff = fs0+(fsn-fs0)*(0:20)/20;
zf = exp(%i*ff'*2*%pi*ts);
qf = exp(horner(fz,zf));
at =[ff',abs(qf),imag(log(qf))/%pi*180];  //  prepare data in the file

if MSDOS then unix('del foo');
else unix('rm -f foo'); end
write('foo',at)

px = piir(ts,2,'foo',0.01);  // The order of px is 2. Phase error is 0.01

//  compare the result
qa=horner(px,zf);
clf()
bode(ff,[qf.';qa.'])
//

///  discretization from a continuos transfer function.
//      ( not from state space )
z=%z;
s=%s
ts = 2;                // sample time
fns=syslin('c',1/(s+0.2)/(s+0.4));  //target function
fn=syslin('c',fns/s);  //s part of target function*(1-1/z)/s
fs0=0.01*1/ts/2;       // start freq.
fsn=0.99*1/ts/2;       // stop freq.
ff = fs0+(fsn-fs0)*(0:20)/20;
i0=(-20:20)/ts

fd=[]
for i=1:size(ff,'*')
  fh = ff(1,i)+i0
  a0 = repfreq(fn,fh)

  zf = exp(%i*ff(1,i)*2*%pi*ts);
  qf = horner(1-1/z,zf);
  fd = [fd;qf*sum(a0)/ts]
end

at =[ff',abs(fd),imag(log(fd))/%pi*180];  //  prepare data in the file

if MSDOS then unix('del foo');
else unix('rm -f foo'); end
write('foo',at)

px = piir(ts,2,'foo',0.01,1);  // The order of px is 2. Phase error is 0.01

dx=ss2tf(dscr(fns,ts))

clf()

bode([px;dx],0.01,1/ts/2,0.001)

// px is almost same as dx

//  the target function is stable except a predetermined factor
//  but not the minimum phase function
z=%z;
ts = 1;                // sample time
fz=(z+4)*(z-0.2)/(z-1)/(z+0.1);    //target function
fs0=0.01*1/ts/2;       // start freq.
fsn=0.99*1/ts/2;       // stop freq.
ff = fs0+(fsn-fs0)*(0:20)/20;
zf = exp(%i*ff'*2*%pi*ts);
qf = horner(fz,zf);
at =[ff',abs(qf),imag(log(qf))/%pi*180];  //  prepare data in the file

if MSDOS then unix('del foo');
else unix('rm -f foo'); end
write('foo',at)

px = piir(ts,-2,'foo',0.001,0,1/(z-1));  // the function has a predetermined factor z-1. 

// px is almost same as fz