error(number): converting occurrences remaining in all .sce .sci files
[scilab.git] / scilab / modules / cacsd / macros / calfrq.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 1985 - 2016 - INRIA - Serge Steer
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13
14 function [frq, bnds, splitf] = calfrq(h, fmin, fmax)
15     //!
16
17     eps = 1.d-14   //minimum absolute lower frequency
18     k = 0.001      // Minimum relative prediction error in the nyquist plan
19     epss = 0.002   // minimum frequency distance with a singularity
20     nptmax = 5000  //maximum number of discretization points
21     tol = 0.01     // Tolerance for testing pure imaginary numbers
22
23     // Check inputs
24     // ------------
25     if and(typeof(h) <> ["state-space" "rational" "zpk"]) then
26         args=["h", "fmin", "fmax"]
27         ierr=execstr("%"+overloadname(h)+"_calfrq("+strcat(args(1:rhs),",")+")","errcatch")
28         if ierr<>0 then
29             msg = _("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n")
30             error(msprintf(msg, "calfrq", 1))
31         end
32         return
33     end
34
35
36     if typeof(h) == "state-space" then
37         h = ss2tf(h)
38     elseif typeof(h) == "zpk" then
39         h=zpk2tf(h)
40     end
41
42     [m, n] = size(h.num)
43     dom = h("dt")
44     select dom
45     case "d" then
46         dom = 1
47     case [] then
48         msg = _("%s: Argument #%d: Undefined time domain.\n");
49         error(msprintf(msg, "calfrq", 1));
50     case 0 then
51         msg = _("%s: Argument #%d: Undefined time domain.\n");
52         error(msprintf(msg, "calfrq", 1));
53     end;
54
55     if type(dom) == 1 then
56         nyq_frq = 1/2/dom;
57         if fmax > nyq_frq then
58             msg = _("%s: Frequencies beyond Nyquist frequency are ignored.\n")
59             warning(msprintf(msg, "calfrq"));
60             fmax = min(fmax, nyq_frq)
61         end
62         if fmin < -nyq_frq then
63             msg = _("%s: Negative frequencies below Nyquist frequency are ignored.\n")
64             warning(msprintf(msg, "calfrq"));
65             fmin = max(fmin, -nyq_frq)
66         end
67     end
68     // Use symmetry to reduce the range
69     // --------------------------------
70     if fmin < 0 & fmax >= 0 then
71         [frq, bnds, splitf] = calfrq(h, eps, -fmin)
72         ns1 = size(splitf, "*")-1;
73         nsp = size(frq, "*");
74         bnds = [bnds(1), bnds(2), -bnds(4), -bnds(3)];
75
76         if fmax > eps then
77             if fmax == -fmin then
78                 splitf = [1, (nsp+2)*ones(1,ns1)-splitf($:-1:2), nsp*ones(ns1)+splitf(2:$)];
79                 bnds = [bnds(1), bnds(2), min(bnds(3), -bnds(3)), max(bnds(4), -bnds(4))];
80                 frq = [-frq($:-1:1), frq]
81             else
82                 [frq2, bnds2, splitf2] = calfrq(h, eps, fmax);
83                 ns2 = size(splitf2,"*")-1
84                 splitf = [1, (nsp+2)*ones(1,ns1)-splitf($:-1:2), nsp*ones(ns2)+splitf2(2:$)];
85                 bnds = [min(bnds(1), bnds2(1)), max(bnds(2), bnds2(2)),...
86                 min(bnds(3), bnds2(3)), max(bnds(4), bnds2(4))];
87                 frq = [-frq($:-1:1), frq2]
88             end
89             return
90         else
91             frq = -frq($:-1:1);
92             nsp = size(frq, "*");
93
94             splitf = [1, (nsp+2)*ones(1, ns1)-splitf($:-1:2)]
95             bnds = bnds;
96             return;
97         end
98     elseif fmin < 0 & fmax <= 0 then
99         [frq, bnds, splitf] = calfrq(h, -fmax, -fmin)
100         ns1 = size(splitf, "*")-1;
101         frq = -frq($:-1:1);
102         nsp = size(frq, "*");
103         splitf = [1, (nsp+2)*ones(1, ns1)-splitf($:-1:2)]
104         bnds = [bnds(1), bnds(2), -bnds(4), -bnds(3)];
105         return;
106     elseif fmin >= fmax then
107         msg = _("%s: Wrong value for input arguments #%d and #%d: %s < %s expected.\n");
108         error(msprintf(msg, "calfrq", 2, 3, "fmin", "fmax"));
109     end
110
111     // Compute dicretisation over a given range
112     // ----------------------------------------
113
114
115     splitf = []
116     if fmin == 0 then fmin = min(1d-14, fmax/10);end
117     //
118     denh = h("den"); numh = h("num")
119     l10 = log(10)
120
121     // Locate singularities to avoid them
122     // ----------------------------------
123     if dom == "c" then
124         c = 2*%pi;
125         // selection function for singularities in the frequency range
126         deff("f=%sel(r, fmin, fmax, tol)",["f = [],";
127         "if prod(size(r)) == 0 then return, end";
128         "f = imag(r(find((abs(real(r))<=tol*abs(r))&(imag(r)>=0))))";
129         "if f <> [] then f = f(find((f>fmin-tol)&(f<fmax+tol))); end"]);
130     else
131         c = 2*%pi*dom
132         // selection function for singularities in the frequency range
133         deff("[f] = %sel(r, fmin, fmax, dom, tol)",["f = [],";
134         "if prod(size(r)) == 0 then return, end";
135         "f = r(find( ((abs(abs(r)-ones(r)))<=tol)&(imag(r)>=0)))";
136         "if f <> [] then ";
137         "  f = atan(imag(f), real(f)); nf = prod(size(f))";
138         "  for k=1:nf,";
139         "    kk = int((fmax-f(k))/(2*%pi))+1;";
140         "    f = [f; f(1:nf)+2*%pi*kk*ones(nf, 1)];";
141         "  end;"
142         "  f = f(find((f>fmin-tol)&(f<fmax+tol)))";
143         "end"]);
144     end
145
146     sing = [];zers = [];
147     fmin = c*fmin, fmax = c*fmax;
148
149     for i=1:m
150         sing = [sing; %sel(roots(denh(i), "e"), fmin, fmax, tol)];
151     end
152
153     pp = gsort(sing', "g", "i"); npp = size(pp, "*");//'
154
155     // singularities just on the left of the range
156     kinf = find(pp<fmin)
157     if kinf <> [] then
158         fmin = fmin+tol
159         pp(kinf) = []
160     end
161
162     // singularities just on the right of the range
163     ksup = find(pp>=fmax)
164     if ksup <> [] then
165         fmax = fmax-tol
166         pp(ksup) = []
167     end
168
169     // check for nearly multiple singularities
170     if pp <> [] then
171         dpp = pp(2:$)-pp(1:$-1)
172         keq = find(abs(dpp)<2*epss)
173         if keq <> [] then pp(keq) = [], end
174     end
175
176     if pp <> [] then
177         frqs = [fmin real(matrix([(1-epss)*pp; (1+epss)*pp], 2*size(pp, "*"), 1)') fmax]
178         //'
179     else
180         frqs = [fmin fmax]
181     end
182     nfrq = size(frqs, "*");
183
184     // Evaluate bounds of nyquist plot
185     //-------------------------------
186     xt = []; Pas = []
187     for i=1:2:nfrq-1
188         w = logspace(log(frqs(i))/log(10), log(frqs(i+1))/log(10), 100);
189         xt = [xt, w]
190         Pas = [Pas w(2)-w(1)]
191     end
192     if dom == "c" then
193         rf = freq(h("num"), h("den"), %i*xt);
194     else
195         rf = freq(h("num"), h("den"), exp(%i*xt));
196     end
197     //
198     xmin = min(real(rf)); xmax = max(real(rf));
199     ymin = min(imag(rf)); ymax = max(imag(rf));
200     bnds = [xmin xmax ymin ymax];
201     dx = max([xmax-xmin, 1]); dy = max([ymax-ymin, 1]);
202
203     // Compute discretization with a step adaptation method
204     // ----------------------------------------------------
205     frq = [];
206     i = 1;
207     nptr = nptmax; // number of unused discretization points
208     l10last = log10(frqs($));
209     while i<nfrq
210         f0 = frqs(i); fmax = frqs(i+1);
211         while f0==fmax do
212             i = i+2;
213             f = frqs(i); fmax = frqs(i+1);
214         end
215         frq = [frq, f0];
216         pas = Pas(floor(i/2)+1)
217         splitf = [splitf size(frq, "*")];
218
219         f = min(f0+pas, fmax);
220
221         if dom == "c" then // Continuous case
222             while f0<fmax
223                 rf0 = freq(h("num"), h("den"), (%i*f0));
224                 rfc = freq(h("num"), h("den"), %i*f);
225                 // compute prediction error
226                 epsd = pas/100; // epsd = 1.d-8
227
228                 rfd = (freq(h("num"), h("den"), %i*(f0+epsd))-rf0)/(epsd);
229                 rfp = rf0+pas*rfd;
230
231                 e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
232                 if e > k then rf0 = freq(h("num"), h("den"), (%i*f0));
233                     rfc = freq(h("num"), h("den"), %i*f);
234                     // compute prediction error
235                     epsd = pas/100; // epsd = 1.d-8
236
237                     rfd = (freq(h("num"), h("den"), %i*(f0+epsd))-rf0)/(epsd);
238                     rfp = rf0+pas*rfd;
239
240                     e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
241                     // compute minimum frequency logarithmic step to ensure a maximum
242                     //of nptmax points to discretize
243                     pasmin = f0*(10^((l10last-log10(f0))/(nptr+1))-1)
244                     pas = pas/2
245                     if pas < pasmin then
246                         pas = pasmin
247                         frq = [frq, f]; nptr = max([1, nptr-1])
248                         f0 = f; f = min(f0+pas, fmax)
249                     else
250                         f = min(f0+pas, fmax)
251                     end
252                 elseif e < k/2 then
253                     pas = 2*pas
254                     frq = [frq, f]; nptr = max([1, nptr-1])
255                     f0 = f; f = min(f0+pas, fmax),
256                 else
257                     frq = [frq, f];nptr = max([1, nptr-1])
258                     f0 = f; f = min(f0+pas, fmax),
259                 end
260             end
261         else  // Discrete case
262             pas = pas/dom
263             while f0<fmax
264                 rf0 = freq(h("num"), h("den"), exp(%i*f0))
265                 rfd = dom*(freq(h("num"), h("den"), exp(%i*(f0+pas/100)))-rf0)/(pas/100);
266                 rfp = rf0+pas*rfd
267                 rfc = freq(h("num"), h("den"), exp(%i*f));
268                 e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
269                 if e > k then
270                     pasmin = f0*(10^((l10last-log10(f0))/(nptr+1))-1)
271                     pas = pas/2
272                     if pas < pasmin then
273                         pas = pasmin
274                         frq = [frq, f]; nptr = max([1, nptr-1])
275                         f0 = f; f = min(f0+pas, fmax)
276                     else
277                         f = min(f0+pas, fmax)
278                     end
279                 elseif e < k/2 then
280                     pas = 2*pas
281                     frq = [frq, f]; nptr = max([1, nptr-1])
282                     f0 = f; f = min(f0+pas, fmax),
283                 else
284                     frq = [frq, f]; nptr = max([1, nptr-1])
285                     f0 = f; f = min(f0+pas, fmax),
286                 end
287             end
288         end
289         i = i+2
290     end
291     frq(size(frq, "*")) = fmax
292     frq = frq/c;
293 endfunction