Hi Juan,
See attached file. Please note that although I feel the algorithms used for the
transform are sound, the implementation in Octave isn't the most efficient. The
intention was to write a reference implementation that proves the concept and
can serve as a basis for improvements.
Regards,
Rudy.
________________________________
From: Juan Manuel Parrilla Gutiérrez <joanma...@gmail.com>
To: reschauz...@yahoo.com
Cc: octave-dev@lists.sourceforge.net
Sent: Wednesday, September 28, 2011 3:19 PM
Subject: Impulse invariance transform implementation
Hi Rudy,
I see here: http://comments.gmane.org/gmane.comp.gnu.octave.sources/116 that
you implemented "impulse invariance transform" for octave.
I'd like to use your implementation with octave. Also the octave community
would be interested in adding your code.
Please, reply as "reply all" in this e-mail so that your answer also arrives to
the list of octave.
Thank you very much.
--
Juan Manuel Parrilla Gutiérrez
1;
% Copyright 2007 R.G.H. Eschauzier
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
function [b_out a_out]=c2d_imp(b_in,a_in,ts)
% Impulse invariant conversion from s to z domain
tol=0.0001;
[r_in p_in k_in]=residue(b_in,a_in); % partial fraction expansion
n=length(r_in); % Number of poles/residues
if (length(k_in)>0) % Greater than zero means we cannot do impulse invariance
error("Order numerator >= order denominator");
endif
r_out=zeros(1,n); % Residues of H(z)
p_out=zeros(1,n); % Poles of H(z)
k_out=0; % Contstant term of H(z)
i=1;
while (i<=n)
m=1;
first_pole=p_in(i); % Pole in the s-domain
while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
i++; %Next residue
m++; % Next multiplicity
endwhile
[r p k]=z_res(r_in(i-m+1:i),first_pole,ts); % Find z-domain residues
k_out+=k; % Add direct term to output
p_out(i-m+1:i)=p; % Copy z-domain pole(s) to output
r_out(i-m+1:i)=r; % Copy z-domain residue(s) to output
i++; % Next s-domain residue/pole
endwhile
[b_out a_out]=inv_residue(r_out,p_out,k_out);
a_out=to_real(a_out); % Get rid of spurious imaginary part
b_out=to_real(b_out);
b_out=polyreduce(b_out);
endfunction
function [b_out a_out]=d2c_imp(b_in,a_in,ts)
% Impulse invariant conversion from s to z domain
tol=0.0001;
[r_in p_in k_in]=residue(b_in,a_in); % partial fraction expansion
n=length(r_in); % Number of poles/residues
if (length(k_in)>1) % Greater than one means we cannot do impulse invariance
error("Order numerator > order denominator");
endif
r_out=zeros(1,n); % Residues of H(s)
sm_out=zeros(1,n); % Poles of H(s)
i=1;
while (i<=n)
m=1;
first_pole=p_in(i); % Pole in the z-domain
while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
i++; %Next residue
m++; % Next multiplicity
endwhile
[r sm k]=inv_z_res(r_in(i-m+1:i),first_pole,ts); % Find s-domain residues
k_in-=k; % Just to check, should end up zero for physical system
sm_out(i-m+1:i)=sm; % Copy s-domain pole(s) to output
r_out(i-m+1:i)=r; % Copy s-domain residue(s) to output
i++; % Next z-domain residue/pole
endwhile
[b_out a_out]=inv_residue(r_out,sm_out,0);
a_out=to_real(a_out); % Get rid of spurious imaginary part
b_out=to_real(b_out);
b_out=polyreduce(b_out);
endfunction
function b=h1_deriv(n,ts)
% Find (z^n)*(d/dz)^n*H1(z), where H1(z)=ts*z/(z-p), ts=sampling period,
% p=exp(sm*ts) and sm is the s-domain pole with multiplicity n+1.
% 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)),
% where b(i) is the binomial coefficient bincoeff(n,i) times n!. Works for n>0.
b=factorial(n)*bincoeff(n,0:n); % Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc.
b*=(-1)^n;
endfunction
function b=h1_z_deriv(n,p,ts)
% Find {-zd/dz}^n*H1(z). I.e., first multiply by -z, then diffentiate, then multiply by -z, 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.
% 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+=pad_poly(polyderiv(d),i+1); % Add the derivative
endfor
% Build output vector
b=zeros(1,n+1);
for i=(1:n)
b+=d(i)*pad_poly(h1_deriv(n-i+1,ts),n+1);
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
endfunction
function p_out=polyrev(p_in)
% Reverse the coefficients of a polynomial
n=length(p_in);
for i=(1:n)
p_out(i)=p_in(n-i+1);
endfor
endfunction
function [r_out p_out k_out]=z_res(r_in,sm,ts);
% 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.
p_out=exp(ts*sm); % z-domain pole
n=length(r_in); % Multiplicity of the pole
r_out=zeros(1,n); % Residue vector
% First pole (no multiplicity)
k_out=r_in(1)*ts; % PFE of z/(z-p) = p/(z-p)+1; direct part
r_out(1)=r_in(1)*ts*p_out; % pole part of PFE
for i=(2:n) % Go through s-domain residues for multiple pole
r_out(1:i)+=r_in(i)*polyrev(h1_z_deriv(i-1,p_out,ts)); % Add z-domain residues
endfor
endfunction
function [r_out sm_out k_out]=inv_z_res(r_in,p_in,ts)
% Inverse function of z_res
n=length(r_in); % multiplicity of the pole
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--;
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 p_out=pad_poly(p_in,n)
% Pad polynomial to length n
l=length(p_in); % Length of input polynomial
p_out(n-l+1:n)=p_in(1:l); % Right shift for proper length
p_out(1:n-l)=0; % Set first elements to zero
endfunction
function [b_out a_out]=inv_residue(r_in,p_in,k_in)
% Inverse of Octave residue function
tol=0.0001;
n=length(r_in); % Number of poles/residues
k=0; % Capture contstant term
if (length(k_in)==1) % A single direct term (order N = order D)
k=k_in(1); % Capture constant term
elseif (length(k_in)>1) % Greater than one means non-physical system
error("Order numerator > order denominator");
endif
a_out=1; % Calculate denominator
for i=(1:n)
a_out=conv(a_out,[1 -p_in(i)]);
endfor
b_out=zeros(1,n+1);
b_out+=k*a_out; % Constant term: add k times denominator to numerator
i=1;
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=pad_poly(p,n+1); % Pad for proper length
b_out+=p;
m=1;
mterm=term;
first_pole=p_in(i);
while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
i++; %Next residue
m++;
mterm=conv(mterm,term); % Next multiplicity to be factored out
p=r_in(i)*deconv(a_out,mterm); % Resulting polynomial
p=pad_poly(p,n+1); % Pad for proper length
b_out+=p;
endwhile
i++;
endwhile
b_out=polyreduce(b_out);
endfunction
function p_out=to_real(p_in)
% Round complex number to nearest real number
p_out=abs(p_in).*sign(real(p_in));
endfunction------------------------------------------------------------------------------
All the data continuously generated in your IT infrastructure contains a
definitive record of customers, application performance, security
threats, fraudulent activity and more. Splunk takes this data and makes
sense of it. Business sense. IT sense. Common sense.
http://p.sf.net/sfu/splunk-d2dcopy1
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
