* Bug 15368 fixed: freson() failed on some continuous systems 30/19730/3
Samuel GOUGEON [Wed, 7 Feb 2018 19:47:24 +0000 (20:47 +0100)]
  http://bugzilla.scilab.org/15368

  test_run cacsd freson // passes

Change-Id: I744b533bb57cb3e3a1810ae1ac8ed45cf7467343

scilab/CHANGES.md
scilab/modules/cacsd/macros/freson.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_15368.tst [new file with mode: 0644]

index 1307626..5317a31 100644 (file)
@@ -534,6 +534,7 @@ Known issues
 * [#15354](http://bugzilla.scilab.org/show_bug.cgi?id=15354): `invr(%s^2)` returned 0 instead of `1/%s^2`. For any scalar number, polynomial or rational `a`, `coffg(a)` returned `0` instead of `1`. `coffg([])` yielded an error. `invr` and `coffg` had no unitary tests. The `coffg` help page was inaccurate and unclear.
 * [#15360](http://bugzilla.scilab.org/show_bug.cgi?id=15360): `numer()` and `denom()` were almost useless, unused, and with more handy replacements. They are declared obsolete to be removed in Scilab 6.1.0.
 * [#15363](http://bugzilla.scilab.org/show_bug.cgi?id=15363): `h2norm()` could no longer be applied to undefined time domain systems, and made the residu() right example failing.
+* [#15368](http://bugzilla.scilab.org/show_bug.cgi?id=15368): `freson()` wrongly returned [] (no peak frequency detected) for some continuous time linear systems.
 * [#15370](http://bugzilla.scilab.org/show_bug.cgi?id=15370): `bezout()` mishandled its output arguments.
 * [#15375](http://bugzilla.scilab.org/show_bug.cgi?id=15375): A .zcos file opened as a palette was greyed out.
 * [#15395](http://bugzilla.scilab.org/show_bug.cgi?id=15395): `ones(2,3,2) / %z` yielded an error..
index 1a15f4c..adff61a 100644 (file)
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function fr=freson(h)
-    [lhs,rhs]=argn(0);
+function fr = freson(h)
+    [lhs,rhs] = argn(0);
 
     if argn(2) < 1 then
-        error(msprintf(_("%s: Wrong number of input argument(s): %d expected.\n"),"freson",1));
+        msg = gettext("%s: Wrong number of input argument(s): %d expected.\n")
+        error(msprintf(msg, "freson", 1));
     end
     if and(typeof(h)<>["state-space","rational","zpk"]) then
-        ierr=execstr("[gm,fr]=%"+overloadname(sys)+"_freson(h)","errcatch")
+        ierr = execstr("[gm,fr]=%"+overloadname(sys)+"_freson(h)","errcatch")
         if ierr<>0 then
-            error(msprintf(gettext("%s: Wrong type for input argument: Linear dynamical system expected.\n"),"freson",1))
+            msg = gettext("%s: Wrong type for input argument: Linear dynamical system expected.\n")
+            error(msprintf(msg, "freson", 1))
         end
         return
     end
     if size(h,"*")<>1 then
-        error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"freson",1))
+        msg = gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n")
+        error(msprintf(msg, "freson", 1))
     end
 
     if typeof(h)=="state-space" then
@@ -32,35 +35,36 @@ function fr=freson(h)
     elseif typeof(h)=="zpk" then
         h=zpk2tf(h)
     end
-    dt=h.dt
-    [n,d]=h(["num","den"]);
+    dt = h.dt
+    [n,d] = h(["num","den"]);
     if type(n)==1 then
         n=poly(n,varn(d),"c");
     end
     if coeff(d,0)==0 then
-        error(msprintf(_("%s: Wrong value for input argument #%d: infinite gain at zero frequency.\n"),"freson",1))
+        msg = gettext("%s: Wrong value for input argument #%d: infinite gain at zero frequency.\n")
+        error(msprintf(msg, "freson", 1))
     end
 
     //look for  omega such that derivative of magn. is zero
     if dt=="c" then
         //find frequencies which zeros the magnitude derivative
-        hh=h*horner(h,-%s)
-        r=roots(derivat(hh).num)
-        k=find(imag(r)>0&abs(real(r))<%eps*abs(r));
-        fr=imag(r(k))/(2*%pi)
+        hh = h*horner(h,-%s)
+        r = roots(derivat(hh).num)
+        // k=find(imag(r)>0&abs(real(r))<%eps*abs(r)); // http://bugzilla.scilab.org/15368
+        k = find(imag(r)>0 & abs(real(r))<1.e-14*abs(r));
+        fr = imag(r(k))/(2*%pi)
     else
-        if dt=="d"|dt==[] then dt=1;end
+        if dt=="d" | dt==[] then dt = 1; end
         //find frequencies which zeros the magnitude derivative
-        hh=h*horner(h,1/%z)
-        r=roots(derivat(hh).num);
-        k=find(abs(abs(r)-1)<=sqrt(%eps)*abs(r));
-        r=imag(log(r(k)));
-        fr=r(r>0&r<0.999*%pi)/(2*%pi*dt)
+        hh = h*horner(h,1/%z)
+        r = roots(derivat(hh).num);
+        k = find(abs(abs(r)-1)<=sqrt(%eps)*abs(r));
+        r = imag(log(r(k)));
+        fr = r(r>0&r<0.999*%pi)/(2*%pi*dt)
     end
     if fr==[] then return;end
-    //find frequencies that correspond to a magnitude peak
-    k=find(abs(repfreq(h,fr))-abs(repfreq(h,fr*0.999))>0)
-    fr=fr(k)
-    if fr==[] then return;end
-    fr=gsort(fr,"g","d");
+
+    //find frequencies that correspond to a magnitude peak and sort them:
+    k = find(abs(repfreq(h,fr))-abs(repfreq(h,fr*0.999))>0)
+    fr = gsort(fr(k),"g","d");
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_15368.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_15368.tst
new file mode 100644 (file)
index 0000000..3cd4b65
--- /dev/null
@@ -0,0 +1,33 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2018 - Lucien POVY
+// Copyright (C) 2018 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+//
+// <-- Non-regression test for bug 15368 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/15368
+//
+// <-- Short Description -->
+// For some continuous time systems, freson() wrongly returned []
+
+s = %s;
+sl1 = syslin("c",1,s*(s+1)*(1+s/3))
+[K,P] = kpure(sl1)
+sl = K*sl1
+taud = 1/imag(P)
+sltaud=sl*(1+taud*s)//first open loop
+sltaud2=sl*(1+2*taud*s)//second open loop
+ff = freson(sltaud/(1+sltaud))//freson of first close loop
+assert_checkalmostequal(ff, 0.32447218051387, 1e-14);
+
+SLBF = sltaud2/(1+sltaud2); //THE CLOSE LOOP OF sltau2
+fr6 = freson(SLBF);
+assert_checkalmostequal(fr6, 0.48257117694410, 1e-14);
+