[Octave-cvsupdate] SF.net SVN: octave:[8728] trunk/octave-forge/main/signal

[<car...@us...>]

Revision: 8728
          http://octave.svn.sourceforge.net/octave/?rev=8728&view=rev
Author:   carandraug
Date:     2011-10-11 14:10:16 +0000 (Tue, 11 Oct 2011)
Log Message:
-----------
impinvar: improvements contributed by Rudy Eschauzier <res...@ya...>
  * test cases for both impinvar and invimpvar
  * bug fixes
  * vectorization for faster code
  * tests make package dependent on control

NOTE: values below tolerance are no longer truncated but converted to zero

Modified Paths:
--------------
    trunk/octave-forge/main/signal/DESCRIPTION
    trunk/octave-forge/main/signal/inst/impinvar.m
    trunk/octave-forge/main/signal/inst/invimpinvar.m
    trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
    trunk/octave-forge/main/signal/inst/private/inv_residue.m

Modified: trunk/octave-forge/main/signal/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/signal/DESCRIPTION	2011-10-10 16:19:34 UTC (rev 8727)
+++ trunk/octave-forge/main/signal/DESCRIPTION	2011-10-11 14:10:16 UTC (rev 8728)
@@ -5,7 +5,7 @@
 Maintainer: The Octave Community
 Title: Signal Processing.
 Description: Signal processing tools, including filtering, windowing and display functions.
-Depends: octave (> 2.9.9), optim (>= 1.0.0), specfun
+Depends: octave (> 2.9.9), optim (>= 1.0.0), specfun, control
 Autoload: yes
 License: GPL version 2 or later
 Url: http://octave.sf.net

Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m	2011-10-10 16:19:34 UTC (rev 8727)
+++ trunk/octave-forge/main/signal/inst/impinvar.m	2011-10-11 14:10:16 UTC (rev 8728)
@@ -41,7 +41,7 @@
 ## @seealso{bilinear, invimpinvar}
 ## @end deftypefn
 
-function [b_out, a_out] = impinvar (b_in, a_in, ts = 1, tol = 0.0001)
+function [b_out, a_out] = impinvar (b_in, a_in, fs = 1, tol = 0.0001)
 
   if (nargin <2)
     print_usage;
@@ -49,10 +49,10 @@
 
   ## to be compatible with the matlab implementation where an empty vector can
   ## be used to get the default
-  if (isempty(ts))
+  if (isempty(fs))
     ts = 1;
   else
