bug 7095 fix, the margins unit_test fails, g_margin, p_margin and unit_test improved 35/1635/8
Serge Steer [Thu, 19 Aug 2010 14:59:06 +0000 (16:59 +0200)]
Change-Id: I9cd32ba419a794877e100b0d4d1291fcc9c984b9

scilab/CHANGES_5.3.X
scilab/modules/cacsd/macros/g_margin.sci
scilab/modules/cacsd/macros/p_margin.sci
scilab/modules/cacsd/tests/unit_tests/margins.dia.ref
scilab/modules/cacsd/tests/unit_tests/margins.tst

index 9c8a110..3bfe7f3 100644 (file)
@@ -149,6 +149,8 @@ Bug Fixes:
 * bug 6767 fixed - The "uicontrol" help page said the "Style" property of an
                    uicontrol can be set after creation but it does not.
 
+* bug 7095 fixed - the margins unit_test failed.
+
 * bug 7163 fixed - tree2code wrongly replaced every (:) occurence in a function
                    definition with (eye()).
 
index a435d40..13d9bc0 100644 (file)
@@ -47,8 +47,9 @@ function [gm,fr]=g_margin(h)
   if k==[] then gm=%inf,fr=[],return,end
   mingain=abs(mingain(k));
   ws=abs(ws(k))// select positive frequency
-  //disp([ws,1/mingain])
+
   gm=-20*log(mingain)/log(10) //tranform into Db
-  [gm,k]=min(gm);ws=ws(k)//select the minimum
-  fr=ws/(2*%pi) 
+  [gm,k]=min(gm);ws=ws(k);//select the minimum
+
+  fr=ws/(2*%pi) //transform in Hz
 endfunction
index 2bfe18f..439c970 100644 (file)
@@ -60,7 +60,7 @@ function [phm,fr]=p_margin(h)
     f=horner(h,exp(%i*ws*dt));
   end
   phm=atand(imag(f),real(f));// phase of the frequency response in [-180 180]
-  phm=pmodulo(phm-180,360);
-  [phm,k]=min(phm)
+  phm=pmodulo(phm,360)-180; 
+  [phm,k]=min(phm);
   fr=ws(k)/(2*%pi);
 endfunction
index 5164349..5e4cc23 100644 (file)
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- ENGLISH IMPOSED -->
 //gain margin
 //-----------
+function ok=check_gmargin(h,fref,gref)
+  eps=1e-10
+  [g,f]=g_margin(h)
+  if fref==0 then
+     ok=abs(f - fref)<eps
+  else
+    ok=abs(f - fref)/abs(fref)<eps
+  end
+  if gref==0 then
+     ok=ok&abs(g - gref)<eps
+  else
+    ok=ok&abs(g - gref)/abs(gref)<eps
+  end
+endfunction
 //discrete time case
-h = syslin(0.1,0.04798*%z+0.0464,%z^2-1.81*%z+0.9048);//pk
-[g ,f]=g_margin(h);
-if norm(g - 6.2423793565916954)>1d-5  then bugmes();quit;end
-if norm(f - 0.86539118158628991)>1d-5  then bugmes();quit;end
-h=syslin(0.1,0.086-0.161*%z+0.078*%z^2,-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[g ,f]=g_margin(h);
-if norm(g- 12.3470513198133762)>1d-5  then bugmes();quit;end
-if norm(f-0.26109869168929700)>1d-5  then bugmes();quit;end
-h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[g ,f]=g_margin(h);
-if norm(g-2.80462622542011131)>1d-5  then bugmes();quit;end
-if norm(f-0.26109869168929828)>1d-5  then bugmes();quit;end
+z=poly(0,'z');
+//the references solutions are computed using the following Maple instructions
+//Digits:=50:
+//z:=exp(I*w):
+//assume(w, 'real', w > 0)
+//h:= rational fraction in z
+//M := solve(Im(h) = 0, w)
+//evalf(M/(2*Pi))
+//evalf(eval(20*log10(1/abs(h)), w = M[k]))
+h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2));
+f_ref=atan(sqrt(9003479136639)/4963519)/(0.1*2*%pi);
+g_ref=6.242379356591869534116;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
+f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
+g_ref=12.3470513198103944037870;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));  
+f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
+g_ref=2.8046262254171456578864840;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h=syslin(1,((21/500)+(-3933/100000)*z+(-15407/100000)*z^2+(9259/50000)*z^3+(6939/(10^12))*z^4)/((-21/500)+(11/125)*z+(57/500)*z^2+(-9/25)*z^3+(1/5)*z^4));  
+f_ref=1/2;
+g_ref=6.48227782514616029706190;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
 //continuuous time case
