* Bug #12794 fixed - Cacsd: standardize calfrq.sci 56/12356/3
Samuel Gougeon [Tue, 27 Aug 2013 10:38:29 +0000 (12:38 +0200)]
The bug is actually invalid, so just standardize calfrq()'s code.

Change-Id: I41d1f47ef9e191b0ac68d4b7cf416a0c31f5811f

scilab/CHANGES_5.5.X
scilab/modules/cacsd/macros/calfrq.sci

index 9793db7..e800e32 100644 (file)
@@ -544,7 +544,9 @@ Bug fixes
 
 * Bug #12790 fixed - Links to zcos files in doc were broken.
 
-* Bug #12795 fixed - Typo fixes in Cacsd doc.
+* Bug #12794 fixed - calfrq.sci code did not follow Scilab standard.
+
+* Bug #12795 fixed - Typo fixes in CACSD doc.
 
 * Bug #12800 fixed - Typo fix in Polynomials doc.
 
index bc845d5..524bf49 100644 (file)
 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 
 
-function [frq,bnds,splitf]=calfrq(h,fmin,fmax)
+function [frq, bnds, splitf] = calfrq(h, fmin, fmax)
     //!
 
-    eps=1.d-14   //minimum absolute lower frequency
-    k=0.001;     // Minimum relative prediction error in the nyquist plan
-    epss=0.002   // minimum frequency distance with a singularity
-    nptmax=5000  //maximum number of discretization points
-    tol=0.01     // Tolerance for testing pure imaginary numbers
+    eps = 1.d-14   //minimum absolute lower frequency
+    k = 0.001      // Minimum relative prediction error in the nyquist plan
+    epss = 0.002   // minimum frequency distance with a singularity
+    nptmax = 5000  //maximum number of discretization points
+    tol = 0.01     // Tolerance for testing pure imaginary numbers
 
     // Check inputs
     // ------------
-    if and(typeof(h)<>["state-space" "rational"])
-        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"calfrq",1))
+    if and(typeof(h) <> ["state-space" "rational"])
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"), "calfrq", 1))
     end
-    if typeof(h)=="state-space" then
-        h=ss2tf(h)
+    if typeof(h) == "state-space" then
+        h = ss2tf(h)
     end
 
-    [m,n]=size(h.num)
-    dom=h("dt")
+    [m, n] = size(h.num)
+    dom = h("dt")
     select dom
     case "d" then
-        dom=1
+        dom = 1
     case [] then
-        error(96,1)
+        error(96, 1)
     case 0 then
-        error(96,1)
+        error(96, 1)
     end;
 
-    if type(dom)==1 then
-        nyq_frq=1/2/dom;
-        if fmax>nyq_frq then
-            warning(msprintf(gettext("%s: Frequencies beyond Nyquist frequency are ignored.\n"),"calfrq"));
-            fmax=min(fmax,nyq_frq)
+    if type(dom) == 1 then
+        nyq_frq = 1/2/dom;
+        if fmax > nyq_frq then
+            warning(msprintf(gettext("%s: Frequencies beyond Nyquist frequency are ignored.\n"), "calfrq"));
+            fmax = min(fmax, nyq_frq)
         end
-        if fmin<-nyq_frq then
-            warning(msprintf(gettext("%s: Negative frequencies below Nyquist frequency are ignored.\n"),"calfrq"));
-            fmin=max(fmin,-nyq_frq)
+        if fmin < -nyq_frq then
+            warning(msprintf(gettext("%s: Negative frequencies below Nyquist frequency are ignored.\n"), "calfrq"));
+            fmin = max(fmin, -nyq_frq)
         end
     end
     // Use symmetry to reduce the range
     // --------------------------------
-    if fmin<0&fmax>=0 then
-        [frq,bnds,splitf]=calfrq(h,eps,-fmin)
-        ns1=size(splitf,"*")-1;
-        nsp=size(frq,"*");
-        bnds=[bnds(1),bnds(2),-bnds(4),-bnds(3)];
+    if fmin < 0 & fmax >= 0 then
+        [frq, bnds, splitf] = calfrq(h, eps, -fmin)
+        ns1 = size(splitf, "*")-1;
+        nsp = size(frq, "*");
+        bnds = [bnds(1), bnds(2), -bnds(4), -bnds(3)];
 
