diff options
Diffstat (limited to 'macros/impinvar.sci')
| -rw-r--r-- | macros/impinvar.sci | 187 |
1 files changed, 145 insertions, 42 deletions
diff --git a/macros/impinvar.sci b/macros/impinvar.sci index cc9aa6f..57fdb88 100644 --- a/macros/impinvar.sci +++ b/macros/impinvar.sci @@ -1,43 +1,146 @@ -function [b_out, a_out] = impinvar (b, a, fs, tol) -//This function converts analog filter with coefficients b and a to digital, conserving impulse response. -//Calling Sequence -//[b, a] = impinvar (b, a) -//[b, a] = impinvar (b, a, fs) -//[b, a] = impinvar (b, a, fs, tol) -//Parameters -//b: real or complex valued scalar or vector -//a: real or complex valued scalar or vector, order should be greater than b -//fs: real or complex value, default value 1Hz -//tol: real or complex value, default value 0.0001 -//Description -//This is an Octave function. -//This function converts analog filter with coefficients b and a to digital, conserving impulse response. -//This function does the inverse of impinvar. -//Examples -//b = 0.0081000 -//a = [2.0000000, 0.56435378, 0.4572792, 0.00705544, 0.091000] -//[ay,by] = impinvar(b,a,10) -//ay = -// 0.0000e+00 7.5293e-08 2.9902e-07 7.4238e-08 -//by = -// 1.00000 -3.96992 5.91203 -3.91428 0.97218 - -funcprot(0); -rhs = argn(2) -if(rhs<2) -error("Wrong number of input arguments.") -end - - - select(rhs) - case 2 then -// [b, a] = callOctave("impinvar",b,a) - [b_out, a_out] = callOctave("impinvar",b,a) - case 3 then -// [b, a] = callOctave("impinvar",b,a,fs) - [b_out, a_out] = callOctave("impinvar",b,a,fs) - case 4 then -// [b, a] = callOctave("impinvar",b,a,fs,tol) - [b_out, a_out] = callOctave("impinvar",b,a,fs,tol) - end +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh , SOE CUSAT +// Last Modified on : Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +/* +Calling Sequence +[b_out, a_out] = impinvar (b, a, fs, tol) +[b_out, a_out] = impinvar (b, a, fs) +[b_out, a_out] = impinvar (b, a) +Converts analog filter with coefficients b and a to digital, conserving impulse response. + +If fs is not specified, or is an empty vector, it defaults to 1Hz. + +If 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 a and b. + +invimpinvar implements the reverse of this function. + +[b, a] = impinvar (b, a); +[b, a] = invimpinvar (b, a); +Reference: Thomas J. Cavicchi (1996) “Impulse invariance and multiple-order poles”. IEEE transactions on signal processing, Vol 44 (9): 2344–2347 +Dependencies + residue + + +*/ +function [b_out, a_out] = impinvar (b_in, a_in, fs , tol) + error("impinvar: Missing functionality .This function is not implemented yet. Will available in next release"); +endfunction +/* + if (nargin <2) + error ("impinvar: Insufficient input arguments"); + end + if nargin < 3 + fs = 1; + end + if nargin < 4 + tol = 0.0001; + end + // to be compatible with the matlab implementation where an empty vector can + // be used to get the default + if (isempty(fs)) + ts = 1; + else + ts = 1/fs; // we should be using sampling frequencies to be compatible with Matlab + end + + [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"); + end + + r_out = zeros(1,n); // Residues of H(z) + p_out = zeros(1,n); // Poles of H(z) + k_out = 0; // Constant 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=i+1; // Next residue + m=m+1; // Next multiplicity + end + [r, p, k] = z_res(r_in(i-m+1:i), first_pole, ts); // Find z-domain residues + k_out = 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=i+1; // Next s-domain residue/pole + end + + [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); + + // Shift results right to account for calculating in z instead of z^-1 + 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_out(1:i) + r_in(i) * polyrev(h1_z_deriv(i-1, p_out, ts)); // Add z-domain residues + end + +endfunction +function p_out = polyrev (p_in) + + p_out = p_in($:-1:1); + + +endfunction +function p_out = to_real(p_in) + + p_out = abs(p_in) .* sign(real(p_in)); + + endfunction +/* +//test passed +[b_out,a_out]=impinvar([1],[1 1],100); +assert_checkalmostequal(b_out,0.01,%eps,1e-4) +assert_checkalmostequal(a_out,[1 -0.99],%eps,1e-4) + +//test passed +[b_out,a_out]=impinvar([1],[1 2 1],100) +assert_checkalmostequal(b_out,[0 9.9005e-5],%eps,1e-4) +assert_checkalmostequal(a_out,[1 -1.9801 0.9802],%eps,1e-4) + +[b_out, a_out] = impinvar([1], [1 3 3 1], 100) // test passed + +[b_out, a_out] = impinvar([1 1], [1 3 3 1], 100) // test passed + +[b_out, a_out] = impinvar([1 1 1], [1 3 3 1], 100) // test passed + +// FIXME : builtin filter doesn't accepts complex parameters +// These test cases will through errors +// [b_out, a_out] = impinvar([1], [1 0 1], 100) +// [b_out, a_out] = impinvar([1 1], [1 0 1], 100) +// [b_out, a_out] = impinvar([1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1 1], [1 0 2 0 1], 100) +// [b_out, a_out] = impinvar([1 1 1 1], [1 0 2 0 1], 100) + +*/ |