1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA -
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution. The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
11 function [frq,bnds,splitf]=calfrq(h,fmin,fmax)
14 eps=1.d-14 //minimum absolute lower frequency
15 k=0.001; // Minimum relative prediction error in the nyquist plan
16 epss=0.002 // minimum frequency distance with a singularity
17 nptmax=5000 //maximum number of discretization points
18 tol=0.01 // Tolerance for testing pure imaginary numbers
22 if and(typeof(h)<>["state-space" "rational"])
23 error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"calfrq",1))
25 if typeof(h)=="state-space" then
43 warning(msprintf(gettext("%s: Frequencies beyond Nyquist frequency are ignored.\n"),"calfrq"));
44 fmax=min(fmax,nyq_frq)
47 warning(msprintf(gettext("%s: Negative frequencies below Nyquist frequency are ignored.\n"),"calfrq"));
48 fmin=max(fmin,-nyq_frq)
51 // Use symmetry to reduce the range
52 // --------------------------------
53 if fmin<0&fmax>=0 then
54 [frq,bnds,splitf]=calfrq(h,eps,-fmin)
55 ns1=size(splitf,"*")-1;
57 bnds=[bnds(1),bnds(2),-bnds(4),-bnds(3)];
61 splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2),nsp*ones(ns1)+splitf(2:$)];
62 bnds=[bnds(1),bnds(2),min(bnds(3),-bnds(3)),max(bnds(4),-bnds(4))];
63 frq=[-frq($:-1:1),frq]
65 [frq2,bnds2,splitf2]=calfrq(h,eps,fmax);
66 ns2=size(splitf2,"*")-1
67 splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2),nsp*ones(ns2)+splitf2(2:$)];
68 bnds=[min(bnds(1),bnds2(1)),max(bnds(2),bnds2(2)),...
69 min(bnds(3),bnds2(3)),max(bnds(4),bnds2(4))];
70 frq=[-frq($:-1:1),frq2]
77 splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2)]
81 elseif fmin<0&fmax<=0 then
82 [frq,bnds,splitf]=calfrq(h,-fmax,-fmin)
83 ns1=size(splitf,"*")-1;
86 splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2)]
87 bnds=[bnds(1),bnds(2),-bnds(4),-bnds(3)];
89 elseif fmin >= fmax then
90 error(msprintf(gettext("%s: Wrong value for input arguments #%d and #%d: %s < %s expected.\n"),..
91 "calfrq",2,3,"fmin","fmax"));
94 // Compute dicretisation over a given range
95 // ----------------------------------------
99 if fmin==0 then fmin=min(1d-14,fmax/10);end
101 denh=h("den");numh=h("num")
104 // Locate singularities to avoid them
105 // ----------------------------------
108 //selection function for singularities in the frequency range
109 deff("f=%sel(r,fmin,fmax,tol)",["f=[],";
110 "if prod(size(r))==0 then return,end";
111 "f=imag(r(find((abs(real(r))<=tol*abs(r))&(imag(r)>=0))))";
112 "if f<>[] then f=f(find((f>fmin-tol)&(f<fmax+tol)));end"]);
115 //selection function for singularities in the frequency range
116 deff("[f]=%sel(r,fmin,fmax,dom,tol)",["f=[],";
117 "if prod(size(r))==0 then return,end";
118 "f=r(find( ((abs(abs(r)-ones(r)))<=tol)&(imag(r)>=0)))";
120 " f=atan(imag(f),real(f));nf=prod(size(f))";
122 " kk=int((fmax-f(k))/(2*%pi))+1;";
123 " f=[f;f(1:nf)+2*%pi*kk*ones(nf,1)];";
125 " f=f(find((f>fmin-tol)&(f<fmax+tol)))";
130 fmin=c*fmin,fmax=c*fmax;
133 sing=[sing;%sel(roots(denh(i),"e"),fmin,fmax,tol)];
136 pp=gsort(sing',"g","i");npp=size(pp,"*");//'
138 // singularities just on the left of the range
145 // singularities just on the right of the range
152 //check for nearly multiple singularities
154 dpp=pp(2:$)-pp(1:$-1)
155 keq=find(abs(dpp)<2*epss)
156 if keq<>[] then pp(keq)=[],end
160 frqs=[fmin real(matrix([(1-epss)*pp;(1+epss)*pp],2*size(pp,"*"),1)') fmax]
167 // Evaluate bounds of nyquist plot
168 //-------------------------------
171 w=logspace(log(frqs(i))/log(10),log(frqs(i+1))/log(10),100);
176 rf=freq(h("num"),h("den"),%i*xt);
178 rf=freq(h("num"),h("den"),exp(%i*xt));
181 xmin=min(real(rf));xmax=max(real(rf));
182 ymin=min(imag(rf));ymax=max(imag(rf));
183 bnds=[xmin xmax ymin ymax];
184 dx=max([xmax-xmin,1]);dy=max([ymax-ymin,1]);
186 // Compute discretization with a step adaptation method
187 // ----------------------------------------------------
190 nptr=nptmax; // number of unused discretization points
191 l10last=log10(frqs($));
193 f0=frqs(i);fmax=frqs(i+1);
196 f=frqs(i);fmax=frqs(i+1);
199 pas=Pas(floor(i/2)+1)
200 splitf=[splitf size(frq,"*")];
204 if dom=="c" then //cas continu
206 rf0=freq(h("num"),h("den"),(%i*f0));
207 rfc=freq(h("num"),h("den"),%i*f);
208 // compute prediction error
209 epsd=pas/100;//epsd=1.d-8
211 rfd=(freq(h("num"),h("den"),%i*(f0+epsd))-rf0)/(epsd);
214 e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
215 if (e>k) then rf0=freq(h("num"),h("den"),(%i*f0));
216 rfc=freq(h("num"),h("den"),%i*f);
217 // compute prediction error
218 epsd=pas/100;//epsd=1.d-8
220 rfd=(freq(h("num"),h("den"),%i*(f0+epsd))-rf0)/(epsd);
223 e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
224 // compute minimum frequency logarithmic step to ensure a maximum
225 //of nptmax points to discretize
226 pasmin=f0*(10^((l10last-log10(f0))/(nptr+1))-1)
230 frq=[frq,f];nptr=max([1,nptr-1])
231 f0=f;f=min(f0+pas,fmax)
237 frq=[frq,f];nptr=max([1,nptr-1])
238 f0=f;f=min(f0+pas,fmax),
240 frq=[frq,f];nptr=max([1,nptr-1])
241 f0=f;f=min(f0+pas,fmax),
247 rf0=freq(h("num"),h("den"),exp(%i*f0))
248 rfd=dom*(freq(h("num"),h("den"),exp(%i*(f0+pas/100)))-rf0)/(pas/100);
250 rfc=freq(h("num"),h("den"),exp(%i*f));
251 e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
253 pasmin=f0*(10^((l10last-log10(f0))/(nptr+1))-1)
257 frq=[frq,f];nptr=max([1,nptr-1])
258 f0=f;f=min(f0+pas,fmax)
264 frq=[frq,f];nptr=max([1,nptr-1])
265 f0=f;f=min(f0+pas,fmax),
267 frq=[frq,f];nptr=max([1,nptr-1])
268 f0=f;f=min(f0+pas,fmax),
274 frq( size(frq,"*") )=fmax