function definition and help file made coherent + unit test added
[scilab.git] / scilab / modules / cacsd / tests / unit_tests / calfrq.tst
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
13 h=syslin('c',1/%s);n=1;
14
15
16 [f,bnds,split]=calfrq(h,0.01,100);
17 if split<>1 then pause,end
18 rf=freq(h.num,h.den,2*%pi*%i*f);
19 //finite difference derivative estimate
20 ddf=diff(f)/pas;
21 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
22 //error between computed and extrapolated value
23 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
24 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
25
26 h=syslin('c',[1/(%s+0.5);1/(%s+0.3)]);n=2;
27 [f,bnds,split]=calfrq(h,0.01,100);
28 if split<>1 then pause,end
29 rf=freq(h.num,h.den,2*%pi*%i*f);
30 //finite difference derivative estimate
31 ddf=diff(f)/pas;
32 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
33 //error between computed and extrapolated value
34 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
35 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
36
37
38
39 h=syslin('c',(%s^2+2*0.9*10*%s+100)/(%s^2+2*0.3*10.1*%s+102.01));
40 //h=h*syslin('c',(%s)/(%s^2+81)) ;
41 n=1;
42 [f,bnds,split]=calfrq(h,0.01,100);
43 if size(split,1)<>1 then pause,end
44
45 rf=freq(h.num,h.den,2*%pi*%i*f);
46 //finite difference derivative estimate
47 ddf=diff(f)/100;
48 drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
49 //error between computed and extrapolated value
50 e=rf(:,2:$)-(rf(:,1:$-1)+drf*100);
51 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
52
53
54
55
56
57
58 //case with singularity
59 h=syslin('c',(%s)/(%s^2+81)) ;
60 sing=9/(2*%pi);
61 n=1;
62 [f,bnds,split]=calfrq(h,0.01,100);
63 if size(split,2)<>2 then pause,end
64 ks=split(2);
65 if abs(f(ks-1)-f(ks))*2*%pi<Epss then pause,end
66 if f(ks-1)>sing|f(ks)<sing then pause,end
67
68 //finite difference derivative estimate
69 f1=f(1:ks-15);//remove points near singularity for which pasmin constraint may be active
70 rf=freq(h.num,h.den,2*%pi*%i*f1);
71 ddf=diff(f1)/pas;
72 drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
73 //error between computed and extrapolated value
74 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
75 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
76 f1=f(ks+15:$);//remove points near singularity for which pasmin constraint may be active
77 rf=freq(h.num,h.den,2*%pi*%i*f1);
78 ddf=diff(f1)/pas;
79 drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
80 //error between computed and extrapolated value
81 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
82 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
83
84
85 //state space case
86 sl=tf2ss(h);
87 n=1;
88 [f1,bnds1,split1]=calfrq(sl,0.01,100);
89 if norm(f-f1)>1d-13 then pause,end
90 if or(split<>split1) then pause,end
91
92 //discrete case
93 h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
94 hd=ss2tf(dscr(h,0.01));
95 [f,bnds,split]=calfrq(hd,0.00001,3);
96 if size(split,1)<>1 then pause,end
97
98 rf=freq(h.num,h.den,exp(2*%pi*%i*f));
99 //finite difference derivative estimate
100 ddf=diff(f)/pas;
101 drf=(freq(h.num,h.den,exp(2*%pi*%i*(f(1:$-1)+ddf)))-rf(:,1:$-1));
102 //error between computed and extrapolated value
103 e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
104 if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end