Mercurial > hg > octave-lyh
changeset 7125:f084ba47812b
[project @ 2007-11-08 02:29:23 by jwe]
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,36 @@ +2007-11-07 Muthiah Annamalai <muthuspost@gmail.com> + + * control/base/bode_bounds.m, control/base/dgram.m, + control/base/dlyap.m, control/base/freqchkw.m, + control/base/gram.m, control/base/place.m, + control/hinf/hinf_ctr.m, control/hinf/hinfsyn_chk.m, + control/hinf/hinfsyn_ric.m, control/system/is_sample.m, + control/system/is_signal_list.m, control/system/ss2tf.m, + control/system/sys2fir.m, control/system/sysgettsam.m, + control/system/sysgettype.m, control/system/sysreorder.m, + control/system/tf2sys.m, control/system/zp2tf.m, + control/util/axis2dlim.m, control/util/swap.m, + control/util/zgfmul.m, control/util/zgfslv.m, + control/util/zginit.m, control/util/zgreduce.m, + control/util/zgrownorm.m, control/util/zgscal.m, + control/util/zgsgiv.m, control/util/zgshsr.m, general/isa.m, + geometry/inpolygon.m, linear-algebra/housh.m, + miscellaneous/compare_versions.m, miscellaneous/inputname.m, + miscellaneous/run.m, quaternion/qconj.m, + quaternion/qcoordinate_plot.m, quaternion/qderiv.m, + quaternion/qderivmat.m, quaternion/qinv.m, quaternion/qmult.m, + quaternion/qtrans.m, quaternion/qtransvmat.m, signal/fractdiff.m, + signal/freqz_plot.m, signal/periodogram.m, signal/rectangle_lw.m, + signal/rectangle_sw.m, signal/sinc.m, signal/triangle_lw.m, + signal/triangle_sw.m, signal/yulewalker.m, sparse/colperm.m, + sparse/etreeplot.m, sparse/nonzeros.m, sparse/spalloc.m, + sparse/spones.m, sparse/spy.m, specfun/isprime.m, + statistics/distributions/empirical_cdf.m, + statistics/distributions/empirical_inv.m, + statistics/distributions/empirical_pdf.m, + statistics/models/logistic_regression_derivatives.m, + statistics/models/logistic_regression_likelihood.m: Check nargin. + 2007-11-07 David Bateman <dbateman@free.fr> * general/gradient.m: Correctly convert deltax and deltay scalar
--- a/scripts/control/base/bode_bounds.m +++ b/scripts/control/base/bode_bounds.m @@ -36,32 +36,51 @@ function [wmin, wmax] = bode_bounds (zer, pol, DIGITAL, tsam) + if (nargin != 4) + print_usage (); + endif + ## make sure zer,pol are row vectors - if(!isempty(pol)) pol = reshape(pol,1,length(pol)); endif - if(!isempty(zer)) zer = reshape(zer,1,length(zer)); endif + if (! isempty (pol)) + pol = reshape (pol, 1, length (pol)); + endif + if (! isempty (zer)) + zer = reshape (zer, 1, length (zer)); + endif ## check for natural frequencies away from omega = 0 if (DIGITAL) ## The 2nd conditions prevents log(0) in the next log command - iiz = find(abs(zer - 1) > norm(zer) * eps && abs(zer) > norm(zer) * eps); - iip = find(abs(pol - 1) > norm(pol) * eps && abs(pol) > norm(pol) * eps); + iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps); + iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps); ## avoid dividing empty matrices, it would work but looks nasty - if (!isempty(iiz)) czer = log(zer(iiz))/tsam; - else czer = []; endif + if (! isempty (iiz)) + czer = log (zer(iiz))/tsam; + else + czer = []; + endif - if (!isempty(iip)) cpol = log(pol(iip))/tsam; - else cpol = []; endif - + if (! isempty (iip)) + cpol = log (pol(iip))/tsam; + else + cpol = []; + endif else ## continuous - iip = find((abs(pol)) > (norm(pol) * eps)); - iiz = find((abs(zer)) > (norm(zer) * eps)); + iip = find (abs(pol) > norm(pol)*eps)); + iiz = find (abs(zer) > norm(zer)*eps)); - if(!isempty(zer)) czer = zer(iiz); - else czer = []; endif - if(!isempty(pol)) cpol = pol(iip); - else cpol = []; endif + if (! isempty (zer)) + czer = zer(iiz); + else + czer = []; + endif + if (! isempty (pol)) + cpol = pol(iip); + else + cpol = []; + endif endif if (isempty (iip) && isempty (iiz)) @@ -69,8 +88,8 @@ wmin = -1; wmax = 3; else - wmin = floor(log10(min(abs([cpol,czer])))); - wmax = ceil(log10(max(abs([cpol,czer])))); + wmin = floor (log10 (min (abs ([cpol, czer])))); + wmax = ceil (log10 (max (abs ([cpol, czer])))); endif ## expand to show the entirety of the "interesting" portion of the plot @@ -78,6 +97,8 @@ wmax++; ## run digital frequency all the way to pi - if (DIGITAL) wmax = log10(pi/tsam); endif + if (DIGITAL) + wmax = log10 (pi/tsam); + endif endfunction
--- a/scripts/control/base/dgram.m +++ b/scripts/control/base/dgram.m @@ -61,6 +61,10 @@ function m = dgram (a, b) + if (nargin != 2) + print_usage(); + endif + ## let dlyap do the error checking... m = dlyap (a, b*b');
--- a/scripts/control/base/dlyap.m +++ b/scripts/control/base/dlyap.m @@ -90,6 +90,10 @@ function x = dlyap (a, b) + if (nargin != 2) + print_usage (); + endif + if ((n = issquare (a)) == 0) warning ("dlyap: a must be square"); endif @@ -129,19 +133,19 @@ blksiz = 1; endif - Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n); + Ajj = kron (s(j:j1,j:j1), s) - eye (blksiz*n); - rhs = reshape (b (:, j:j1), blksiz*n, 1); + rhs = reshape (b (:,j:j1), blksiz*n, 1); if (j1 < n) - rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)'); + rhs2 = s*(x(:,(j1+1):n) * s(j:j1,(j1+1):n)'); rhs = rhs + reshape (rhs2, blksiz*n, 1); endif v = - Ajj\rhs; - x (:, j) = v (1:n); + x(:,j) = v (1:n); - if(blksiz == 2) + if (blksiz == 2) x (:, j1) = v ((n+1):blksiz*n); endif
--- a/scripts/control/base/freqchkw.m +++ b/scripts/control/base/freqchkw.m @@ -29,6 +29,10 @@ function USEW = freqchkw (w) + if (nargin != 1) + print_usage (); + endif + if (isempty (w)) USEW = 0; elseif (! isvector (w))
--- a/scripts/control/base/gram.m +++ b/scripts/control/base/gram.m @@ -29,6 +29,10 @@ function m = gram (a, b) + if (nargin != 2) + print_usage (); + endif + ## Let lyap do the error checking... m = lyap (a, b*b');
--- a/scripts/control/base/place.m +++ b/scripts/control/base/place.m @@ -45,78 +45,81 @@ function K = place (sys, P) + if (nargin != 2) + print_usage (); + endif + ## check arguments - if(!isstruct(sys)) - error("sys must be in system data structure format (see ss)"); + if(! isstruct (sys)) + error ("sys must be in system data structure format (see ss)"); endif - sys = sysupdate(sys,"ss"); # make sure it has state space form up to date - if(!is_controllable(sys)) - error("sys is not controllable."); - elseif( min(size(P)) != 1) - error("P must be a vector") + sys = sysupdate (sys, "ss"); # make sure it has state space form up to date + if (! is_controllable (sys)) + error ("sys is not controllable"); + elseif (min (size (P)) != 1) + error ("P must be a vector") else - P = reshape(P,length(P),1); # make P a column vector + P = P(:); # make P a column vector endif ## system must be purely continuous or discrete - is_digital(sys); - [n,nz,m,p] = sysdimensions(sys); + is_digital (sys); + [n, nz, m, p] = sysdimensions (sys); nx = n+nz; # already checked that it's not a mixed system. if(m != 1) - error(["sys has ", num2str(m)," inputs; need only 1"]); + error ("sys has %d inputs; need only 1", m); endif ## takes the A and B matrix from the system representation - [A,B]=sys2ss(sys); - sp = length(P); - if(nx == 0) - error("place: A matrix is empty (0x0)"); - elseif(nx != length(P)) - error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ... - "entries."]) + [A, B] = sys2ss (sys); + sp = length (P); + if (nx == 0) + error ("place: A matrix is empty (0x0)"); + elseif (nx != length (P)) + error ("A=(%dx%d), P has %d entries", nx, nx, length (P)) endif ## arguments appear to be compatible; let's give it a try! ## The second step is the calculation of the characteristic polynomial ofA - PC=poly(A); + PC = poly (A); ## Third step: Calculate the transformation matrix T that transforms the state ## equation in the controllable canonical form. ## first we must calculate the controllability matrix M: - M=B; - AA=A; + M = B; + AA = A; for n = 2:nx - M(:,n)=AA*B; - AA=AA*A; + M(:,n) = AA*B; + AA = AA*A; endfor ## second, construct the matrix W - PCO=PC(nx:-1:1); - PC1=PCO; # Matrix to shift and create W row by row + PCO = PC(nx:-1:1); + PC1 = PCO; # Matrix to shift and create W row by row for n = 1:nx W(n,:) = PC1; - PC1=[PCO(n+1:nx),zeros(1,n)]; + PC1 = [PCO(n+1:nx), zeros(1,n)]; endfor - T=M*W; + T = M*W; ## finaly the matrix K is calculated - PD = poly(P); # The desired characteristic polynomial + PD = poly (P); # The desired characteristic polynomial PD = PD(nx+1:-1:2); PC = PC(nx+1:-1:2); K = (PD-PC)/T; ## Check if the eigenvalues of (A-BK) are the same specified in P - Pcalc = eig(A-B*K); + Pcalc = eig (A-B*K); - Pcalc = sortcom(Pcalc); - P = sortcom(P); + Pcalc = sortcom (Pcalc); + P = sortcom (P); - if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1) - disp("Place: Pole placed at more than 10% relative error from specified"); + if (max ((abs(Pcalc)-abs(P))./abs(P) ) > 0.1) + warning ("place: Pole placed at more than 10% relative error from specified"); endif endfunction
--- a/scripts/control/hinf/hinf_ctr.m +++ b/scripts/control/hinf/hinf_ctr.m @@ -57,6 +57,10 @@ function K = hinf_ctr (dgs, F, H, Z, g) + if (nargin != 5) + print_usage (); + endif + nw = dgs.nw; nu = dgs.nu; nz = dgs.nz; @@ -78,7 +82,7 @@ nout = nz + ny; nin = nw + nu; - nstates = size(A, 1); + nstates = size (A, 1); F11 = F(1:(nw-ny),:); F12 = F((nw-ny+1):nw,:); @@ -94,33 +98,33 @@ D1122 = D11((nz-nu+1):nz,(nw-ny+1):nw); ## D11ik may be the empty matrix, don't calculate with empty matrices - [nd1111,md1111] = size(D1111); - md1112 = length(D1112); - md1121 = length(D1121); + [nd1111, md1111] = size (D1111); + md1112 = length (D1112); + md1121 = length (D1121); - if ((nd1111 == 0) || (md1112 == 0)) + if (nd1111 == 0) || md1112 == 0) d11hat = -D1122; else - xx = inv(g*g*eye(nz-nu) - D1111*D1111'); + xx = inv (g*g*eye(nz-nu) - D1111*D1111'); d11hat = -D1121*D1111'*xx*D1112 - D1122; endif if (md1112 == 0) - d21hat = eye(ny); + d21hat = eye (ny); elseif (nd1111 == 0) - d21hat = chol(eye(ny) - D1112'*D1112/g/g); + d21hat = chol (eye(ny) - D1112'*D1112/g/g); else - xx = inv(g*g*eye(nz-nu) - D1111*D1111'); - xx = eye(ny) - D1112'*xx*D1112; - d21hat = chol(xx); + xx = inv (g*g*eye(nz-nu) - D1111*D1111'); + xx = eye (ny) - D1112'*xx*D1112; + d21hat = chol (xx); endif if (md1121 == 0) - d12hat = eye(nu); + d12hat = eye (nu); elseif (md1111 == 0) - d12hat = chol(eye(nu) - D1121*D1121'/g/g)'; + d12hat = chol (eye(nu) - D1121*D1121'/g/g)'; else - xx = inv(g*g*eye(nw-ny) - D1111'*D1111); - xx = eye(nu)-D1121*xx*D1121'; - d12hat = chol(xx)'; + xx = inv (g*g*eye(nw-ny) - D1111'*D1111); + xx = eye (nu)-D1121*xx*D1121'; + d12hat = chol (xx)'; endif b2hat = (B2+H12)*d12hat; @@ -138,13 +142,13 @@ ## non-zero D22 is a special case if (d22nz) - if (rank(eye(nu) + d11hat*D22) < nu) + if (rank (eye(nu) + d11hat*D22) < nu) error(" *** cannot compute controller for D22 non-zero."); endif d22new = [D22, zeros(ny,ny); zeros(nu,nu), 0*D22']; - xx = inv(eye(nu+ny) + d22new*dhat); - mhat = inv(eye(nu+ny) + dhat*d22new); + xx = inv (eye(nu+ny) + d22new*dhat); + mhat = inv (eye(nu+ny) + dhat*d22new); ahat = ahat - bhat*((eye(nu+ny)-xx)/dhat)*chat; bhat = bhat*xx; chat = mhat*chat; @@ -152,6 +156,6 @@ endif - K = ss(ahat,bhat(:,1:ny),chat(1:nu,:),dhat(1:nu,1:ny)); + K = ss (ahat, bhat(:,1:ny), chat(1:nu,:), dhat(1:nu,1:ny)); endfunction
--- a/scripts/control/hinf/hinfsyn_chk.m +++ b/scripts/control/hinf/hinfsyn_chk.m @@ -80,6 +80,10 @@ function [retval, Pc, Pf] = hinfsyn_chk (A, B1, B2, C1, C2, D12, D21, g, ptol) + if (nargin != 9) + print_usage (); + endif + Pc = Pf = []; ## Construct the two Hamiltonians @@ -88,39 +92,39 @@ Hf = [ A' , g2*C1'*C1 - C2'*C2; -B1*B1' , -A]; ## check if Hc, Hf are in dom(Ric) - Hcminval = min(abs(real(eig(Hc)))); - Hfminval = min(abs(real(eig(Hf)))); - if(Hcminval < ptol); - disp("hinfsyn_chk: Hc is not in dom(Ric)"); + Hcminval = min (abs (real (eig (Hc)))); + Hfminval = min (abs (real (eig (Hf)))); + if (Hcminval < ptol); + warning ("hinfsyn_chk: Hc is not in dom(Ric)"); retval = 0; return endif if(Hfminval < ptol) - disp("hinfsyn_chk: Hf is not in dom(Ric)"); + warning ("hinfsyn_chk: Hf is not in dom(Ric)"); retval = 0; return endif ## Solve ARE's - Pc = are(A, B2*B2'-g2*B1*B1',C1'*C1); - Pf = are(A',C2'*C2-g2*C1'*C1,B1*B1'); + Pc = are (A, B2*B2'-g2*B1*B1', C1'*C1); + Pf = are (A', C2'*C2-g2*C1'*C1, B1*B1'); - Pceig = eig(Pc); - Pfeig = eig(Pf); - Pcfeig = eig(Pc*Pf); + Pceig = eig (Pc); + Pfeig = eig (Pf); + Pcfeig = eig (Pc*Pf); - if(min(Pceig) < -ptol) - disp("hinfsyn_chk: Pc is not >= 0"); + if (min (Pceig) < -ptol) + warning ("hinfsyn_chk: Pc is not >= 0"); retval = 0; return endif - if(min(Pfeig) < -ptol) - disp("hinfsyn_chk: Pf is not >= 0"); + if (min (Pfeig) < -ptol) + warning ("hinfsyn_chk: Pf is not >= 0"); retval = 0; return endif - if(max(abs(Pcfeig)) >= g*g) - disp("hinfsyn_chk: rho(Pf*Pc) is not < g^2"); + if (max (abs (Pcfeig)) >= g*g) + warning ("hinfsyn_chk: rho(Pf*Pc) is not < g^2"); retval = 0; return endif
--- a/scripts/control/hinf/hinfsyn_ric.m +++ b/scripts/control/hinf/hinfsyn_ric.m @@ -47,39 +47,44 @@ function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol) + if (nargin != 6) + print_usage (); + endif + x_ha_err = 0; # assume success Xinf = []; # default return value - n = issquare(A); - nw = issquare(R); - if(rank(R) != nw) x_ha_err = 6; + n = issquare (A); + nw = issquare (R); + if (rank (R) != nw) + x_ha_err = 6; else # build hamiltonian Ha for X_inf xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB']; Ha = [A, 0*A; -C1'*C1, -A'] - xx; x_ha_err = 0; - [d, Ha] = balance(Ha); - [u, s] = schur(Ha, "A"); - rev = real(eig(s)); + [d, Ha] = balance (Ha); + [u, s] = schur (Ha, "A"); + rev = real (eig (s)); - if (any(abs(rev) <= ptol)) # eigenvalues near the imaginary axis + if (any (abs (rev) <= ptol)) # eigenvalues near the imaginary axis x_ha_err = 1; - elseif (sum(rev > 0) != sum(rev < 0)) + elseif (sum (rev > 0) != sum (rev < 0)) ## unequal number of positive and negative eigenvalues x_ha_err = 2; else ## compute positive Riccati equation solution u = d * u; Xinf = u(n+1:2*n,1:n) / u(1:n,1:n); - if (!all(all(finite(Xinf)))) + if (! all (all (finite (Xinf)))) x_ha_err = 3; - elseif (norm(Xinf-Xinf') >= 10*ptol) + elseif (norm (Xinf-Xinf') >= 10*ptol) ## solution not symmetric x_ha_err = 4; else ## positive semidefinite? ## force symmetry (faster, avoids some convergence problems) Xinf = (Xinf + Xinf')/2; - rev = eig(Xinf); - if (any(rev <= -ptol)) + rev = eig (Xinf); + if (any (rev <= -ptol)) x_ha_err = 5; endif endif
--- a/scripts/control/system/is_sample.m +++ b/scripts/control/system/is_sample.m @@ -28,6 +28,10 @@ function out = is_sample (Ts) + if (nargin != 1) + print_usage (); + endif + out = (isscalar (Ts) && (Ts == abs (Ts)) && (Ts != 0)); endfunction
--- a/scripts/control/system/is_signal_list.m +++ b/scripts/control/system/is_signal_list.m @@ -24,9 +24,13 @@ function flg = is_signal_list (mylist) + if (nargin != 1) + print_usage (); + endif + flg = iscell (mylist); - if(flg) - flg = (rows(mylist) == 1 | columns(mylist) == 1); + if (flg) + flg = (rows (mylist) == 1 || columns (mylist) == 1); end if (flg) for ii = 1:length (mylist)
--- a/scripts/control/system/ss2tf.m +++ b/scripts/control/system/ss2tf.m @@ -59,16 +59,20 @@ function [num, den] = ss2tf (a, b, c, d) + if (nargin != 4) + print_usage (); + endif + ## Check args - [n,m,p] = abcddim(a,b,c,d); + [n, m, p] = abcddim (a, b, c, d); if (n == -1) num = []; den = []; error("ss2tf: Non compatible matrix arguments"); - elseif ( (m != 1) | (p != 1)) + elseif (m != 1 || p != 1) num = []; den = []; - error(["ss2tf: not SISO system: m=",num2str(m)," p=",num2str(p)]); + error ("ss2tf: not SISO system: m=%d, p=%d", m, p); endif if(n == 0) @@ -77,21 +81,21 @@ den = 1; else ## First, get the denominator coefficients - den = poly(a); + den = poly (a); ## Get the zeros of the system - [zz,g] = tzero(a,b,c,d); + [zz, g] = tzero (a, b, c, d); ## Form the Numerator (and include the gain) - if (!isempty(zz)) - num = g * poly(zz); + if (! isempty (zz)) + num = g * poly (zz); else num = g; endif ## the coefficients must be real - den = real(den); - num = real(num); + den = real (den); + num = real (num); endif + endfunction -
--- a/scripts/control/system/sys2fir.m +++ b/scripts/control/system/sys2fir.m @@ -30,6 +30,10 @@ function [c, tsam, inname, outname] = sys2fir (sys) + if (nargin != 1) + print_usage (); + endif + ## let sys2tf do most of the work [num, den, tsam, inname, outname] = sys2tf (sys);
--- a/scripts/control/system/sysgettsam.m +++ b/scripts/control/system/sysgettsam.m @@ -24,8 +24,12 @@ function T = sysgettsam (sys) + if (nargin != 1) + print_usage (); + endif + if (! isstruct (sys)) - print_usage (); + error ("sysgettsam: expecting argument to be system structure"); endif T = sys.tsam;
--- a/scripts/control/system/sysgettype.m +++ b/scripts/control/system/sysgettype.m @@ -39,11 +39,15 @@ function systype = sysgettype (sys) + if (nargin != 1) + print_usage (); + endif + if (! isstruct (sys)) error ("sysgettype: input sys is not a structure"); endif typestr = {"tf", "zp", "ss"}; - systype = typestr{ sys.sys(1) + 1}; + systype = typestr{sys.sys(1) + 1}; endfunction
--- a/scripts/control/system/sysreorder.m +++ b/scripts/control/system/sysreorder.m @@ -44,6 +44,10 @@ function pv = sysreorder (vlen, list) + if (nargin != 2) + print_usage (); + endif + ## disp('sysreorder: entry') pv = 1:vlen;
--- a/scripts/control/system/tf2sys.m +++ b/scripts/control/system/tf2sys.m @@ -63,7 +63,7 @@ ## name changed to tf Feb 2004 function outsys = tf2sys (varargin) - + warning("tf2sys is deprecated. Use tf() instead."); outsys = tf(varargin{:});
--- a/scripts/control/system/zp2tf.m +++ b/scripts/control/system/zp2tf.m @@ -37,18 +37,22 @@ function [num, den] = zp2tf (zer, pol, k) + if (nargin != 3) + print_usage (); + endif + ## Find out whether data was entered as a row or a column vector and ## convert to a column vector if necessary. - [rp,cp] = size(pol); - [rz,cz] = size(zer); + [rp, cp] = size (pol); + [rz, cz] = size (zer); - if(!(isvector(zer) | isempty(zer)) ) - error(sprintf("zer(%dx%d) must be a vector",rz,cz)); - elseif(!(isvector(pol) | isempty(pol)) ) - error(sprintf("pol(%dx%d) must be a vector",rp,cp)); - elseif(length(zer) > length(pol)) - error(sprintf("zer(%dx%d) longer than pol(%dx%d)",rz,cz,rp,cp)); + if (! (isvector (zer) || isempty (zer))) + error ("zer(%dx%d) must be a vector", rz, cz); + elseif(! (isvector (pol) || isempty (pol))) + error ("pol(%dx%d) must be a vector", rp, cp); + elseif (length (zer) > length (pol)) + error ("zer(%dx%d) longer than pol(%dx%d)", rz, cz, rp, cp); endif ## initialize converted polynomials @@ -58,18 +62,24 @@ ## call __zp2ssg2__ if there are complex conjugate pairs left, otherwise ## construct real zeros one by one. Repeat for poles. - while(!isempty(zer)) - if( max(abs(imag(zer))) ) [poly, zer] = __zp2ssg2__ (zer); - else poly = [1, -zer(1)]; - zer = zer(2:length(zer)); endif - num = conv(num,poly); + while(! isempty (zer)) + if (max (abs (imag (zer)))) + [poly, zer] = __zp2ssg2__ (zer); + else + poly = [1, -zer(1)]; + zer = zer(2:length(zer)); + endif + num = conv (num, poly); endwhile - while(!isempty(pol)) - if( max(abs(imag(pol))) ) [poly, pol] = __zp2ssg2__ (pol); - else poly = [1, -pol(1)]; - pol = pol(2:length(pol)); endif - den = conv(den,poly); + while (! isempty (pol)) + if (max (abs (imag (pol)))) + [poly, pol] = __zp2ssg2__ (pol); + else + poly = [1, -pol(1)]; + pol = pol(2:length(pol)); + endif + den = conv (den, poly); endwhile endfunction
--- a/scripts/control/util/axis2dlim.m +++ b/scripts/control/util/axis2dlim.m @@ -39,13 +39,13 @@ function axvec = axis2dlim (axdata) - if(isempty(axdata)) + if (nargin < 1 || isempty (axdata)) axdata = 0; endif ## compute axis limits - minv = min(axdata); - maxv = max(axdata); + minv = min (axdata); + maxv = max (axdata); delv = (maxv-minv)/2; # breadth of the plot midv = (minv + maxv)/2; # midpoint of the plot axmid = [midv(1), midv(1), midv(2), midv(2)];
--- a/scripts/control/util/swap.m +++ b/scripts/control/util/swap.m @@ -31,6 +31,10 @@ function [a1, b1] = swap (a, b) + if (nargin != 2) + print_usage (); + endif + a1 = b; b1 = a;
--- a/scripts/control/util/zgfmul.m +++ b/scripts/control/util/zgfmul.m @@ -33,6 +33,10 @@ function y = zgfmul (a, b, c, d, x) + if (nargin != 5) + print_usage (); + endif + [n,m] = size(b); [p,m1] = size(c); nm = n+m;
--- a/scripts/control/util/zgfslv.m +++ b/scripts/control/util/zgfslv.m @@ -27,6 +27,10 @@ function x = zgfslv (n, m, p, b) + if (nargin != 4) + print_usage (); + endif + nmp = n+m+p; gam1 = (2*n)+m+p; gam2 = n+p; gam3 = n+m;
--- a/scripts/control/util/zginit.m +++ b/scripts/control/util/zginit.m @@ -34,6 +34,10 @@ function zz = zginit (a, b, c, d) + if (nargin != 4) + print_usage (); + endif + [nn,mm] = size(b); [pp,mm] = size(d);
--- a/scripts/control/util/zgreduce.m +++ b/scripts/control/util/zgreduce.m @@ -25,6 +25,10 @@ function retsys = zgreduce (Asys, meps) + if (nargin != 2) + print_usage (); + endif + ## SYS_INTERNAL accesses members of system data structure is_digital(Asys); # make sure it's pure digital/continuous
--- a/scripts/control/util/zgrownorm.m +++ b/scripts/control/util/zgrownorm.m @@ -26,6 +26,10 @@ function [sig, tau] = zgrownorm (mat, meps) + if (nargin != 2) + print_usage (); + endif + rownorm = []; for ii = 1:rows (mat) rownorm(ii) = norm (mat(ii,:));
--- a/scripts/control/util/zgscal.m +++ b/scripts/control/util/zgscal.m @@ -34,6 +34,10 @@ function x = zgscal (a, b, c, d, z, n, m, p) + if (nargin != 8) + print_usage (); + endif + ## initialize parameters: ## Givens rotations, diagonalized 2x2 block of F, gcg vector initialization
--- a/scripts/control/util/zgsgiv.m +++ b/scripts/control/util/zgsgiv.m @@ -30,6 +30,10 @@ function [a, b] = zgsgiv (c, s, a, b) + if (nargin != 4) + print_usage (); + endif + t1 = c*a + s*b; t2 = -s*a + c*b; a = t1;
--- a/scripts/control/util/zgshsr.m +++ b/scripts/control/util/zgshsr.m @@ -38,8 +38,12 @@ function x = zgshsr (y) + if (nargin != 1) + print_usage (); + endif + if (! isvector (y)) - error (sprintf ("y(%dx%d) must be a vector", rows (y), columns (y))); + error ("y(%dx%d) must be a vector", rows (y), columns (y)); endif x = vec (y); m = length (x);
--- a/scripts/general/isa.m +++ b/scripts/general/isa.m @@ -25,5 +25,11 @@ ## Adapted-by: jwe function retval = isa (x, cname) + + if (nargin != 2) + print_usage (); + endif + retval = strcmp (class (x), cname); + endfunction
--- a/scripts/geometry/inpolygon.m +++ b/scripts/geometry/inpolygon.m @@ -38,15 +38,19 @@ function [IN, ON] = inpolygon (X, Y, xv, yv) - if ( !(isreal(X) && isreal(Y) && ismatrix(Y) && ... - ismatrix(Y) && size_equal(X,Y)) ) + if (nargin != 4) + print_usage (); + endif + + if (! (isreal (X) && isreal (Y) && ismatrix (Y) && ismatrix (Y) + && size_equal (X, Y))) error ("inpolygon: first two arguments must be real matrices of same size"); - elseif ( !(isreal(xv) && isreal(yv) && isvector(xv) && ... - isvector(yv) && size_equal(xv,yv)) ) + elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv) + && size_equal (xv, yv))) error ("inpolygon: last two arguments must be real vectors of same size"); endif - npol = length(xv); + npol = length (xv); do_boundary = (nargout >= 2); IN = zeros (size(X), "logical"); @@ -76,6 +80,7 @@ endif j = i; endfor + endfunction %!demo
--- a/scripts/linear-algebra/housh.m +++ b/scripts/linear-algebra/housh.m @@ -35,8 +35,12 @@ function [housv, beta, zer] = housh (x, j, z) - ## check for valid inputs - if (!isvector (x) && ! isscalar (x)) + if (nargin != 3) + print_usage (); + endif + + ## Check for valid inputs. + if (! isvector (x) && ! isscalar (x)) error ("housh: first input must be a vector") elseif (! isscalar(j)) error ("housh: second argment must be an integer scalar")
--- a/scripts/miscellaneous/compare_versions.m +++ b/scripts/miscellaneous/compare_versions.m @@ -70,6 +70,10 @@ function out = compare_versions (v1, v2, operator) + if (nargin != 3) + print_usage (); + endif + ## Make sure that the version numbers are valid. if (! (ischar (v1) && ischar (v2))) error ("compare_versions: both version numbers must be strings"); @@ -184,6 +188,7 @@ ## Reverse the output if not is given. out = xor (not_op, out); + endfunction ## tests
--- a/scripts/miscellaneous/inputname.m +++ b/scripts/miscellaneous/inputname.m @@ -25,7 +25,13 @@ ## @end deftypefn function s = inputname (n) + + if (nargin != 1) + print_usage (); + endif + s = evalin ("caller", sprintf ("deblank (argn(%d,:));", n)); + endfunction ## Warning: heap big magic in the following tests!!!
--- a/scripts/miscellaneous/run.m +++ b/scripts/miscellaneous/run.m @@ -29,6 +29,11 @@ ## PKG_ADD: mark_as_command run function run (s) + + if (nargin != 1) + print_usage (); + endif + [d, f, ext] = fileparts (s); if (! isempty (d)) if (exist (d, "dir"))
--- a/scripts/quaternion/qconj.m +++ b/scripts/quaternion/qconj.m @@ -32,6 +32,10 @@ function retval = qconj (q) + if (nargin != 1 ) + print_usage (); + endif + [a, b, c, d] = quaternion (q); retval = quaternion (-a, -b, -c, d);
--- a/scripts/quaternion/qcoordinate_plot.m +++ b/scripts/quaternion/qcoordinate_plot.m @@ -38,6 +38,10 @@ function qcoordinate_plot (qf, qb, qv) + if (nargin != 3 ) + print_usage (); + endif + degrees = pi / 180; d180 = 180 * degrees;
--- a/scripts/quaternion/qderiv.m +++ b/scripts/quaternion/qderiv.m @@ -43,6 +43,10 @@ function Dmat = qderivmat (Omega) + if (nargin != 1) + print_usage (); + endif + Omega = vec (Omega); if (length (Omega) != 3)
--- a/scripts/quaternion/qderivmat.m +++ b/scripts/quaternion/qderivmat.m @@ -43,6 +43,10 @@ function Dmat = qderivmat (Omega) + if (nargin != 1 ) + print_usage (); + endif + Omega = vec (Omega); if (length (Omega) != 3)
--- a/scripts/quaternion/qinv.m +++ b/scripts/quaternion/qinv.m @@ -32,6 +32,10 @@ function retval = qinv (q) + if (nargin != 1) + print_usage (); + endif + if (norm (q) != 0) retval = qconj (q) / sum (q .* q); else
--- a/scripts/quaternion/qmult.m +++ b/scripts/quaternion/qmult.m @@ -40,7 +40,11 @@ ## Adapted-By: jwe function retval = qmult (a, b) - + + if (nargin != 2 ) + print_usage (); + endif + [a1, b1, c1, d1] = quaternion (a); [a2, b2, c2, d2] = quaternion (b);
--- a/scripts/quaternion/qtrans.m +++ b/scripts/quaternion/qtrans.m @@ -28,6 +28,10 @@ function v = qtrans (v, q) + if (nargin != 2) + print_usage (); + endif + if (! isvector (v) || length (v) != 4) error ("qtrans: v(%d,%d) must be a quaternion", rows (v), columns (v)); elseif (! isvector (q) || length (q) != 4)
--- a/scripts/quaternion/qtransvmat.m +++ b/scripts/quaternion/qtransvmat.m @@ -29,6 +29,10 @@ function Aib = qtransvmat (qib) + if (nargin != 1) + print_usage (); + endif + if (! isvector(qib) || length (qib) != 4) error ("qtransvmat: q(%d,%d) must be a quaternion", rows (qib), \ columns (qib));
--- a/scripts/signal/fractdiff.m +++ b/scripts/signal/fractdiff.m @@ -28,6 +28,10 @@ function retval = fractdiff (x, d) + if (nargin != 2) + print_usage (); + endif + N = 100; if (! isvector (x))
--- a/scripts/signal/freqz_plot.m +++ b/scripts/signal/freqz_plot.m @@ -25,38 +25,42 @@ function freqz_plot (w, h) - n = length (w); + if (nargin != 2) + print_usage (); + endif - ## ## exclude zero-frequency - ## h = h (2 : length (h)); - ## w = w (2 : length (w)); - ## n = n-1; + n = length (w); - mag = 20 * log10 (abs (h)); - phase = unwrap (arg (h)); - maxmag = max (mag); + ## ## exclude zero-frequency + ## h = h (2 : length (h)); + ## w = w (2 : length (w)); + ## n = n-1; - subplot (3, 1, 1); - plot (w, mag); - grid ("on"); - legend("Pass band (dB)"); - axis ([w(1), w(n), maxmag-3, maxmag], "labely"); + mag = 20 * log10 (abs (h)); + phase = unwrap (arg (h)); + maxmag = max (mag); - subplot (3, 1, 2); - plot (w, mag); - grid ("on"); - legend ("Stop band (dB)"); - if (maxmag - min (mag) > 100) - axis ([w(1), w(n), maxmag-100, maxmag], "labely"); - else - axis ("autoy", "labely"); - endif + subplot (3, 1, 1); + plot (w, mag); + grid ("on"); + legend("Pass band (dB)"); + axis ([w(1), w(n), maxmag-3, maxmag], "labely"); - subplot (3, 1, 3); - plot (w, phase*360/(2*pi)); - grid ("on"); - legend ("Phase (degrees)"); - xlabel ("Frequency"); - axis ([w(1), w(n)], "autoy", "label"); + subplot (3, 1, 2); + plot (w, mag); + grid ("on"); + legend ("Stop band (dB)"); + if (maxmag - min (mag) > 100) + axis ([w(1), w(n), maxmag-100, maxmag], "labely"); + else + axis ("autoy", "labely"); + endif + + subplot (3, 1, 3); + plot (w, phase*360/(2*pi)); + grid ("on"); + legend ("Phase (degrees)"); + xlabel ("Frequency"); + axis ([w(1), w(n)], "autoy", "label"); endfunction
--- a/scripts/signal/periodogram.m +++ b/scripts/signal/periodogram.m @@ -28,6 +28,10 @@ function retval = periodogram (x) + if (nargin != 1) + print_usage (); + endif + [r, c] = size(x); if (r == 1)
--- a/scripts/signal/rectangle_lw.m +++ b/scripts/signal/rectangle_lw.m @@ -28,6 +28,10 @@ function retval = rectangle_lw (n, b) + if (nargin != 2) + print_usage (); + endif + retval = zeros (n, 1); t = floor (1 / b);
--- a/scripts/signal/rectangle_sw.m +++ b/scripts/signal/rectangle_sw.m @@ -28,6 +28,10 @@ function retval = rectangle_sw (n, b) + if (nargin != 2) + print_usage (); + endif + retval = zeros (n, 1); retval(1) = 2 / b + 1;
--- a/scripts/signal/sinc.m +++ b/scripts/signal/sinc.m @@ -34,6 +34,10 @@ function result = sinc (x) + if (nargin != 1) + print_usage (); + endif + result = ones (size (x)); i = (x != 0);
--- a/scripts/signal/triangle_lw.m +++ b/scripts/signal/triangle_lw.m @@ -28,6 +28,10 @@ function retval = triangle_lw (n, b) + if (nargin != 2) + print_usage (); + endif + retval = 1 - (0 : n-1)' * b; retval = max ([retval'; (zeros (1, n))])';
--- a/scripts/signal/triangle_sw.m +++ b/scripts/signal/triangle_sw.m @@ -28,6 +28,10 @@ function retval = triangle_sw (n, b) + if (nargin != 2) + print_usage (); + endif + retval = zeros(n,1); retval(1) = 1 / b;
--- a/scripts/signal/yulewalker.m +++ b/scripts/signal/yulewalker.m @@ -31,6 +31,10 @@ function [a, v] = yulewalker (c) + if (nargin != 1) + print_usage (); + endif + p = length (c) - 1; if (columns (c) > 1)
--- a/scripts/sparse/colperm.m +++ b/scripts/sparse/colperm.m @@ -26,6 +26,11 @@ ## @end deftypefn function p = colperm (s) + + if (nargin != 1) + print_usage (); + endif + [i, j] = spfind (s); idx = find (diff ([j; Inf]) != 0); [dummy, p] = sort (idx - [0; idx(1:(end-1))]);
--- a/scripts/sparse/etreeplot.m +++ b/scripts/sparse/etreeplot.m @@ -27,5 +27,10 @@ ## @end deftypefn function etreeplot (s, varargin) + + if (nargin < 1) + print_usage (); + endif + treeplot (etree (s+s'), varargin{:}); endfunction
--- a/scripts/sparse/nonzeros.m +++ b/scripts/sparse/nonzeros.m @@ -22,6 +22,11 @@ ## @end deftypefn function t = nonzeros (s) + + if (nargin != 1) + print_usage (); + endif + if (issparse (s)) [i, j, t] = spfind (s); else
--- a/scripts/sparse/spalloc.m +++ b/scripts/sparse/spalloc.m @@ -41,5 +41,10 @@ ## @end deftypefn function s = spalloc (r, c, nz) + + if (nargin < 2) + print_usage (); + endif + s = sparse (r, c); endfunction
--- a/scripts/sparse/spones.m +++ b/scripts/sparse/spones.m @@ -23,6 +23,11 @@ ## @end deftypefn function s = spones (s) + + if (nargin != 1) + print_usage (); + endif + if (issparse (s)) [i, j, v, m, n] = spfind (s); else
--- a/scripts/sparse/spy.m +++ b/scripts/sparse/spy.m @@ -30,6 +30,10 @@ function spy (S, varargin) + if (nargin < 1) + print_usage (); + endif + markersize = NaN; if (numel (i) < 1000) LineSpec = "*";
--- a/scripts/specfun/isprime.m +++ b/scripts/specfun/isprime.m @@ -35,6 +35,11 @@ ## @end deftypefn function t = isprime (n) + + if (nargin < 1) + print_usage (); + endif + if (! isscalar (n)) nel = numel (n); t = n;
--- a/scripts/statistics/distributions/empirical_cdf.m +++ b/scripts/statistics/distributions/empirical_cdf.m @@ -28,6 +28,10 @@ function cdf = empirical_cdf (x, data) + if (nargin != 2) + print_usage (); + endif + if (! isvector (data)) error ("empirical_cdf: data must be a vector"); endif
--- a/scripts/statistics/distributions/empirical_inv.m +++ b/scripts/statistics/distributions/empirical_inv.m @@ -28,6 +28,10 @@ function inv = empirical_inv (x, data) + if (nargin != 2) + print_usage (); + endif + if (! isvector (data)) error ("empirical_inv: data must be a vector"); endif
--- a/scripts/statistics/distributions/empirical_pdf.m +++ b/scripts/statistics/distributions/empirical_pdf.m @@ -28,6 +28,10 @@ function pdf = empirical_pdf (x, data) + if (nargin != 2) + print_usage (); + endif + if (! isvector (data)) error ("empirical_pdf: data must be a vector"); endif
--- a/scripts/statistics/models/logistic_regression_derivatives.m +++ b/scripts/statistics/models/logistic_regression_derivatives.m @@ -29,6 +29,10 @@ function [dl, d2l] = logistic_regression_derivatives (x, z, z1, g, g1, p) + if (nargin != 6) + print_usage (); + endif + ## first derivative v = g .* (1 - g) ./ p; v1 = g1 .* (1 - g1) ./ p; dlogp = [(dmult (v, z) - dmult (v1, z1)), (dmult (v - v1, x))]; @@ -39,4 +43,4 @@ d2l = [z, x]' * dmult (w, [z, x]) - [z1, x]' * dmult (w1, [z1, x]) ... - dlogp' * dlogp; -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/models/logistic_regression_likelihood.m +++ b/scripts/statistics/models/logistic_regression_likelihood.m @@ -29,6 +29,10 @@ function [g, g1, p, dev] = logistic_regression_likelihood (y, x, beta, z, z1) + if (nargin != 5) + print_usage (); + endif + e = exp ([z, x] * beta); e1 = exp ([z1, x] * beta); g = e ./ (1 + e); g1 = e1 ./ (1 + e1); g = max (y == max (y), g); g1 = min (y > min(y), g1);