-h=syslin('c',-1+%s,3+2*%s+%s^2);//ok
-[g,f]=g_margin(h);
-if norm(g-20*log(3)/log(10))>1d-5  then bugmes();quit;end
-if norm(f)>1d-5  then bugmes();quit;end
- h = syslin('c',0.2+0.8*%s+0.3*%s^3,0.0409+0.1827*%s+1.28225*%s^2+3.1909*%s^3+2.56*%s^4+%s^5);
-[g,f]=g_margin(h);
-if norm(g + 4.91687406933815385)>1d-5  then bugmes();quit;end
-if norm(f-0.07145552582020068)>1d-5  then bugmes();quit;end
+s=poly(0,'s');
+//the reference solutions are computed using  the following Maple instructions
+//s:=I*w:
+//assume(w, 'real', w > 0)
+//h:= rational fraction in z
+//M := solve(Im(h) = 0, w)
+//evalf(M/(2*Pi))
+//evalf(eval(20*log10(1/abs(h)), w = M[k]))
+h=syslin('c',(-1+s)/(3+2*s+s^2));  
+f_ref=0;
+g_ref=9.5424250943932487459005580;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h = syslin('c',((1/5)+(4/5)*s+(0/1)*s^2+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4 +s^5)); 
+f_ref=0.0714555258202006740373134;
+g_ref=-4.91687406933815400242335;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h=syslin('c',485000/(10000*s+200*s^2+s^3));  
+f_ref=100/(2*%pi);
+g_ref=12.305765141234350772862319;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
+h = syslin('c',1/(s+2*s^2+s^3));
+f_ref=1/(2*%pi);
+g_ref=6.0205999132796239042747779;
+if ~check_gmargin(h,f_ref,g_ref) then bugmes();quit;end
 //phase margin
 //-----------
+function ok=check_pmargin(h,f_ref,p_ref)
+  eps=1e-9
+  [p,f]=p_margin(h)
+  if f_ref==0 then
+     ok=abs(f - f_ref)<eps
+  else
+    ok=abs(f - f_ref)/abs(f_ref)<eps
+  end
+  if p_ref==0 then
+     ok=ok&abs(p - p_ref)<eps
+  else
+    ok=ok&abs(p - p_ref)/abs(p_ref)<eps
+  end
+endfunction
 //discrete time case
-h = syslin(0.1,0.04798*%z+0.0464,%z^2-1.81*%z+0.9048);//ok
-[p ,f]=p_margin(h);
-if norm(p-13.5711556361199754)>1d-5  then bugmes();quit;end
-if norm(f-0.69301660031536083)>1d-5  then bugmes();quit;end
-h=syslin(0.1,0.086-0.161*%z+0.078*%z^2,-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
+//the reference solutions are computed using  the following Maple instructions
+//z:=exp(I*w):
+//assume(w, 'real', w > 0)
+//P := solve(abs(h) = 1, w)
+//evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
+h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2)); 
+f_ref=0.693016600315284442350578876;
+p_ref=13.57115563612946355428439468;
+if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
+h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
 [p ,f]=p_margin(h);
 if p<>%inf  then bugmes();quit;end
 if f<>[]  then bugmes();quit;end
-h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[p ,f]=p_margin(h);
-if norm(p-52.9496741629204308 )>1d-5  then bugmes();quit;end
-if norm(f-0.21233648894721618)>1d-5  then bugmes();quit;end
+h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);
+h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
+f_ref=0.212336488950669705771059018;
+p_ref=52.94967415965772478856630911;
+if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
 //continuous case
-h=syslin('c',1.1+2.4*%s+0.7*%s^2,3+2*%s+%s^2);
-[p,f]=p_margin(h);
-if norm(p-(-148.547076202317413))>1d-5  then bugmes();quit;end
-if norm(f-0.18945852512301298)>1d-5  then bugmes();quit;end
- h = syslin('c',0.2+0.8*%s+0.3*%s^3,0.0409+0.1827*%s+1.28225*%s^2+3.1909*%s^3+2.56*%s^4+%s^5);
-[p,f]=p_margin(h);
-if norm(p + 13.1128497150069734 )>1d-5  then bugmes();quit;end
-if norm(f-0.09144216563554157  )>1d-5  then bugmes();quit;end
+//the reference solutions are computed using  the following Maple instructions
+//z:=I*w:
+//assume(w, 'real', w > 0)
+//P := solve(abs(h) = 1, w)
+//evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
+h=syslin('c',((11/10)+(12/5)*s+(7/10)*s^2)/(3+2*s+s^2));
+f_ref=(1/51)*sqrt(15861-204*sqrt(3562))/(2*%pi);
+p_ref=-148.547076202317410601324666;
+if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
+h = syslin('c',((1/5)+(4/5)*s+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4+s^5));  
+f_ref=0.09144216563554157543991;
+p_ref=-13.1128497150069802772313;
+if  ~check_pmargin(h,f_ref,p_ref) then bugmes();quit;end
index e937db3..286120c 100644 (file)
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- ENGLISH IMPOSED -->
+
 //gain margin
 //-----------
