dicrete stable IIR filter
fp=piir(t0,nf,file1,tol,[r,[nfx]])
The program finds a nf-order stable discrete IIR filter having specified phase. With log_fir() you can design a stable filter (minimum phase) with arbitrary gain.
sample time
order of the IIR filter. ( in case minus value, gain fitting is performed if the target function is stable but not the minimum phase function.)
an external file name. The file which describes a phase spec.
max. phase error at specified frequencies.(degree)
degree of denominator of fp minus degree of numerator of fp.
The result IIR filter fp has the factor nfx. fp = nfx*xx
dicrete minimum phase IIR filter.
If you like to have a filter fp whoes phase is p1(degree) at f1(Hz), p2(degree) at f2(Hz),.. you create a file "file1". In the file
f1 g1 p1
f2 g2 p2
..
The program uses LMI. If LMI stops, you may need to increase nf. If message ' no n order answer, type abort' appears, you may select nf from even number.
The phase of fp minus the target is less than tol.
See
M.Kisaka: Phase design of an IIR filter and its application to an inverse system, Proceeding of the 2006 IEEE international conference on control system,Munich, Germany,Oct.2619/2622(2006)
///// 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 | ![]() | ![]() |