1 files changed, 42 insertions, 0 deletions
diff --git a/modules/cacsd/macros/freson.sci b/modules/cacsd/macros/freson.sci
new file mode 100755
index 000000000..a27394a73
--- /dev/null
+++ b/modules/cacsd/macros/freson.sci
@@ -0,0 +1,42 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA -
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function fr=freson(h)
+    [lhs,rhs]=argn(0);
+    [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))
+    end
+    //look for  omega such that derivative of magn. is zero
+    niw=horner(n,%i*poly(0,"w"));
+    diw=horner(d,%i*poly(0,"w"));
+    niw=real(niw*conj(niw));diw=real(diw*conj(diw));
+    modul_d=derivat(niw/diw);w=roots(modul_d.num);
+
+    // get extreme points
+    k=find(imag(w)==0&real(w)>=0);
+    if k==[] then
+        fr=[];
+        g=[];
+        return
+    end
+    w=gsort(real(w(k)),"g","i");
+
+    //find maximums
+    wx=max(w)+0.5;
+    if horner(modul_d,wx)<0 then
+        w=w($:-2:1);
+    else
+        w=w($-1:-2:1);
+    end
+    fr=w/(2*%pi);
+endfunction