+function ok=check_gmargin(h,fref,gref)
+  eps=1e-10
+  [g,f]=g_margin(h)
+  if fref==0 then
+     ok=abs(f - fref)<eps
+  else
+    ok=abs(f - fref)/abs(fref)<eps
+  end
+  if gref==0 then
+     ok=ok&abs(g - gref)<eps
+  else
+    ok=ok&abs(g - gref)/abs(gref)<eps
+  end
+endfunction
+
 //discrete time case
-h = syslin(0.1,0.04798*%z+0.0464,%z^2-1.81*%z+0.9048);//pk
-[g ,f]=g_margin(h);
-if norm(g - 6.2423793565916954)>1d-5  then pause,end
-if norm(f - 0.86539118158628991)>1d-5  then pause,end
+z=poly(0,'z');
+//the references solutions are computed using the following Maple instructions
+//Digits:=50:
+//z:=exp(I*w):
+//assume(w, 'real', w > 0)
+//h:= rational fraction in z
+//M := solve(Im(h) = 0, w)
+//evalf(M/(2*Pi))
+//evalf(eval(20*log10(1/abs(h)), w = M[k]))
+
+h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2));
+f_ref=atan(sqrt(9003479136639)/4963519)/(0.1*2*%pi);
+g_ref=6.242379356591869534116;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
+
 
 
-h=syslin(0.1,0.086-0.161*%z+0.078*%z^2,-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[g ,f]=g_margin(h);
-if norm(g- 12.3470513198133762)>1d-5  then pause,end
-if norm(f-0.26109869168929700)>1d-5  then pause,end
+h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));
+f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
+g_ref=12.3470513198103944037870;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
 
-h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[g ,f]=g_margin(h);
-if norm(g-2.80462622542011131)>1d-5  then pause,end
-if norm(f-0.26109869168929828)>1d-5  then pause,end
+h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3));  
+f_ref=atan((1/215)*sqrt(1136805-8330*sqrt(5970))/(833/43+(1/215)*sqrt(5970)))/(0.1*2*%pi);
+g_ref=2.8046262254171456578864840;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
 
 
+h=syslin(1,((21/500)+(-3933/100000)*z+(-15407/100000)*z^2+(9259/50000)*z^3+(6939/(10^12))*z^4)/((-21/500)+(11/125)*z+(57/500)*z^2+(-9/25)*z^3+(1/5)*z^4));  
+f_ref=1/2;
+g_ref=6.48227782514616029706190;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
+
 //continuuous time case
-h=syslin('c',-1+%s,3+2*%s+%s^2);//ok
-[g,f]=g_margin(h);
-if norm(g-20*log(3)/log(10))>1d-5  then pause,end
-if norm(f)>1d-5  then pause,end
+s=poly(0,'s');
+//the reference solutions are computed using  the following Maple instructions
+//s:=I*w:
+//assume(w, 'real', w > 0)
+//h:= rational fraction in z
+//M := solve(Im(h) = 0, w)
+//evalf(M/(2*Pi))
+//evalf(eval(20*log10(1/abs(h)), w = M[k]))
+
+h=syslin('c',(-1+s)/(3+2*s+s^2));  
+f_ref=0;
+g_ref=9.5424250943932487459005580;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
+
+
+h = syslin('c',((1/5)+(4/5)*s+(0/1)*s^2+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4 +s^5)); 
+f_ref=0.0714555258202006740373134;
+g_ref=-4.91687406933815400242335;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
+
 
+h=syslin('c',485000/(10000*s+200*s^2+s^3));  
+f_ref=100/(2*%pi);
+g_ref=12.305765141234350772862319;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
 
