summaryrefslogtreecommitdiffstats
path: root/scilab/modules
diff options
context:
space:
mode:
authorSerge Steer <serge.steer@scilab.org>2011-04-18 14:34:08 +0200
committerVincent COUVERT <vincent.couvert@scilab.org>2011-05-06 10:29:42 +0200
commit7e9b83d4b50a6bbdab3fb7b740fd6d968d3140a2 (patch)
tree165af4a3b8997788abb6be24eaac675139057715 /scilab/modules
parent6346cbd2e5274deef40efa7c1472320e27db5ad1 (diff)
downloadscilab-7e9b83d4b50a6bbdab3fb7b740fd6d968d3140a2.zip
scilab-7e9b83d4b50a6bbdab3fb7b740fd6d968d3140a2.tar.gz
bug 6744 fixed (again);p_margin() returned an erroneous result
Change-Id: If9ce487a1639ae2b2a8ddf28af4353a1e2ec0695
Diffstat (limited to 'scilab/modules')
-rw-r--r--scilab/modules/cacsd/macros/p_margin.sci26
1 files changed, 16 insertions, 10 deletions
diff --git a/scilab/modules/cacsd/macros/p_margin.sci b/scilab/modules/cacsd/macros/p_margin.sci
index 439c970..8a06636 100644
--- a/scilab/modules/cacsd/macros/p_margin.sci
+++ b/scilab/modules/cacsd/macros/p_margin.sci
@@ -1,5 +1,5 @@
1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) 2010 - INRIA - Serge Steer 2// Copyright (C) 2010-2011 - INRIA - Serge Steer
3// 3//
4// This file must be used under the terms of the CeCILL. 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 5// This source file is licensed as described in the file COPYING, which
@@ -8,7 +8,7 @@
8// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt 8// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9 9
10function [phm,fr]=p_margin(h) 10function [phm,fr]=p_margin(h)
11 //compute the phase margin of a SISO transfer function 11//compute the phase margin of a SISO transfer function
12 select typeof(h) 12 select typeof(h)
13 case 'rational' then 13 case 'rational' then
14 case 'state-space' then 14 case 'state-space' then
@@ -17,7 +17,7 @@ function [phm,fr]=p_margin(h)
17 error(97,1) 17 error(97,1)
18 end 18 end
19 if or(size(h)<>[1 1]) then 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"),"p_margin",1)) 20 error(msprintf(_("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"p_margin",1))
21 end 21 end
22 eps=1.e-7;// threshold used for testing if complex numbers are real or pure imaginary 22 eps=1.e-7;// threshold used for testing if complex numbers are real or pure imaginary
23 23
@@ -29,16 +29,16 @@ function [phm,fr]=p_margin(h)
29 w=roots(real(niw*conj(niw)-diw*conj(diw)),"e"); 29 w=roots(real(niw*conj(niw)-diw*conj(diw)),"e");
30 //select positive real roots 30 //select positive real roots
31 ws=real(w(find((abs(imag(w))<eps)&(real(w)>0)))); //frequency points with unitary modulus 31 ws=real(w(find((abs(imag(w))<eps)&(real(w)>0)))); //frequency points with unitary modulus
32 if ws==[] then 32 if ws==[] then
33 phm=[]; 33 phm=[];
34 fr=[]; 34 fr=[];
35 return 35 return
36 end 36 end
37 f=horner(h,%i*ws); 37 f=horner(h,%i*ws);
38 else //discrete time case 38 else //discrete time case
39 if h.dt=='d' then 39 if h.dt=='d' then
40 dt=1; 40 dt=1;
41 else 41 else
42 dt=h.dt; 42 dt=h.dt;
43 end 43 end
44 // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt)) 44 // |h(e^(i*w*dt))|=1 <-- h(e^(i*w*dt))*h(e^(-i*w*dt))
@@ -52,15 +52,21 @@ function [phm,fr]=p_margin(h)
52 z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1 52 z(abs(abs(z)-1)>eps)=[];// retain only roots with modulus equal to 1
53 w=log(z)/(%i*dt); 53 w=log(z)/(%i*dt);
54 ws=real(w(abs(imag(w))<eps&real(w)>0)); //frequency points with unitary modulus 54 ws=real(w(abs(imag(w))<eps&real(w)>0)); //frequency points with unitary modulus
55 if ws==[] then 55 if ws==[] then
56 phm=%inf; 56 phm=%inf;
57 fr=[]; 57 fr=[];
58 return 58 return
59 end 59 end
60 f=horner(h,exp(%i*ws*dt)); 60 f=horner(h,exp(%i*ws*dt));
61 end 61 end
62 phm=atand(imag(f),real(f));// phase of the frequency response in [-180 180] 62 phi=atand(imag(f),real(f));// phase of the frequency response (in [-180 180])
63 phm=pmodulo(phm,360)-180; 63
64 [phm,k]=min(phm); 64 //avoid near 0 negative phases that will give phm=180 instead of -180
65 phi(phi>-1e-12&phi<0)=0;
66 //compute the margins
67 phm=pmodulo(phi,360)-180;
68 //select the min value together with associated frequency in Hz
69 [w,k]=min(abs(phm));
70 phm=phm(k)
65 fr=ws(k)/(2*%pi); 71 fr=ws(k)/(2*%pi);
66endfunction 72endfunction