Mercurial > forge
changeset 8454:bc4257a29d2a octave-forge
impinvar/invimpinvar: new functions for signal package submitted by Rudy Eschauzier <reschauzier@yahoo.com>
| author | carandraug |
|---|---|
| date | Fri, 30 Sep 2011 13:42:44 +0000 |
| parents | c8750827ad19 |
| children | 062e063e7161 |
| files | main/signal/inst/impinvar.m main/signal/inst/invimpinvar.m main/signal/inst/private/h1_z_deriv.m main/signal/inst/private/inv_residue.m main/signal/inst/private/pad_poly.m main/signal/inst/private/polyrev.m main/signal/inst/private/to_real.m |
| diffstat | 7 files changed, 390 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/impinvar.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,103 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts}, @var{tol}) +## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts}) +## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}) +## Converts analog filter with coefficients @var{b} and @var{a} to digital, +## conserving impulse response. +## +## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz. +## +## If @var{tol} is not specified, it defaults to 0.0001 (0.1%) +## This function does the inverse of impinvar so that the following example should +## restore the original values of @var{a} and @var{b}. +## +## @command{invimpinvar} implements the reverse of this function. +## @example +## [b, a] = impinvar (b, a); +## [b, a] = invimpinvar (b, a); +## @end example +## +## @seealso{bilinear, invimpinvar} +## @end deftypefn + +function [b_out, a_out] = impinvar (b_in, a_in, ts = 1, tol = 0.0001) + + if (nargin <2) + print_usage; + endif + + ## to be compatible with the matlab implementation where an empty vector can + ## be used to get the default + if (isempty(ts)) + ts = 1; + endif + + [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 + +## 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); + + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/invimpinvar.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,98 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts}, @var{tol}) +## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts}) +## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}) +## Converts digital filter with coefficients @var{b} and @var{a} to analog, +## conserving impulse response. +## +## This function does the inverse of impinvar so that the following example should +## restore the original values of @var{a} and @var{b}. +## @example +## [b, a] = impinvar (b, a); +## [b, a] = invimpinvar (b, a); +## @end example +## +## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz. +## +## If @var{tol} is not specified, it defaults to 0.0001 (0.1%) +## +## @seealso{bilinear, impinvar} +## @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) + + if (nargin <2) + print_usage; + endif + + [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 + +## Inverse function of z_res (see impinvar source) +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 + + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/private/h1_z_deriv.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,53 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## 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. +## 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) + + %% 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 + +## 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. +function b = h1_deriv(n,ts) + b = factorial(n)*bincoeff(n,0:n); % Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc. + b *= (-1)^n; +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/private/inv_residue.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,62 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## This function is necessary for impinvar and invimpinvar of the signal package + +## Inverse of Octave residue function +function [b_out, a_out] = inv_residue(r_in, p_in, k_in) + + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/private/pad_poly.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,25 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## This function is necessary for impinvar and invimpinvar of the signal package + +## Pad polynomial to length n +function p_out = pad_poly(p_in,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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/private/polyrev.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,26 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## This function is necessary for impinvar and invimpinvar of the signal package + +## Reverse the coefficients of a polynomial +function p_out = polyrev (p_in) + n = length(p_in); + for i=(1:n) + p_out(i) = p_in(n-i+1); + endfor +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/private/to_real.m Fri Sep 30 13:42:44 2011 +0000 @@ -0,0 +1,23 @@ +## Copyright 2007 R.G.H. Eschauzier <reschauzier@yahoo.com> +## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com> +## +## 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 + +## This function is necessary for impinvar and invimpinvar of the signal package + +## Round complex number to nearest real number +function p_out = to_real(p_in) + p_out = abs(p_in) .* sign(real(p_in)); +endfunction