function definition and help file made coherent + unit test added
[scilab.git] / scilab / modules / cacsd / tests / unit_tests / calfrq.dia.ref
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - rINRIA - Serge Steer
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 K=0.001;     // Minimum relative prediction error in the nyquist plan
8 Epss=0.002;   // minimum frequency distance with a singularity
9 nptmax=5000;  //maximum number of discretisation points
10 pas=100/(2*%pi);
11 s=%s;
12 h=syslin('c',1/%s);n=1;
13 [f,bnds,split]=calfrq(h,0.01,100);
14 if split<>1 then bugmes();quit;end
15 rf=freq(h.num,h.den,2*%pi*%i*f);
16 //finite difference derivative estimate
17 ddf=diff(f)/pas;
18 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
19 //error between computed and extrapolated value
20 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
21 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
22 h=syslin('c',[1/(%s+0.5);1/(%s+0.3)]);n=2;
23 [f,bnds,split]=calfrq(h,0.01,100);
24 if split<>1 then bugmes();quit;end
25 rf=freq(h.num,h.den,2*%pi*%i*f);
26 //finite difference derivative estimate
27 ddf=diff(f)/pas;
28 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
29 //error between computed and extrapolated value
30 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
31 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
32 h=syslin('c',(%s^2+2*0.9*10*%s+100)/(%s^2+2*0.3*10.1*%s+102.01));
33 //h=h*syslin('c',(%s)/(%s^2+81)) ;
34 n=1;
35 [f,bnds,split]=calfrq(h,0.01,100);
36 if size(split,1)<>1 then bugmes();quit;end
37 rf=freq(h.num,h.den,2*%pi*%i*f);
38 //finite difference derivative estimate
39 ddf=diff(f)/100;
40 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
41 //error between computed and extrapolated value
42 e=rf(:,2:$)-(rf(:,1:$-1)+drf*100);
43 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
44 //case with singularity
45 h=syslin('c',(%s)/(%s^2+81)) ;
46 sing=9/(2*%pi);
47 n=1;
48 [f,bnds,split]=calfrq(h,0.01,100);
49 if size(split,2)<>2 then bugmes();quit;end
50 ks=split(2);
51 if abs(f(ks-1)-f(ks))*2*%pi<Epss then bugmes();quit;end
52 if f(ks-1)>sing|f(ks)<sing then bugmes();quit;end
53 //finite difference derivative estimate
54 f1=f(1:ks-15);//remove points near singularity for which pasmin constraint may be active
55 rf=freq(h.num,h.den,2*%pi*%i*f1);
56 ddf=diff(f1)/pas;
57 drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
58 //error between computed and extrapolated value
59 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
60 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
61 f1=f(ks+15:$);//remove points near singularity for which pasmin constraint may be active
62 rf=freq(h.num,h.den,2*%pi*%i*f1);
63 ddf=diff(f1)/pas;
64 drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
65 //error between computed and extrapolated value
66 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
67 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
68 //state space case
69 sl=tf2ss(h);
70 n=1;
71 [f1,bnds1,split1]=calfrq(sl,0.01,100);
72 if norm(f-f1)>1d-13 then bugmes();quit;end
73 if or(split<>split1) then bugmes();quit;end
74 //discrete case
75 h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
76 hd=ss2tf(dscr(h,0.01));
77 [f,bnds,split]=calfrq(hd,0.00001,3);
78 if size(split,1)<>1 then bugmes();quit;end
79 rf=freq(h.num,h.den,exp(2*%pi*%i*f));
80 //finite difference derivative estimate
81 ddf=diff(f)/pas;
82 drf=(freq(h.num,h.den,exp(2*%pi*%i*(f(1:$-1)+ddf)))-rf(:,1:$-1));
83 //error between computed and extrapolated value
84 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
85 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end