diff -urN impinvar-ORG/impinvar.m impinvar/impinvar.m --- impinvar-ORG/impinvar.m 2011-10-10 19:41:21.000000000 +0200 +++ impinvar/impinvar.m 2011-10-10 19:34:06.000000000 +0200 @@ -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)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); diff -urN impinvar-ORG/private/h1_z_deriv.m impinvar/private/h1_z_deriv.m --- impinvar-ORG/private/h1_z_deriv.m 2011-10-10 19:41:39.000000000 +0200 +++ impinvar/private/h1_z_deriv.m 2011-10-10 19:34:25.000000000 +0200 @@ -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) @@ -26,20 +26,20 @@ %% 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 + 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 += 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 diff -urN impinvar-ORG/private/inv_residue.m impinvar/private/inv_residue.m --- impinvar-ORG/private/inv_residue.m 2011-10-10 19:41:39.000000000 +0200 +++ impinvar/private/inv_residue.m 2011-10-10 19:34:25.000000000 +0200 @@ -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