-        if fmax>eps then
-            if fmax==-fmin then
-                splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2),nsp*ones(ns1)+splitf(2:$)];
-                bnds=[bnds(1),bnds(2),min(bnds(3),-bnds(3)),max(bnds(4),-bnds(4))];
-                frq=[-frq($:-1:1),frq]
+        if fmax > eps then
+            if fmax == -fmin then
+                splitf = [1, (nsp+2)*ones(1,ns1)-splitf($:-1:2), nsp*ones(ns1)+splitf(2:$)];
+                bnds = [bnds(1), bnds(2), min(bnds(3), -bnds(3)), max(bnds(4), -bnds(4))];
+                frq = [-frq($:-1:1), frq]
             else
-                [frq2,bnds2,splitf2]=calfrq(h,eps,fmax);
-                ns2=size(splitf2,"*")-1
-                splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2),nsp*ones(ns2)+splitf2(2:$)];
-                bnds=[min(bnds(1),bnds2(1)),max(bnds(2),bnds2(2)),...
-                min(bnds(3),bnds2(3)),max(bnds(4),bnds2(4))];
-                frq=[-frq($:-1:1),frq2]
+                [frq2, bnds2, splitf2] = calfrq(h, eps, fmax);
+                ns2 = size(splitf2,"*")-1
+                splitf = [1, (nsp+2)*ones(1,ns1)-splitf($:-1:2), nsp*ones(ns2)+splitf2(2:$)];
+                bnds = [min(bnds(1), bnds2(1)), max(bnds(2), bnds2(2)),...
+                min(bnds(3), bnds2(3)), max(bnds(4), bnds2(4))];
+                frq = [-frq($:-1:1), frq2]
             end
             return
         else
-            frq=-frq($:-1:1);
-            nsp=size(frq,"*");
+            frq = -frq($:-1:1);
+            nsp = size(frq, "*");
 
-            splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2)]
-            bnds=bnds;
+            splitf = [1, (nsp+2)*ones(1, ns1)-splitf($:-1:2)]
+            bnds = bnds;
             return;
         end
-    elseif fmin<0&fmax<=0 then
-        [frq,bnds,splitf]=calfrq(h,-fmax,-fmin)
-        ns1=size(splitf,"*")-1;
-        frq=-frq($:-1:1);
-        nsp=size(frq,"*");
-        splitf=[1, (nsp+2)*ones(1,ns1)-splitf($:-1:2)]
-        bnds=[bnds(1),bnds(2),-bnds(4),-bnds(3)];
+    elseif fmin < 0 & fmax <= 0 then
+        [frq, bnds, splitf] = calfrq(h, -fmax, -fmin)
+        ns1 = size(splitf, "*")-1;
+        frq = -frq($:-1:1);
+        nsp = size(frq, "*");
+        splitf = [1, (nsp+2)*ones(1, ns1)-splitf($:-1:2)]
+        bnds = [bnds(1), bnds(2), -bnds(4), -bnds(3)];
         return;
     elseif fmin >= fmax then
         error(msprintf(gettext("%s: Wrong value for input arguments #%d and #%d: %s < %s expected.\n"),..
-        "calfrq",2,3,"fmin","fmax"));
+        "calfrq", 2, 3, "fmin", "fmax"));
     end
 
     // Compute dicretisation over a given range
     // ----------------------------------------
 
 
-    splitf=[]
-    if fmin==0 then fmin=min(1d-14,fmax/10);end
+    splitf = []
+    if fmin == 0 then fmin = min(1d-14, fmax/10);end
     //
-    denh=h("den");numh=h("num")
-    l10=log(10)
+    denh = h("den"); numh = h("num")
+    l10 = log(10)
 
     // Locate singularities to avoid them
     // ----------------------------------