-    ts = 1/ts; # we should be using sampling frequencies to be compatible with Matlab
+    ts = 1/fs; # we should be using sampling frequencies to be compatible with Matlab
   endif
 
   [r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
@@ -86,17 +86,19 @@
   [b_out, a_out] = inv_residue(r_out, p_out, k_out, tol);
   a_out          = to_real(a_out); % Get rid of spurious imaginary part
   b_out          = to_real(b_out);
-  b_out          = polyreduce(b_out);
 
+  % Shift results right to account for calculating in z instead of z^-1
+  b_out(end)=[];
+
   ## respect the required tolerance values
-  b_out(abs(b_out)<tol) = [];
-  a_out(abs(a_out)<tol) = [];
+  b_out(abs(b_out)<tol) = 0;
+  a_out(abs(a_out)<tol) = 0;
 
 endfunction
 
 ## Convert residue vector for single and multiple poles in s-domain (located at sm) to
 ## residue vector in z-domain. The variable k is the direct term of the result.
-function [r_out, p_out, k_out] = z_res (r_in, sm, ts);
+function [r_out, p_out, k_out] = z_res (r_in, sm, ts)
 
   p_out = exp(ts * sm); % z-domain pole
   n     = length(r_in); % Multiplicity of the pole
@@ -111,3 +113,39 @@
   endfor
 
 endfunction
+
+
+%!function err = stozerr(bs,as,fs)
+%!
+%!  % number of time steps
+%!  n=10;
+%!
+%!  % impulse invariant transform to z-domain
+%!  [bz az]=impinvar(bs,as,fs);
+%!
+%!  % create sys object of transfer function
+%!  s=tf(bs,as);
+%!
+%!  % calculate impulse response of continuous time system
+%!  % at discrete time intervals 1/fs
+%!  ys=impulse(s,1,(n-1)/fs,n);
+%!
+%!  % impulse response of discrete time system
+%!  yz=filter(bz,az,[1 zeros(1,n-1)]);
+%!
+%!  % find rms error
+%!  err=sqrt(sum((yz*fs.-ys).^2)/length(ys));
+%!  endfunction
+%!
+%!assert(stozerr([1],[1 1],0.1),0,0.0001);
+%!assert(stozerr([1],[1 2 1],0.1),0,0.0001);
+%!assert(stozerr([1 1],[1 2 1],0.1),0,0.0001);
+%!assert(stozerr([1],[1 3 3 1],0.1),0,0.0001);
+%!assert(stozerr([1 1],[1 3 3 1],0.1),0,0.0001);
+%!assert(stozerr([1 1 1],[1 3 3 1],0.1),0,0.0001);
+%!assert(stozerr([1],[1 0 1],0.1),0,0.0001);
+%!assert(stozerr([1 1],[1 0 1],0.1),0,0.0001);
+%!assert(stozerr([1],[1 0 2 0 1],0.1),0,0.0001);
+%!assert(stozerr([1 1],[1 0 2 0 1],0.1),0,0.0001);
+%!assert(stozerr([1 1 1],[1 0 2 0 1],0.1),0,0.0001);
+%!assert(stozerr([1 1 1 1],[1 0 2 0 1],0.1),0,0.0001);

Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m	2011-10-10 16:19:34 UTC (rev 8727)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m	2011-10-11 14:10:16 UTC (rev 8728)
@@ -40,7 +40,7 @@
 ## @end deftypefn
 
 ## Impulse invariant conversion from s to z domain
-function [b_out, a_out] = invimpinvar (b_in, a_in, ts = 1, tol = 0.0001)
+function [b_out, a_out] = invimpinvar (b_in, a_in, fs = 1, tol = 0.0001)
 
   if (nargin <2)
     print_usage;
@@ -48,12 +48,14 @@
 
   ## to be compatible with the matlab implementation where an empty vector can
   ## be used to get the default
-  if (isempty(ts))
+  if (isempty(fs))
     ts = 1;
   else
-    ts = 1/ts; # we should be using sampling frequencies to be compatible with Matlab
+    ts = 1/fs; # we should be using sampling frequencies to be compatible with Matlab
   endif
 
+  b_in = [b_in 0]; %so we can calculate in z instead of z^-1
+
   [r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
 
   n = length(r_in); % Number of poles/residues
@@ -83,6 +85,11 @@
   [b_out, a_out] = inv_residue(r_out, sm_out , 0, tol);
   a_out          = to_real(a_out);      % Get rid of spurious imaginary part
   b_out          = to_real(b_out);
+
+  ## respect the required tolerance values
+  b_out(abs(b_out)<tol) = 0;
+  a_out(abs(a_out)<tol) = 0;
+
   b_out          = polyreduce(b_out);
 
 endfunction
@@ -91,19 +98,57 @@
 function [r_out sm_out k_out] = inv_z_res (r_in,p_in,ts)
 
   n    = length(r_in); % multiplicity of the pole
-  r_in = r_in';        % From column vector to row vector
+  r_in = r_in.';       % From column vector to row vector
 
-  i=n;
-  while (i>1) % Go through residues starting from highest order down
-    r_out(i)   = r_in(i) / ((ts * p_in)^i);                     % Back to binomial coefficient for highest order (always 1)
-    r_in(1:i) -= r_out(i) * polyrev(h1_z_deriv(i-1,p_in,ts)); % Subtract highest order result, leaving r_in(i) zero
-    i--;
+  j=n;
+  while (j>1) % Go through residues starting from highest order down
+    r_out(j)   = r_in(j) / ((ts * p_in)^j);                   % Back to binomial coefficient for highest order (always 1)
+    r_in(1:j) -= r_out(j) * polyrev(h1_z_deriv(j-1,p_in,ts)); % Subtract highest order result, leaving r_in(j) zero
+    j--;
   endwhile
 
   %% Single pole (no multiplicity)
   r_out(1) = r_in(1) / ((ts * p_in));
   k_out    = r_in(1) / p_in;
-
   sm_out   = log(p_in) / ts;
 
 endfunction
+
+
+%!function err = ztoserr(bz,az,fs)
+%!
+%!  % number of time steps
+%!  n=10;
+%!
+%!  % make sure system is realizable (no delays)
+%!  bz=prepad(bz,length(az)-1,0,2);
+%!
+%!  % inverse impulse invariant transform to s-domain
+%!  [bs as]=invimpinvar(bz,az,fs);
+%!
+%!  % create sys object of transfer function
+%!  s=tf(bs,as);
+%!
+%!  % calculate impulse response of continuous time system
+%!  % at discrete time intervals 1/fs
+%!  ys=impulse(s,1,(n-1)/fs,n);
+%!
+%!  % impulse response of discrete time system
+%!  yz=filter(bz,az,[1 zeros(1,n-1)]);
+%!
+%!  % find rms error
+%!  err=sqrt(sum((yz*fs.-ys).^2)/length(ys));
+%!  endfunction
+%!
+%!assert(ztoserr([1],[1 -0.5],10),0,0.0001);
+%!assert(ztoserr([1],[1 -1 0.25],10),0,0.0001);
+%!assert(ztoserr([1 1],[1 -1 0.25],10),0,0.0001);
+%!assert(ztoserr([1],[1 -1.5 0.75 -0.125],10),0,0.0001);
+%!assert(ztoserr([1 1],[1 -1.5 0.75 -0.125],10),0,0.0001);
+%!assert(ztoserr([1 1 1],[1 -1.5 0.75 -0.125],10),0,0.0001);
+%!assert(ztoserr([1],[1 0 0.25],10),0,0.0001);
+%!assert(ztoserr([1 1],[1 0 0.25],10),0,0.0001);
+%!assert(ztoserr([1],[1 0 0.5 0 0.0625],10),0,0.0001);
+%!assert(ztoserr([1 1],[1 0 0.5 0 0.0625],10),0,0.0001);
+%!assert(ztoserr([1 1 1],[1 0 0.5 0 0.0625],10),0,0.0001);
+%!assert(ztoserr([1 1 1 1],[1 0 0.5 0 0.0625],10),0,0.0001);

Modified: trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m	2011-10-10 16:19:34 UTC (rev 8727)
+++ trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m	2011-10-11 14:10:16 UTC (rev 8728)
@@ -17,7 +17,7 @@
 
 ## This function is necessary for impinvar and invimpinvar of the signal package
 
-## Find {-zd/dz}^n*H1(z). I.e., first multiply by -z, then diffentiate, then multiply by -z, etc.
+## Find {-zd/dz}^n*H1(z). I.e., first differentiate, then multiply by -z, then differentiate, etc.
 ## The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)).
 ## Works for n>0.
 function b = h1_z_deriv(n, p, ts)
@@ -25,22 +25,22 @@
   %% Build the vector d that holds coefficients for all the derivatives of H1(z)
   %% The results reads d(n)*z^(1)*(d/dz)^(1)*H1(z) + d(n-1)*z^(2)*(d/dz)^(2)*H1(z) +...+ d(1)*z^(n)*(d/dz)^(n)*H1(z)
   d = (-1)^n; % Vector with the derivatives of H1(z)
-  for i=(1:n-1)
-    d  = conv(d,[1 0]);               % Shift result right (multiply by -z)
-    d += prepad(polyderiv(d), i+1, 0) % Add the derivative
+  for i= (1:n-1)
+    d  = [d 0];                           % Shift result right (multiply by -z)
+    d += prepad(polyderiv(d), i+1, 0, 2); % Add the derivative
   endfor
 
   %% Build output vector
-  b = zeros(1,n+1);
-  for i=(1:n)
-    b += d(i) * prepad(h1_deriv(n-i+1), n+1, 0);
+  b = zeros (1, n + 1);
+  for i = (1:n)
+    b += d(i) * prepad(h1_deriv(n-i+1), n+1, 0, 2);
   endfor
 
   b *= ts^(n+1)/factorial(n);
-  for i=(1:n+1)
-    b(n-i+2) *= p^i; % Multiply coefficients with p^i, where i is the index of the coeff.
-  endfor
 
+  %% Multiply coefficients with p^i, where i is the index of the coeff.
+  b.*=p.^(n+1:-1:1);
+
 endfunction
 
 ## Find (z^n)*(d/dz)^n*H1(z), where H1(z)=ts*z/(z-p), ts=sampling period,

Modified: trunk/octave-forge/main/signal/inst/private/inv_residue.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/inv_residue.m	2011-10-10 16:19:34 UTC (rev 8727)
+++ trunk/octave-forge/main/signal/inst/private/inv_residue.m	2011-10-11 14:10:16 UTC (rev 8728)
@@ -37,7 +37,7 @@
   while (i<=n)
     term   = [1 -p_in(i)];               % Term to be factored out
     p      = r_in(i)*deconv(a_out,term); % Residue times resulting polynomial
-    p      = prepad(p, n+1, 0);          % Pad for proper length
+    p      = prepad(p, n+1, 0, 2);       % Pad for proper length
     b_out += p;
 
     m          = 1;
@@ -48,10 +48,9 @@
        m++;
        mterm  = conv(mterm, term);              % Next multiplicity to be factored out
        p      = r_in(i) * deconv(a_out, mterm); % Resulting polynomial
-       p      = prepad(p, n+1, 0);              % Pad for proper length
+       p      = prepad(p, n+1, 0, 2);           % Pad for proper length
        b_out += p;
     endwhile
   i++;
   endwhile
-  b_out = polyreduce(b_out);
 endfunction

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.