CACSD: fix freson() Scilab 5 behavior
[scilab.git] / scilab / modules / cacsd / macros / freson.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2000 - 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 function fr=freson(h)
14     [lhs,rhs]=argn(0);
15
16     if argn(2) < 1 then
17         error(msprintf(_("%s: Wrong number of input argument(s): %d expected.\n"),"freson",1));
18     end
19     if and(typeof(h)<>["state-space","rational","zpk"]) then
20         ierr=execstr("[gm,fr]=%"+overloadname(sys)+"_freson(h)","errcatch")
21         if ierr<>0 then
22             error(msprintf(gettext("%s: Wrong type for input argument: Linear dynamical system expected.\n"),"freson",1))
23         end
24         return
25     end
26     if size(h,"*")<>1 then
27         error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"freson",1))
28     end
29
30     if typeof(h)=="state-space" then
31         h=ss2tf(h)
32     elseif typeof(h)=="zpk" then
33         h=zpk2tf(h)
34     end
35     dt=h.dt
36     [n,d]=h(["num","den"]);
37     if type(n)==1 then
38         n=poly(n,varn(d),"c");
39     end
40     if coeff(d,0)==0 then
41         error(msprintf(_("%s: Wrong value for input argument #%d: infinite gain at zero frequency.\n"),"freson",1))
42     end
43
44     //look for  omega such that derivative of magn. is zero
45     if dt=="c" then
46         //find frequencies which zeros the magnitude derivative
47         hh=h*horner(h,-%s)
48         r=roots(derivat(hh).num)
49         k=find(imag(r)>0&abs(real(r))<%eps*abs(r));
50         fr=imag(r(k))/(2*%pi)
51     else
52         if dt=="d"|dt==[] then dt=1;end
53         //find frequencies which zeros the magnitude derivative
54         hh=h*horner(h,1/%z)
55         r=roots(derivat(hh).num);
56         k=find(abs(abs(r)-1)<=sqrt(%eps)*abs(r));
57         r=imag(log(r(k)));
58         fr=r(r>0&r<0.999*%pi)/(2*%pi*dt)
59     end
60     if fr==[] then return;end
61     //find frequencies that correspond to a magnitude peak
62     k=find(abs(repfreq(h,fr))-abs(repfreq(h,fr*0.999))>0)
63     fr=fr(k)
64     if fr==[] then return;end
65     fr=gsort(fr,"g","d");
66 endfunction