-    if dom=="c" then
-        c=2*%pi;
-        //selection function for singularities in the frequency range
-        deff("f=%sel(r,fmin,fmax,tol)",["f=[],";
-        "if prod(size(r))==0 then return,end";
-        "f=imag(r(find((abs(real(r))<=tol*abs(r))&(imag(r)>=0))))";
-        "if f<>[] then  f=f(find((f>fmin-tol)&(f<fmax+tol)));end"]);
+    if dom == "c" then
+        c = 2*%pi;
+        // selection function for singularities in the frequency range
+        deff("f=%sel(r, fmin, fmax, tol)",["f = [],";
+        "if prod(size(r)) == 0 then return, end";
+        "f = imag(r(find((abs(real(r))<=tol*abs(r))&(imag(r)>=0))))";
+        "if f <> [] then f = f(find((f>fmin-tol)&(f<fmax+tol))); end"]);
     else
-        c=2*%pi*dom
-        //selection function for singularities in the frequency range
-        deff("[f]=%sel(r,fmin,fmax,dom,tol)",["f=[],";
-        "if prod(size(r))==0 then return,end";
-        "f=r(find( ((abs(abs(r)-ones(r)))<=tol)&(imag(r)>=0)))";
-        "if f<>[] then ";
-        "  f=atan(imag(f),real(f));nf=prod(size(f))";
-        "  for k=1:nf ,";
-        "    kk=int((fmax-f(k))/(2*%pi))+1;";
-        "    f=[f;f(1:nf)+2*%pi*kk*ones(nf,1)];";
+        c = 2*%pi*dom
+        // selection function for singularities in the frequency range
+        deff("[f] = %sel(r, fmin, fmax, dom, tol)",["f = [],";
+        "if prod(size(r)) == 0 then return, end";
+        "f = r(find( ((abs(abs(r)-ones(r)))<=tol)&(imag(r)>=0)))";
+        "if f <> [] then ";
+        "  f = atan(imag(f), real(f)); nf = prod(size(f))";
+        "  for k=1:nf,";
+        "    kk = int((fmax-f(k))/(2*%pi))+1;";
+        "    f = [f; f(1:nf)+2*%pi*kk*ones(nf, 1)];";
         "  end;"
-        "  f=f(find((f>fmin-tol)&(f<fmax+tol)))";
+        "  f = f(find((f>fmin-tol)&(f<fmax+tol)))";
         "end"]);
     end
 
-    sing=[];zers=[];
-    fmin=c*fmin,fmax=c*fmax;
+    sing = [];zers = [];
+    fmin = c*fmin, fmax = c*fmax;
 
     for i=1:m
-        sing=[sing;%sel(roots(denh(i),"e"),fmin,fmax,tol)];
+        sing = [sing; %sel(roots(denh(i), "e"), fmin, fmax, tol)];
     end
 
-    pp=gsort(sing',"g","i");npp=size(pp,"*");//'
+    pp = gsort(sing', "g", "i"); npp = size(pp, "*");//'
 
     // singularities just on the left of the range
-    kinf=find(pp<fmin)
-    if kinf<>[] then
-        fmin=fmin+tol
-        pp(kinf)=[]
+    kinf = find(pp<fmin)
+    if kinf <> [] then
+        fmin = fmin+tol
+        pp(kinf) = []
     end
 
     // singularities just on the right of the range
-    ksup=find(pp>=fmax)
-    if ksup<>[] then
-        fmax=fmax-tol
-        pp(ksup)=[]
+    ksup = find(pp>=fmax)
+    if ksup <> [] then
+        fmax = fmax-tol
+        pp(ksup) = []
     end
 
-    //check for nearly multiple singularities
-    if pp<>[] then
-        dpp=pp(2:$)-pp(1:$-1)
-        keq=find(abs(dpp)<2*epss)
-        if keq<>[] then pp(keq)=[],end
+    // check for nearly multiple singularities
+    if pp <> [] then
+        dpp = pp(2:$)-pp(1:$-1)
+        keq = find(abs(dpp)<2*epss)
+        if keq <> [] then pp(keq) = [], end
     end
 
-    if pp<>[] then
-        frqs=[fmin real(matrix([(1-epss)*pp;(1+epss)*pp],2*size(pp,"*"),1)') fmax]
+    if pp <> [] then
+        frqs = [fmin real(matrix([(1-epss)*pp; (1+epss)*pp], 2*size(pp, "*"), 1)') fmax]
         //'
     else
-        frqs=[fmin fmax]
+        frqs = [fmin fmax]
     end
-    nfrq=size(frqs,"*");
+    nfrq = size(frqs, "*");
 
     // Evaluate bounds of nyquist plot
     //-------------------------------
-    xt=[];Pas=[]
+    xt = []; Pas = []
     for i=1:2:nfrq-1
-        w=logspace(log(frqs(i))/log(10),log(frqs(i+1))/log(10),100);
-        xt=[xt,w]
-        Pas=[Pas w(2)-w(1)]
+        w = logspace(log(frqs(i))/log(10), log(frqs(i+1))/log(10), 100);
+        xt = [xt, w]
+        Pas = [Pas w(2)-w(1)]
     end
-    if dom=="c" then
-        rf=freq(h("num"),h("den"),%i*xt);
+    if dom == "c" then
+        rf = freq(h("num"), h("den"), %i*xt);
     else
-        rf=freq(h("num"),h("den"),exp(%i*xt));
+        rf = freq(h("num"), h("den"), exp(%i*xt));
     end
     //
-    xmin=min(real(rf));xmax=max(real(rf));
-    ymin=min(imag(rf));ymax=max(imag(rf));
-    bnds=[xmin xmax ymin ymax];
-    dx=max([xmax-xmin,1]);dy=max([ymax-ymin,1]);
+    xmin = min(real(rf)); xmax = max(real(rf));
+    ymin = min(imag(rf)); ymax = max(imag(rf));
+    bnds = [xmin xmax ymin ymax];
+    dx = max([xmax-xmin, 1]); dy = max([ymax-ymin, 1]);
 
     // Compute discretization with a step adaptation method
     // ----------------------------------------------------
-    frq=[];
-    i=1;
-    nptr=nptmax; // number of unused discretization points
-    l10last=log10(frqs($));
+    frq = [];
+    i = 1;
+    nptr = nptmax; // number of unused discretization points
+    l10last = log10(frqs($));
     while i<nfrq
-        f0=frqs(i);fmax=frqs(i+1);
+        f0 = frqs(i); fmax = frqs(i+1);
         while f0==fmax do
-            i=i+2;
-            f=frqs(i);fmax=frqs(i+1);
+            i = i+2;
+            f = frqs(i); fmax = frqs(i+1);
         end
-        frq=[frq,f0];
-        pas=Pas(floor(i/2)+1)
-        splitf=[splitf size(frq,"*")];
+        frq = [frq, f0];
+        pas = Pas(floor(i/2)+1)
+        splitf = [splitf size(frq, "*")];
 
-        f=min(f0+pas,fmax);
+        f = min(f0+pas, fmax);
 
-        if dom=="c" then //cas continu
+        if dom == "c" then // Continuous case
             while f0<fmax
-                rf0=freq(h("num"),h("den"),(%i*f0));
-                rfc=freq(h("num"),h("den"),%i*f);
+                rf0 = freq(h("num"), h("den"), (%i*f0));
+                rfc = freq(h("num"), h("den"), %i*f);
                 // compute prediction error
-                epsd=pas/100;//epsd=1.d-8
+                epsd = pas/100; // epsd = 1.d-8
 
-                rfd=(freq(h("num"),h("den"),%i*(f0+epsd))-rf0)/(epsd);
-                rfp=rf0+pas*rfd;
+                rfd = (freq(h("num"), h("den"), %i*(f0+epsd))-rf0)/(epsd);
+                rfp = rf0+pas*rfd;
 
-                e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
-                if (e>k) then rf0=freq(h("num"),h("den"),(%i*f0));
-                    rfc=freq(h("num"),h("den"),%i*f);
+                e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
+                if e > k then rf0 = freq(h("num"), h("den"), (%i*f0));
+                    rfc = freq(h("num"), h("den"), %i*f);
                     // compute prediction error
-                    epsd=pas/100;//epsd=1.d-8
+                    epsd = pas/100; // epsd = 1.d-8
 
-                    rfd=(freq(h("num"),h("den"),%i*(f0+epsd))-rf0)/(epsd);
-                    rfp=rf0+pas*rfd;
+                    rfd = (freq(h("num"), h("den"), %i*(f0+epsd))-rf0)/(epsd);
+                    rfp = rf0+pas*rfd;
 
-                    e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
+                    e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
                     // compute minimum frequency logarithmic step to ensure a maximum
                     //of nptmax points to discretize
-                    pasmin=f0*(10^((l10last-log10(f0))/(nptr+1))-1)
-                    pas=pas/2
-                    if pas<pasmin then
-                        pas=pasmin
-                        frq=[frq,f];nptr=max([1,nptr-1])
-                        f0=f;f=min(f0+pas,fmax)
+                    pasmin = f0*(10^((l10last-log10(f0))/(nptr+1))-1)
+                    pas = pas/2
+                    if pas < pasmin then
+                        pas = pasmin
+                        frq = [frq, f]; nptr = max([1, nptr-1])
+                        f0 = f; f = min(f0+pas, fmax)
                     else
-                        f=min(f0+pas,fmax)
+                        f = min(f0+pas, fmax)
                     end
-                elseif e<k/2 then
-                    pas=2*pas
-                    frq=[frq,f];nptr=max([1,nptr-1])
-                    f0=f;f=min(f0+pas,fmax),
+                elseif e < k/2 then
+                    pas = 2*pas
+                    frq = [frq, f]; nptr = max([1, nptr-1])
+                    f0 = f; f = min(f0+pas, fmax),
                 else
-                    frq=[frq,f];nptr=max([1,nptr-1])
-                    f0=f;f=min(f0+pas,fmax),
+                    frq = [frq, f];nptr = max([1, nptr-1])
+                    f0 = f; f = min(f0+pas, fmax),
                 end
             end
-        else  //cas discret
-            pas=pas/dom
+        else  // Discrete case
+            pas = pas/dom
             while f0<fmax
-                rf0=freq(h("num"),h("den"),exp(%i*f0))
-                rfd=dom*(freq(h("num"),h("den"),exp(%i*(f0+pas/100)))-rf0)/(pas/100);
-                rfp=rf0+pas*rfd
-                rfc=freq(h("num"),h("den"),exp(%i*f));
-                e=max([abs(imag(rfp-rfc))/dy;abs(real(rfp-rfc))/dx])
-                if (e>k) then
-                    pasmin=f0*(10^((l10last-log10(f0))/(nptr+1))-1)
-                    pas=pas/2
-                    if pas<pasmin then
-                        pas=pasmin
-                        frq=[frq,f];nptr=max([1,nptr-1])
-                        f0=f;f=min(f0+pas,fmax)
+                rf0 = freq(h("num"), h("den"), exp(%i*f0))
+                rfd = dom*(freq(h("num"), h("den"), exp(%i*(f0+pas/100)))-rf0)/(pas/100);
+                rfp = rf0+pas*rfd
+                rfc = freq(h("num"), h("den"), exp(%i*f));
+                e = max([abs(imag(rfp-rfc))/dy; abs(real(rfp-rfc))/dx])
+                if e > k then
+                    pasmin = f0*(10^((l10last-log10(f0))/(nptr+1))-1)
+                    pas = pas/2
+                    if pas < pasmin then
+                        pas = pasmin
+                        frq = [frq, f]; nptr = max([1, nptr-1])
+                        f0 = f; f = min(f0+pas, fmax)
                     else
-                        f=min(f0+pas,fmax)
+                        f = min(f0+pas, fmax)
                     end
-                elseif e<k/2 then
-                    pas=2*pas
-                    frq=[frq,f];nptr=max([1,nptr-1])
-                    f0=f;f=min(f0+pas,fmax),
+                elseif e < k/2 then
+                    pas = 2*pas
+                    frq = [frq, f]; nptr = max([1, nptr-1])
+                    f0 = f; f = min(f0+pas, fmax),
                 else
-                    frq=[frq,f];nptr=max([1,nptr-1])
-                    f0=f;f=min(f0+pas,fmax),
+                    frq = [frq, f]; nptr = max([1, nptr-1])
+                    f0 = f; f = min(f0+pas, fmax),
                 end
             end
         end
-        i=i+2
+        i = i+2
     end
-    frq( size(frq,"*") )=fmax
-    frq=frq/c;
+    frq(size(frq, "*")) = fmax
+    frq = frq/c;
 endfunction