bug 9285 fix
[scilab.git] / scilab / modules / cacsd / macros / g_margin.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA -  Author: Serge Steer
3 // 
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-en.txt
9
10 function [gm,fr]=g_margin(h)
11   //compute the gain margin of a SISO transfer function
12   select typeof(h)
13   case 'rational' then 
14   case 'state-space' then 
15     h=ss2tf(h)
16   else 
17     error(97,1),
18   end;
19   if or(size(h)<>[1 1]) then 
20     error(msprintf(_("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"g_margin",1))
21   end
22   //
23   epsr=1.e-7;//used for testing if complex numbers are real
24   eps1=1.e-7;//used for testing if complex numbers have a modulus near 1
25   epssing=1e-10; //used for testing if arguments are not singular points of h
26   if h.dt=='c' then  //continuous time case
27     // get s such as h(s)=h(-s) and s=iw 
28      s=%i*poly(0,"w");
29      //compute h(s)-h(-s)=num/den
30      num=imag(horner(h.num,s)*conj(horner(h.den,s)))
31      den=real(horner(h.den,s)*conj(horner(h.den,s)))
32      //necessary condition
33      w=roots(num,"e");
34      ws=real(w(abs(imag(w))<epsr&real(w)<=0)) //points where phase is -180°
35                                              
36      //remove nearly singular points     
37      ws(abs(horner(num,ws))>=epssing*abs(horner(den,ws)))=[]
38      if ws==[] then gm=%inf,fr=[],return,end
39      mingain=real(freq(h.num,h.den,%i*ws))
40   else  //discrete time case
41     if h.dt=='d' then dt=1,else dt=h.dt,end
42     //get z such as h(z)=h(1/z) and z=e^(%i*w*dt)
43     //form hh=h(z)-h(1/z)
44     z=poly(0,varn(h.den));
45     sm=simp_mode();simp_mode(%f);hh=h-horner(h,1/z);simp_mode(sm)
46     //find the numerator roots
47     z=roots(hh.num,"e");
48     z(abs(abs(z)-1)>eps1)=[]// retain only roots with modulus equal to 1
49     
50     //remove nearly singular points                                       
51     z(abs(horner(hh.num,z))>=epssing*abs(horner(hh.den,z)))=[];
52
53     w=log(z)/(%i*dt)
54     ws=real(w(abs(imag(w))<epsr)) //points where phase is -180°
55     if ws==[] then gm=%inf,fr=[],return,end
56     mingain=real(horner(h,exp(%i*ws*dt)))
57   end
58   
59   k=find(mingain<0)
60   if k==[] then gm=%inf,fr=[],return,end
61   mingain=abs(mingain(k));
62   ws=abs(ws(k))// select positive frequency
63
64   gm=-20*log(mingain)/log(10) //tranform into Db
65   [gm,k]=min(gm);ws=ws(k);//select the minimum
66
67   fr=ws/(2*%pi) //transform in Hz
68 endfunction