From 818cbe80162eb20b7f030efbf2ad11bef7ffa4e0 Mon Sep 17 00:00:00 2001 From: Samuel Gougeon Date: Tue, 27 Aug 2013 12:38:29 +0200 Subject: [PATCH] * Bug #12794 fixed - Cacsd: standardize calfrq.sci The bug is actually invalid, so just standardize calfrq()'s code. Change-Id: I41d1f47ef9e191b0ac68d4b7cf416a0c31f5811f --- scilab/CHANGES_5.5.X | 4 +- scilab/modules/cacsd/macros/calfrq.sci | 332 ++++++++++++++++---------------- 2 files changed, 169 insertions(+), 167 deletions(-) diff --git a/scilab/CHANGES_5.5.X b/scilab/CHANGES_5.5.X index 9793db7..e800e32 100644 --- a/scilab/CHANGES_5.5.X +++ b/scilab/CHANGES_5.5.X @@ -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. diff --git a/scilab/modules/cacsd/macros/calfrq.sci b/scilab/modules/cacsd/macros/calfrq.sci index bc845d5..524bf49 100644 --- a/scilab/modules/cacsd/macros/calfrq.sci +++ b/scilab/modules/cacsd/macros/calfrq.sci @@ -8,269 +8,269 @@ // 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=0))))"; + "if f <> [] then f = f(find((f>fmin-tol)&(f=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)&(ffmin-tol)&(f[] then - fmin=fmin+tol - pp(kinf)=[] + kinf = find(pp [] 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 ik) 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 pask) then - pasmin=f0*(10^((l10last-log10(f0))/(nptr+1))-1) - pas=pas/2 - if pas 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