- h = syslin('c',0.2+0.8*%s+0.3*%s^3,0.0409+0.1827*%s+1.28225*%s^2+3.1909*%s^3+2.56*%s^4+%s^5);
-[g,f]=g_margin(h);
-if norm(g + 4.91687406933815385)>1d-5  then pause,end
-if norm(f-0.07145552582020068)>1d-5  then pause,end
+
+h = syslin('c',1/(s+2*s^2+s^3));
+f_ref=1/(2*%pi);
+g_ref=6.0205999132796239042747779;
+if ~check_gmargin(h,f_ref,g_ref) then pause,end
 
 //phase margin
 //-----------
+function ok=check_pmargin(h,f_ref,p_ref)
+  eps=1e-9
+  [p,f]=p_margin(h)
+  if f_ref==0 then
+     ok=abs(f - f_ref)<eps
+  else
+    ok=abs(f - f_ref)/abs(f_ref)<eps
+  end
+  if p_ref==0 then
+     ok=ok&abs(p - p_ref)<eps
+  else
+    ok=ok&abs(p - p_ref)/abs(p_ref)<eps
+  end
+endfunction
+
 //discrete time case
-h = syslin(0.1,0.04798*%z+0.0464,%z^2-1.81*%z+0.9048);//ok
-[p ,f]=p_margin(h);
-if norm(p-13.5711556361199754)>1d-5  then pause,end
-if norm(f-0.69301660031536083)>1d-5  then pause,end
+//the reference solutions are computed using  the following Maple instructions
+//z:=exp(I*w):
+//assume(w, 'real', w > 0)
+//P := solve(abs(h) = 1, w)
+//evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
 
+h = syslin(0.1,((29/625)+(2399/50000)*z)/((1131/1250)+(-181/100)*z+z^2)); 
+f_ref=0.693016600315284442350578876;
+p_ref=13.57115563612946355428439468;
+if  ~check_pmargin(h,f_ref,p_ref) then pause,end
 
-h=syslin(0.1,0.086-0.161*%z+0.078*%z^2,-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
+
+h=syslin(0.1,((43/500)+(-161/1000)*z+(39/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
 [p ,f]=p_margin(h);
 if p<>%inf  then pause,end
 if f<>[]  then pause,end
 
-h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);//ok
-[p ,f]=p_margin(h);
-if norm(p-52.9496741629204308 )>1d-5  then pause,end
-if norm(f-0.21233648894721618)>1d-5  then pause,end
-
+h=syslin(0.1,3*(0.086-0.161*%z+0.078*%z^2),-0.67+2.286*%z-2.61*%z^2+%z^3);
+h=syslin(0.1,((129/500)+(-483/1000)*z+(117/500)*z^2)/((-67/100)+(1143/500)*z+(-261/100)*z^2+z^3)); 
+f_ref=0.212336488950669705771059018;
+p_ref=52.94967415965772478856630911;
+if  ~check_pmargin(h,f_ref,p_ref) then pause,end
 
 //continuous case
-h=syslin('c',1.1+2.4*%s+0.7*%s^2,3+2*%s+%s^2);
-[p,f]=p_margin(h);
-if norm(p-(-148.547076202317413))>1d-5  then pause,end
-if norm(f-0.18945852512301298)>1d-5  then pause,end
-
- h = syslin('c',0.2+0.8*%s+0.3*%s^3,0.0409+0.1827*%s+1.28225*%s^2+3.1909*%s^3+2.56*%s^4+%s^5);
-[p,f]=p_margin(h);
-if norm(p + 13.1128497150069734 )>1d-5  then pause,end
-if norm(f-0.09144216563554157  )>1d-5  then pause,end
+//the reference solutions are computed using  the following Maple instructions
+//z:=I*w:
+//assume(w, 'real', w > 0)
+//P := solve(abs(h) = 1, w)
+//evalf(-(eval(180-180*argument(h)/Pi, w = Re(P[k]))))+360
+
+h=syslin('c',((11/10)+(12/5)*s+(7/10)*s^2)/(3+2*s+s^2));
+f_ref=(1/51)*sqrt(15861-204*sqrt(3562))/(2*%pi);
+p_ref=-148.547076202317410601324666;
+if  ~check_pmargin(h,f_ref,p_ref) then pause,end
+
+
+h = syslin('c',((1/5)+(4/5)*s+(3/10)*s^3)/((409/10000)+(1827/10000)*s+(5129/4000)*s^2+(31909/10000)*s^3+(64/25)*s^4+s^5));  
+f_ref=0.09144216563554157543991;
+p_ref=-13.1128497150069802772313;
+if  ~check_pmargin(h,f_ref,p_ref) then pause,end