FOSSEE Winter Internship 2024 work done by Abinash Singh

13 files changed, 1936 insertions, 886 deletions

diff --git a/macros/deconv.sci b/macros/deconv.sci
index 1c4d8fa..08ea92b 100644
--- a/macros/deconv.sci
+++ b/macros/deconv.sci
@@ -1,104 +1,106 @@
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
 function [b, r] = deconv (y, a)
-// calling sequence:
-// [b,r]= deconv (y, a)
-// Deconvolve two vectors.
-//
-// [b, r] = deconv (y, a) solves for b and r such that
-// y = conv (a, b) + r.
-//
-// If y and a are polynomial coefficient vectors, b will
-// contain the coefficients of the polynomial quotient and r will be
-// a remainder polynomial of lowest order.
 
-//Test cases:
-//1.
-//[b, r] = deconv ([3, 6, 9, 9], [1, 2, 3])
-//Output:
-//b=[3, 0]
-//r=[0, 0, 0, 9]
+  if (nargin ~= 2)
+    error("deconv : Two arguments are required ");
+  end
 
-//2.
-//[b, r] = deconv ([3, 6], [1; 2; 3])
-//Output:
-//b  =  0.
-//r= [- 2. ; 8]
+    if (~isvector(y) && ~isscalar(y))
+        error ("deconv: Y must be vector");
+    end
+if ( ~isvector(a) && ~isscalar(a))
+    error ("deconv: A must be vector");
+end
 
+  // Ensure A is oriented as Y.
+  if ((size(y,1)==1 && size(a,2)==1) || (size(y,2)==1 && size(a,1)==1))
+    a = a.';
+  end
+
+  la = length (a);
+  ly = length (y);
+
+  lb = ly - la + 1;
+
+  if (ly > la)
+    x = zeros (size (y,1) - size (a,1) + 1,size(y,2)-size(a,2)+1);
+    x(1) = 1;
+    [b, r] = filter (y, a, x);
+    r = r * a(1);
+  elseif (ly == la)
+    [b, r] = filter (y, a, 1);
+    r = r * a(1);
+  else
+    b = 0;
+    r = y;
+  end
+
+  if (nargout() > 1)
+    if (ly >= la)
+      r = [zeros(ly - la + 1, 1); r(1:la - 1)];
+      // Respect the orientation of Y
+      r = matrix (r, size (y,1),size(y,2));
+    end
+  end
+endfunction
 
+/*
+ //test
+ [b, r] = deconv ([3, 6, 9, 9], [1, 2, 3]);
+ assert_checkequal (b, [3, 0]);
+ assert_checkequal (r, [0, 0, 0, 9]);
 
+ //test
+ [b, r] = deconv ([3, 6], [1, 2, 3]);
+ assert_checkequal (b, 0);
+ assert_checkequal (r, [3, 6]);
 
+ //test
+ [b, r] = deconv ([3, 6], [1; 2; 3]);
+ assert_checkequal (b, 0);
+ assert_checkequal (r, [3, 6]);
 
-    [nargout,nargin]=argn();
+//test
+ [b,r] = deconv ([3; 6], [1; 2; 3]);
+ assert_checkequal (b, 0);
+ assert_checkequal (r, [3; 6]);
 
-    if (nargin ~= 2)
-        error ("wrong number of input arguments");
-    end
+//test
+ [b, r] = deconv ([3; 6], [1, 2, 3]);
+ assert_checkequal (b, 0);
+ assert_checkequal (r, [3; 6]);
 
-    if (~ (isvector (y) & isvector (a)))
-        error ("deconv: both arguments must be vectors");
-    end
+ assert_checkequal (deconv ((1:3)',[1, 1]), [1; 1])
 
-    la = length (a);
-    ly = length (y);
+// Test input validation
+// error deconv (1)
+// error deconv (1,2,3)
+// error <Y .* must be vector> deconv ([3, 6], [1, 2; 3, 4])
+// error <A must be vector> deconv ([3, 6], [1, 2; 3, 4])
 
-    lb = ly - la + 1;
+//test
+ y = (10:-1:1);
+ a = (4:-1:1);
+ [b, r] = deconv (y, a);
+ assert_checkequal (conv (a, b) + r, y);
 
-    // Ensure A is oriented as Y.
-    if (diff (size (y)) * diff (size (a)) < 0)
-        a = permute (a, [2, 1]);
-    end
 
-    if (ly > la)
-        o=size (y) - size (a) + 1;
-        x = zeros (o(1),o(2));
-        x(1) = 1;
-        b = filter (real(y), real(a), x);
-    elseif (ly == la)
-        b = filter (real(y), real(a), 1);
-    else
-        b = 0;
-    end
+//test
+ [b, r] = deconv ([1, 1], 1);
+ assert_checkequal (r, [0, 0]);
 
-    lc = la + length (b) - 1;
-    if (ly == lc)
-        if (length(a)==length(b) | length(a)>length(b))
-            if isrow(a)
-                q=conv(a,b);
-
-                u=[];
-                for i=1:length(y)
-                    u=[u;q];
-                end
-
-                w=[];
-                if (isrow(y))
-                for i=1:length(q)
-                    w=[w;y];
-                end
-            else
-                for i=1:length(q)
-                    w=[w,y];
-                    end
-            end
-
-                r = w-u;
-
-                r=diag(r);
-            end
-        end
-     r=y-conv(a,b);
-
-
-elseif(la~=lc)
-    // Respect the orientation of Y"
-    if (size (y,"r") <= size (y,"c"))
-        r = [(zeros (1, lc - ly)), y] - conv (a, b);
-    else
-        r = [(zeros (lc - ly, 1)); y] - conv (a, b);
-    end
-    if (ly < la)
-        // Trim the remainder is equal to the length of Y.
-        r = r($-(length(y)-1):$);
-    end
-end
+//test
+ [b, r] = deconv ([1; 1], 1);
+ assert_checkequal (r, [0; 0]);
 
-endfunction
+ */
diff --git a/macros/h1_z_deriv.sci b/macros/h1_z_deriv.sci
new file mode 100644
index 0000000..65e87e8
--- /dev/null
+++ b/macros/h1_z_deriv.sci
@@ -0,0 +1,47 @@
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+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  = [d 0];                           // Shift result right (multiply by -z)
+    d = d + prepad(polyder(d), i+1, 0, 2); // Add the derivative
+  end
+
+  // Build output vector
+  b = zeros (1, n + 1);
+  for i = (1:n)
+    b = b + d(i) * prepad(h1_deriv(n-i+1), n+1, 0, 2);
+  end
+  b  = b * ts^(n+1)/factorial(n);
+  // Multiply coefficients with p^i, where i is the index of the coeff.
+  b = 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,
+// 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)
+  b  = factorial(n)*nchoosek(n,0:n); // Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc.
+  b = b*(-1)^n;
+endfunction
+
+function y = polyder(a)
+    y = poly(flipdim(a,2),'a','coeff')
+    y = derivat(y)
+    y = coeff(y)
+endfunction
+
diff --git a/macros/impinvar.sci b/macros/impinvar.sci
index cc9aa6f..0e8b29b 100644
--- a/macros/impinvar.sci
+++ b/macros/impinvar.sci
@@ -1,43 +1,144 @@
-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)
+
+  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)
+
+*/
diff --git a/macros/impz.sci b/macros/impz.sci
index 1172b81..e4e1ba4 100644
--- a/macros/impz.sci
+++ b/macros/impz.sci
@@ -1,69 +1,125 @@
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+/*
+Calling Sequence
+ [x, t] = impz (b) ¶
+ [x, t] = impz (b, a) ¶
+ [x, t] = impz (b, a, n) ¶
+ [x, t] = impz (b, a, n, fs) ¶
+ impz (…) ¶
+Generate impulse-response characteristics of the filter. 
+The filter coefficients correspond to the the z-plane rational function with numerator b and denominator a. If a is not specified, it defaults to 1. 
+If n is not specified, or specified as [], it will be chosen such that the signal has a chance to die down to -120dB, or to not explode beyond 120dB, or to show five periods if there is no significant damping.
+If no return arguments are requested, plot the results.
+Dependencies
+  fftfilt
+*/
 function [x_r, t_r] = impz(b, a, n, fs)
-// It gives Impulse response of digital filter
 
-//Calling Sequence
-//x_r = impz(b)
-//x_r = impz(b, a)
-//x_r = impz(b, a, n)
-//x_r = impz(b, a, n, fs)
-//[x_r, t_r] = impz(b, a, n, fs)
+  if nargin < 1 || nargin > 4 then 
+    error(" impz : Incorrect number of input arguments ")
+  end
+  if nargin < 2 then a = [1]; end 
+  if nargin < 3 then n = [] ; end
+  if nargin < 4 then fs = 1 ; end
+  
+  if  (~isvector(b) && ~isscalar(b)) || (~isvector(a) && ~isscalar(a) )  then
+    error("impz: B and A must be vectors");
+  end
+  if isempty(n) && length(a) > 1 then 
+    precision = 1e-6;
+    r = roots(a);
+    maxpole = max(abs(r));
+    if (maxpole > 1+precision)     // unstable -- cutoff at 120 dB
+      n = floor(6/log10(maxpole));
+    elseif (maxpole < 1-precision) // stable -- cutoff at -120 dB
+      n = floor(-6/log10(maxpole));
+    else                           // periodic -- cutoff after 5 cycles
+      n = 30;
 
-//Parameters
-//x_r: impz chooses the number of samples and returns the response in the column vector, x_r.
-//t_r : impz returns the sample times in the column vector, t_r
-// b : numerator coefficients of the filter
-// a : denominator coefficients of the filter
-// n : samples of the impulse response   t(by default ,n = length(t) and  is computed automatically.
-// fs : sampling frequency
+      // find longest period less than infinity
+      // cutoff after 5 cycles (w=10*pi)
+      rperiodic = r(find(abs(r)>=1-precision & abs(arg(r))>0));
+      if ~isempty(rperiodic)
+        n_periodic = ceil(10*pi./min(abs(arg(rperiodic))));
+        if (n_periodic > n)
+          n = n_periodic;
+        end
+      end
 
-//Description
-//[x_r,t_r] = impz(b,a) returns the impulse response of the filter with numerator coefficients, b, and denominator coefficients, a. impz chooses the number of samples and returns the response in the column vector, x_r, and the sample times in the column vector, t_r. t_r = [0:n-1]' and n = length(t) is computed automatically.
+      // find most damped pole
+      // cutoff at -60 dB
+      rdamped = r(find(abs(r)<1-precision));
+      if ~isempty(rdamped)
+        n_damped = floor(-3/log10(max(abs(rdamped))));
+        if (n_damped > n)
+          n = n_damped;
+        end
+      end
+    end
+    n = n + length(b);
+  elseif isempty(n)
+    n = length(b);
+  elseif ( ~isscalar(n))
+    // TODO: missing option of having N as a vector of values to
+    //       compute the impulse response.
+    error ("impz: N must be empty or a scalar");
+  end
 
-//Examples
-//[x_r,t_r]=impz([0 1 1],[1 -3 3 -1],10)
-//OUTPUT :
-// t_r = 0.    1.     2.    3.    4.       5.      6.      7.      8.       9
-// x_r=  0.    1.    4.    9.    16.    25.    36.    49.....64......81
+  if length(a) == 1 then 
+    x = fftfilt(b/a, [1, zeros(1,n-1)]');
+  else
+    x = filter(b, a, [1, zeros(1,n-1)]');
+  end
+  t = [0:n-1]/fs;
 
-//[x_r,t_r]=impz(1,[1 1],5)
-//OUTPUT
-// t_r =   0.    1.    2.    3.    4
-//x_r =   1.  - 1.    1.  - 1.    1.
+  if nargout >= 1 x_r = x; end
+  if nargout >= 2 t_r = t; end
+  //if nargout ~= 2 then  //FIXME: fix nargout to detect 0 output arguments . till then plot it always
+      title("Impulse Response");
+      if (fs > 1000)
+        t = t * 1000;
+        xlabel("Time (msec)");
+      else
+        xlabel("Time (sec)");
+      end
+      plot(t, x,);
+  //end
 
-//This function is being called from Octave
+endfunction
+/*
 
+test case 1
+assert_checkequal(size(impz (1, [1 -1 0.9], 100)), [100 1]) // passed
 
-funcprot(0);
-rhs = argn(2)
-lhs = argn(1)
-if(rhs<1 | rhs>4)
-error("Wrong number of input arguments.")
-end
+test case 2
+// 7th order butterworth filter with fc = 0.5 //passed
+B=[0.016565   0.115957   0.347871   0.579785   0.579785   0.347871   0.115957   0.016565];
+A=[1.0000e+00  -5.6205e-16   9.1997e-01  -3.6350e-16   1.9270e-01  -4.3812e-17   7.6835e-03  -4.2652e-19];
+impz(B, A) 
 
-	select(rhs)
-	case 1 then
-	if(lhs==1)
-	[x_r] = callOctave("impz",b)
-	elseif(lhs==2)
-	[x_r,t_r] = callOctave("impz",b)
-	end
-	case 2 then
-	if(lhs==1)
-	[x_r] = callOctave("impz",b,a)
-	elseif(lhs==2)
-	[x_r,t_r] = callOctave("impz",b,a)
-	end
-	case 3 then
-	if(lhs==1)
-	[x_r] = callOctave("impz",b,a,n)
-	elseif(lhs==2)
-	[x_r,t_r] = callOctave("impz",b,a,n)
-	end
-	case 4 then
-	if(lhs==1)
-	[x_r] = callOctave("impz",b,a,n,fs)
-	elseif(lhs==2)
-	[x_r,t_r] = callOctave("impz",b,a,n,fs)
-	end
-	end
-endfunction
+test case 3
+// 
+[x_r,tr]=impz([0.4 0.2 8 2] ,1,10)
+assert_checkalmostequal(x_r,[ 0.4000000 0.2000000 8. 2.0000000 1.943D-16 1.388D-17 0. 0  0. 0.]',%eps,1e-4);
+assert_checkalmostequal(tr,0:9,%eps,1e-4);
+
+//test case 4 
+[xr,tr]=impz([0.4 0.2],[4 5 6],3,10) 
+
+
+// test case 5
+B = [0.021895   0.109474   0.218948   0.218948   0.109474   0.021895];
+A = [1.0000  -1.2210   1.7567  -1.3348   0.7556  -0.2560];
+impz(B, A) 
+
+ */
\ No newline at end of file
diff --git a/macros/inv_residue.sci b/macros/inv_residue.sci
new file mode 100644
index 0000000..e7354a6
--- /dev/null
+++ b/macros/inv_residue.sci
@@ -0,0 +1,54 @@
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+// Inverse of residue function
+function [b_out, a_out] = inv_residue(r_in, p_in, k_in, tol)
+
+  n = length(r_in); // Number of poles/residues
+
+  k = 0; // Capture constant 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");
+  end
+
+  a_out = octave_poly(p_in);
+
+  b_out  = zeros(1,n+1);
+  b_out = 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      = prepad(p, n+1, 0, 2);       // Pad for proper length
+    b_out = 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=i+1; // Next residue
+      m=m+1;
+      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, 2);           // Pad for proper length
+      b_out = b_out + p;
+    end
+  i=i+1;
+  end
+
+endfunction
+function ocpoly = octave_poly(A)
+    p = poly(A, 'x');
+    c = coeff(p);
+    ocpoly = c($:-1:1);
+endfunction
\ No newline at end of file
diff --git a/macros/invfreq.sci b/macros/invfreq.sci
index a720115..0fe4c81 100644
--- a/macros/invfreq.sci
+++ b/macros/invfreq.sci
@@ -1,43 +1,304 @@
-function [B,A] = invfreq(H,F,nB,nA,W,iter,tol, plane)
-// Calculates inverse frequency vectors
-//
-// Calling Sequence
-//[B,A] = invfreq(H,F,nB,nA)
-//[B,A] = invfreq(H,F,nB,nA,W)
-//[B,A] = invfreq(H,F,nB,nA,W,[],[],plane)
-//[B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane)
-// 
-// Parameters
-// H: desired complex frequency response,It is assumed that A and B are real polynomials, hence H is one-sided.
-// F: vector of frequency samples in radians
-// nA: order of denominator polynomial A
-// nB: order of numerator polynomial B
-//
-// Description
-//Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at frequency points F. A and B are real polynomial coefficients of order nA and nB respectively.  Optionally, the fit-errors can be weighted vs frequency according to the weights W. Also, the transform plane can be specified as either 's' for continuous time or 'z' for discrete time. 'z' is chosen by default.  Eventually, Steiglitz-McBride iterations will be specified by iter and tol.
-//
-// Examples
-//  [B,A] = butter(12,1/4);
-//  [H,w] = freqz(B,A,128);
-//  [Bh,Ah] = invfreq(H,F,4,4);
-//  Hh = freqz(Bh,Ah);
-//  disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));
-// 
-funcprot(0);
-lhs= argn(1);
-rhs= argn(2);
-if(rhs < 4 | rhs > 8 | rhs == 6 | rhs == 7 )
-error("Wrong number of input arguments");
-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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+/*
+  : [B,A] = invfreq(H,F,nB,nA,W)
+  : [B,A] = invfreq(H,F,nB,nA,W,[],[],plane) 
+  : [B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane) 
+
+    Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at frequency points F.
+
+    A and B are real polynomial coefficients of order nA and nB respectively. Optionally, the fit-errors can be weighted vs frequency according to the weights W. Also, the transform plane can be specified as either ’s’ for continuous time or ’z’ for discrete time. ’z’ is chosen by default. Eventually, Steiglitz-McBride iterations will be specified by iter and tol.
+
+    H: desired complex frequency response It is assumed that A and B are real polynomials, hence H is one-sided.
+
+    F: vector of frequency samples in radians
+
+    nA: order of denominator polynomial A
+
+    nB: order of numerator polynomial B
+
+    plane=’z’: F on unit circle (discrete-time spectra, z-plane design)
+
+    plane=’s’: F on jw axis (continuous-time spectra, s-plane design)
+
+    H(k) = spectral samples of filter frequency response at points zk, where zk=exp(sqrt(-1)*F(k)) when plane=’z’ (F(k) in [0,.5]) and zk=(sqrt(-1)*F(k)) when plane=’s’ (F(k) nonnegative)
+*/
+// FIXME: implement Steiglitz-McBride iterations
+// FIXME: improve numerical stability for high order filters (matlab is a bit better)
+// FIXME: modify to accept more argument configurations
+function [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, plane,varargin)
+
+  if nargin < 4  then 
+    error("invfreq : Incorrect number of input arguments ")
+  end
 
-select(rhs)
-	case 4 then
-		[B,A]= callOctave("invfreq", H,F,nB,nA);
-	case 5 then
-		[B,A]= callOctave("invfreq", H,F,nB,nA, W);
-	case 8 then
-		[B,A]= callOctave("invfreq", H,F,nB, nA,iter, tol, plane);
+  if ~isvector(H) && ~isscalar(H) then 
+    error("invfreq : H is the desired frequency response , a vector expected")
+  end
+
+  if ~isvector(F) && ~isscalar(F) then 
+    error("invfreq : F is a vector of frequency samples in radians")
+  end
+
+  if max(size(nB)) > 1 then zB = nB(2); nB = nB(1); else zB = 0; end
+  n = max(nA, nB);
+  m = n+1; mA = nA+1; mB = nB+1;
+  nF = max(size(F));
+  if nargin < 5 || isempty(W) then W = ones(1, nF); end
+  if nargin < 6 then iter = []; end
+  if nargin < 7  then tol = []; end
+  if nargin < 8 || isempty(tr) then  tr = ''; end
+  if nargin < 9 then plane = 'z'; end
+  if nargin < 10 then varargin = {}; end
+  if ( strcmp (plane, "s") &&  strcmp (plane, "z"))
+    error (sprintf("invfreq: invealid PLANE argument %s, expected  s  or  z ", plane))
+  end
+
+  fname = ["invfreq", plane];
+
+  if (nF ~= max(size(H))) then
+    error ("%s: Length of H and F must be the same\n", fname)
+  end
+
+  if (~ isempty (iter) || ~ isempty (tol)) then
+    warning (sprintf("%s: iterative algorithm not yet implemented, ", ...
+              "ITER and TOL arguments are ignored\n", fname));
+  end
+
+//////////////////////////////////////////////////////////////
+norm = %f ; // should we normalize freqs to avoid matrices with rank deficiency ?
+method = 'LS'; // by default, use Ordinary Least Square to solve normal equations
+prop = varargin;
+if length(prop) > 0 then 
+  indi = 1;
+  while indi < length(prop)
+      switch prop(indi)
+          case 'norm'
+              if indi < length(prop) && ~(type(prop(indi+1)) == 10)
+                  norm = prop(indi+1);
+                  indi = indi + 2; // Skip the processed element
+              else
+                  norm = %t; // Default true
+                  indi = indi + 1;
+              end
+              
+          case 'method'
+              if indi < length(prop) && type(prop(indi+1)) == 10 && strcmp(prop(indi+1), "norm")
+                  method = prop(indi+1);
+                  indi = indi + 2; // Skip the processed element
+              else
+                  error("invfreq : incorrect/missing method argument");
+                  indi = indi + 1;
+              end
+              
+          otherwise
+              disp("Ignoring unknown or incomplete argument");
+              indi = indi + 1;
+      end
+  end
 end
+
+
+////////////////////////////////////////////////////////////////
+
+
+  Ruu = zeros(mB, mB); Ryy = zeros(nA, nA); Ryu = zeros(nA, mB);
+  Pu = zeros(mB, 1);   Py = zeros(nA,1);
+  if ~strcmp(tr,'trace')
+    disp(' ')
+    disp('Computing nonuniformly sampled, equation-error, rational filter.');
+    disp(['plane = ',plane]);
+    disp(' ')
+  end
+
+  s = sqrt(-1)*F;
+  switch plane
+    case 'z' 
+      if max(F) > %pi || min(F) < 0 then
+        disp('hey, you frequency is outside the range 0 to %pi, making my own')
+        F = linspace(0, %pi, max(size(H)));
+        s = sqrt(-1)*F;
+      end
+      s = exp(-s);
+    case 's'
+      if max(F) > 1e6 && n > 5 then 
+        if ~norm then 
+          disp('Be careful, there are risks of generating singular matrices');
+          disp('Call invfreqs as (..., norm, 1) to avoid it');
+        else
+          Fmax = max(F); s = sqrt(-1)*F/Fmax;
+        end
+      end
+  end
+
+  //////////////////////////////
+  /////////////////////////////
+  for k=1:nF,
+    Zk = (s(k).^[0:n]).';
+    Hk = H(k);
+    aHks = Hk*conj(Hk);
+    Rk = (W(k)*Zk)*Zk';
+    rRk = real(Rk);
+    Ruu = clean(Ruu + rRk(1:mB, 1:mB));
+    Ryy = Ryy + aHks*rRk(2:mA, 2:mA);
+    Ryu = Ryu + real(Hk*Rk(2:mA, 1:mB));
+    Pu = Pu + W(k)*real(conj(Hk)*Zk(1:mB));
+    Py = Py + (W(k)*aHks)*real(Zk(2:mA));
+  end
+  Rr = ones(max(size(s)), mB+nA); Zk = s;
+  for k = 1:min(nA, nB),
+    Rr(:, 1+k) = Zk;
+    Rr(:, mB+k) = -Zk.*H;
+    Zk = Zk.*s;
+  end
+  for k = 1+min(nA, nB):max(nA, nB)-1,
+    if k <= nB, Rr(:, 1+k) = Zk; end
+    if k <= nA, Rr(:, mB+k) = -Zk.*H; end
+    Zk = Zk.*s;
+  end
+  k = k+1;
+  if k <= nB then Rr(:, 1+k) = Zk; end
+  if k <= nA then Rr(:, mB+k) = -Zk.*H; end
+
+  // complex to real equation system -- this ensures real solution
+  Rr = Rr(:, 1+zB:$);
+  Rr = [real(Rr); imag(Rr)]; Pr = [real(H(:)); imag(H(:))];
+  // normal equations -- keep for ref
+  // Rn= [Ruu(1+zB:mB, 1+zB:mB), -Ryu(:, 1+zB:mB)';  -Ryu(:, 1+zB:mB), Ryy];
+  // Pn= [Pu(1+zB:mB); -Py];
+  ////////////////////////////////////////////////
+  switch method
+    case {'ls' 'LS'}
+      // avoid scaling errors with Theta = R\P;
+      // [Q, R] = qr([Rn Pn]); Theta = R(1:$, 1:$-1)\R(1:$, $);
+      [Q, R] = qr([Rr Pr]); Theta = pinv(R(1:$-1, 1:$-1)) * R(1:$-1, $);
+      //////////////////////////////////////////////////
+      // SigN = R($, $-1);
+      SigN = R($, $);
+    case {'tls' 'TLS'}
+      // [U, S, V] = svd([Rn Pn]);
+      // SigN = S($, $-1);
+      // Theta =  -V(1:$-1, $)/V($, $);
+      [U, S, V] = svd([Rr Pr], 0);
+      SigN = S($, $);
+      Theta =  -V(1:$-1, $)/V($, $);
+    case {'mls' 'MLS' 'qr' 'QR'}
+      // [Q, R] = qr([Rn Pn], 0);
+      // solve the noised part -- DO NOT USE ECONOMY SIZE ~
+      // [U, S, V] = svd(R(nA+1:$, nA+1:$));
+      // SigN = S($, $-1);
+      // Theta = -V(1:$-1, $)/V($, $);
+      // unnoised part -- remove B contribution and back-substitute
+      // Theta = [R(1:nA, 1:nA)\(R(1:nA, $) - R(1:nA, nA+1:$-1)*Theta)
+      //         Theta];
+      // solve the noised part -- economy size OK as #rows > #columns
+      [Q, R] = qr([Rr Pr], 0);
+      eB = mB-zB; sA = eB+1;
+      [U, S, V] = svd(R(sA:$, sA:$));
+      // noised (A) coefficients
+      Theta = -V(1:$-1, $)/V($, $);
+      // unnoised (B) part -- remove A contribution and back-substitute
+      Theta = [R(1:eB, 1:eB)\(R(1:eB, $) - R(1:eB, sA:$-1)*Theta)
+               Theta];
+      SigN = S($, $);
+    otherwise
+      error(sprintf("invfreq : unknown method %s", method));
+  end
+
+  B = [zeros(zB, 1); Theta(1:mB-zB)].';
+  A = [1; Theta(mB-zB+(1:nA))].';
+  if ~strcmp(plane,'s')
+    B = B(mB:-1:1);
+    A = A(mA:-1:1);
+    if norm, // Frequencies were normalized -- unscale coefficients
+      Zk = Fmax.^[n:-1:0].';
+      for k = nB:-1:1+zB, B(k) = B(k)/Zk(k); end
+      for k = nA:-1:1, A(k) = A(k)/Zk(k); end
+    end
+  end
 endfunction
 
- 
+/*
+// method - LS
+test case 1 // passed
+
+[B,A,Sign] = invfreq(1,1,1,1,1,[],[],'','z','norm',1,'method','LS')
+assert_checkequal(B,[0.6314 0.3411])
+assert_checkequal(A,[1 -0.3411])
+assert_checkequal(Sign,0)
+
+[B,A,Sign] = invfreq(1,1,1,1,1,[],[],'','s')
+assert_checkequal(B,[0 1])
+assert_checkequal(A,[0 1])
+assert_checkequal(Sign,0)
+
+
+test case 2 // passed 
+order = 6 
+fc = 1/2 
+n = 128 
+B = [    0.029588   0.177529   0.443823   0.591764   0.443823   0.177529   0.029588] ;
+A = [ 1.0000e+00  -6.6613e-16   7.7770e-01  -2.8192e-16   1.1420e-01  -1.4472e-17   1.7509e-03];
+[H,w] = freqz(B,A,n) ; 
+[Bh , Ah] = invfreq(H,w,order,order);
+[Hh,wh] = freqz(Bh,Ah,n);
+plot(w,[abs(H), abs(Hh)])
+xlabel("Frequency (rad/sample)");
+ylabel("Magnitude");
+legend('Original','Measured');
+err = norm(H-Hh);
+disp(sprintf('L2 norm of frequency response error = %f',err));
+
+test case 3 // passed 
+// buttter worth filter of order 12 and fc=1/4
+B = [ 1.1318e-06   1.3582e-05   7.4702e-05   2.4901e-04   5.6026e-04   8.9642e-04   1.0458e-03   8.9642e-04 5.6026e-04   2.4901e-04   7.4702e-05   1.3582e-05   1.1318e-06];
+A = [ 1.0000e+00  -5.9891e+00   1.7337e+01  -3.1687e+01   4.0439e+01  -3.7776e+01   2.6390e+01  -1.3851e+01  5.4089e+00  -1.5296e+00   2.9688e-01  -3.5459e-02   1.9688e-03];
+[H,w] = freqz(B,A,128);
+[Bh,Ah] = invfreq(H,w,4,4);
+[Hh,wh] = freqz(Bh,Ah,128);
+disp(sprintf('||frequency response error||= %f',norm(H-Hh)));
+
+
+method TLS 
+
+test case 1 // passed
+
+B = [ 1.1318e-06   1.3582e-05   7.4702e-05   2.4901e-04   5.6026e-04   8.9642e-04   1.0458e-03   8.9642e-04 5.6026e-04   2.4901e-04   7.4702e-05   1.3582e-05   1.1318e-06];
+A = [ 1.0000e+00  -5.9891e+00   1.7337e+01  -3.1687e+01   4.0439e+01  -3.7776e+01   2.6390e+01  -1.3851e+01  5.4089e+00  -1.5296e+00   2.9688e-01  -3.5459e-02   1.9688e-03];
+[H,w] = freqz(B,A,128);
+[Bh,Ah] = invfreq(H,w,4,4,[],[],[],'','z','norm',1,'method','TLS');
+[Hh,wh] = freqz(Bh,Ah,128);
+disp(sprintf('||frequency response error||= %f',norm(H-Hh)));
+
+
+
+
+method MLS - passed 
+
+
+// elliptic filter with  ellip (5, 1, 90, [.1 .2])
+n = 128 
+B = [ 1.3214e-04  -6.6404e-04   1.4928e-03  -1.9628e-03   1.4428e-03  0  -1.4428e-03   1.9628e-03 -1.4928e-03   6.6404e-04  -1.3214e-04] ;
+A = [ 1.0000    -8.6483    34.6032   -84.2155   137.9276  -158.7598   130.0425   -74.8636    29.0044    -6.8359  0.7456];
+[H,w] = freqz(B,A,n) ; 
+
+[Bh,Ah] = invfreq(H,w,4,4,[],[],[],'','z','norm',1,'method','MLS');
+[Hh,wh] = freqz(Bh,Ah,n);
+plot(w,[abs(H), abs(Hh)])
+xlabel("Frequency (rad/sample)");
+ylabel("Magnitude");
+legend('Original','Measured');
+err = norm(H-Hh);
+disp(sprintf('L2 norm of frequency response error = %f',err));
+
+
+*/
diff --git a/macros/invfreqs.sci b/macros/invfreqs.sci
index 32e732b..99c721b 100644
--- a/macros/invfreqs.sci
+++ b/macros/invfreqs.sci
@@ -1,96 +1,110 @@
-function [B,A,C] = invfreqs(H,F,nB,nA,W,iter,tol,trace)
-//Fit filter B(s)/A(s)to the complex frequency response H at frequency points F.  A and B are real polynomial coefficients of order nA and nB.
-//Calling Sequence
-//[B,A,C] = invfreqs(H,F,nB,nA,W,iter,tol,trace)
-//[B,A,C] = invfreqs(H,F,nB,nA,W)
-//[B,A,C] = invfreqs(H,F,nB,nA)
-//Parameters
-//H: desired complex frequency response.
-//F: frequency (must be same length as H).
-//nB: order of the numerator polynomial B.
-//nA: order of the denominator polynomial A.
-//W: vector of weights (must be same length as F).
-//Description
-//This is an Octave function.
-//Fit filter B(s)/A(s)to the complex frequency response H at frequency points F.  A and B are real polynomial coefficients of order nA and nB.
-//Optionally, the fit-errors can be weighted vs frequency according to the weights W.
-//Note: all the guts are in invfreq.m 
-//Examples
-//B = [1/2 1];
-//A = [1 1];
-//w = linspace(0,4,128);
-//H = freqs(B,A,w);
-//[Bh,Ah, C] = invfreqs(H,w,1,1);
-//Bh =
-//
-//   0.50000   1.00000
-//
-//Ah =
-//
-//   1.0000   1.0000
-//
-//C =   -3.0964e-15
-
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 4 | rhs > 8)
-error("Wrong number of input arguments.")
-end
-
-select(rhs)
-	
-	case 4 then
-		if(lhs==1)
-		B = callOctave("invfreqs",H,F,nB,nA)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqs",H,F,nB,nA)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqs",H,F,nB,nA)
-		else
-		error("Wrong number of output argments.")
-		end
-
-	case 5 then
-		if(lhs==1)
-		B = callOctave("invfreqs",H,F,nB,nA,W)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqs",H,F,nB,nA,W)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqs",H,F,nB,nA,W)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 6 then
-		if(lhs==1)
-		B = callOctave("invfreqs",H,F,nB,nA,W,iter)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 7 then
-		if(lhs==1)
-		B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 8 then
-		if(lhs==1)
-		B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace)
-		else
-		error("Wrong number of output argments.")
-		end
-	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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+// FIXME: check invfreq.sci for todo's
+
+/*
+  :[B,A] = invfreqs(H,F,nB,nA)
+  :[B,A] = invfreqs(H,F,nB,nA,W) 
+  :[B,A] = invfreqs(H,F,nB,nA,W,iter,tol,'trace') 
+
+    Fit filter B(s)/A(s)to the complex frequency response H at frequency points F.
+
+    A and B are real polynomial coefficients of order nA and nB.
+
+    Optionally, the fit-errors can be weighted vs frequency according to the weights W.
+
+    Note: all the guts are in invfreq.m
+
+    H: desired complex frequency response
+
+    F: frequency (must be same length as H)
+
+    nA: order of the denominator polynomial A
+
+    nB: order of the numerator polynomial B
+
+    W: vector of weights (must be same length as F)
+
+*/
+// Dependencies
+// invfreq
+function [B, A, SigN] = invfreqs(H,F,nB,nA,W,iter,tol,tr, varargin)
+
+  if nargin < 9
+    varargin = {};
+    if nargin < 8
+      tr = '';
+      if nargin < 7
+        tol = [];
+        if nargin < 6
+          iter = [];
+          if nargin < 5
+            W = ones(1,length(F));
+          end
+        end
+      end
+    end
+  end
+
+  // now for the real work
+  [B, A, SigN] = invfreq(H, F,nB, nA, W, iter, tol, tr, 's', varargin);
+
 endfunction
+/*
+demo 
+ B = [1 0 0];
+ A = [1 6 15 15]/15;
+ w = linspace(0, 8, 128);
+ [H0 ,_ ] = freqz(B, A, w);
+ Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2);
+ order = length(A) - 1;
+ [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1);
+ [Hh ,_ ] = freqz(Bh,Ah,w);
+ [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS");
+ [HLS,_] = freqz(BLS, ALS, w);
+ [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS");
+ [HTLS,_] = freqz(BTLS, ATLS, w);
+ [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR");
+ [HMLS,_] = freqz(BMLS, AMLS, w);
+ plot(w,[abs(H0); abs(Hh)])
+ xlabel("Frequency (rad/sec)");
+ ylabel("Magnitude");
+ legend('Original','Measured');
+ err = norm(H0-Hh);
+ disp(sprintf('L2 norm of frequency response error = %f',err));
+*/
+/*  Octave version
+ B = [1 0 0];
+ A = [1 6 15 15]/15;
+ w = linspace(0, 8, 128);
+ [H0 ,_ ] = freqz(B, A, w);
+ Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2);
+ order = length(A) - 1;
+ [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1);
+ [Hh ,_ ] = freqz(Bh,Ah,w);
+ [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS");
+ [HLS,_] = freqz(BLS, ALS, w);
+ [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS");
+ [HTLS,_] = freqz(BTLS, ATLS, w);
+ [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR");
+ [HMLS,_] = freqz(BMLS, AMLS, w);
+ plot(w,[abs(H0); abs(Hh)])
+ xlabel("Frequency (rad/sec)");
+ ylabel("Magnitude");
+ legend('Original','Measured');
+ err = norm(H0-Hh);
+ disp(sprintf('L2 norm of frequency response error = %f',err));
+*/
+
+
+
 
diff --git a/macros/invfreqz.sci b/macros/invfreqz.sci
index d051499..433100b 100644
--- a/macros/invfreqz.sci
+++ b/macros/invfreqz.sci
@@ -1,94 +1,107 @@
-function [B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace)
-//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F.  A and B are real polynomial coefficients of order nA and nB.
-//Calling Sequence
-//[B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace)
-//[B,A,C] = invfreqz(H,F,nB,nA,W)
-//[B,A,C] = invfreqz(H,F,nB,nA)
-//Parameters
-//H: desired complex frequency response.
-//F: frequency (must be same length as H).
-//nB: order of the numerator polynomial B.
-//nA: order of the denominator polynomial A.
-//W: vector of weights (must be same length as F).
-//Description
-//This is an Octave function.
-//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F.  A and B are real polynomial coefficients of order nA and nB.
-//Optionally, the fit-errors can be weighted vs frequency according to the weights W.
-//Note: all the guts are in invfreq.m 
-//Examples
-//[B,A] = butter(4,1/4);
-//[H,F] = freqz(B,A);
-//[Bh,Ah,C] = invfreq(H,F,4,4)
-//Bh =
-//
-//   0.010209   0.040838   0.061257   0.040838   0.010209
-//
-//Ah =
-//
-//   1.00000  -1.96843   1.73586  -0.72447   0.12039
-//
-//C =   -7.7065e-15
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+// FIXME: check invfreq.sci for todo's
+/*
+  : [B,A] = invfreqz(H,F,nB,nA) ¶
+  : [B,A] = invfreqz(H,F,nB,nA,W) ¶
+  : [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,'trace') ¶
 
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 4 | rhs > 8)
-error("Wrong number of input arguments.")
-end
+  Fit filter B(z)/A(z)to the complex frequency response H at frequency points F.
 
-select(rhs)
-	
-	case 4 then
-		if(lhs==1)
-		B = callOctave("invfreqz",H,F,nB,nA)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqz",H,F,nB,nA)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqz",H,F,nB,nA)
-		else
-		error("Wrong number of output argments.")
-		end
+  A and B are real polynomial coefficients of order nA and nB. Optionally, the fit-errors can be weighted vs frequency according to the weights W.
 
-	case 5 then
-		if(lhs==1)
-		B = callOctave("invfreqz",H,F,nB,nA,W)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqz",H,F,nB,nA,W)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqz",H,F,nB,nA,W)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 6 then
-		if(lhs==1)
-		B = callOctave("invfreqz",H,F,nB,nA,W,iter)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 7 then
-		if(lhs==1)
-		B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol)
-		else
-		error("Wrong number of output argments.")
-		end
-	case 8 then
-		if(lhs==1)
-		B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace)
-		elseif(lhs==2)
-		[B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace)
-		elseif(lhs==3)
-		[B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace)
-		else
-		error("Wrong number of output argments.")
-		end
-	end
-endfunction
+  Note: all the guts are in invfreq.m
+
+  H: desired complex frequency response
+
+  F: normalized frequency (0 to pi) (must be same length as H)
+
+  nA: order of the denominator polynomial A
+
+  nB: order of the numerator polynomial B
 
+  W: vector of weights (must be same length as F)
+
+*/
+// Dependencies
+// invfreq
+function [B, A, SigN] = invfreqz(H, F, nB, nA, W, iter, tol, tr, varargin)
+
+  if nargin < 9
+    varargin = {};
+    if nargin < 8
+      tr = '';
+      if nargin < 7
+        tol = [];
+        if nargin < 6
+          iter = [];
+          if nargin < 5
+            W = ones(1,length(F));
+          end
+        end
+      end
+    end
+  end
+  // now for the real work
+  [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, 'z', varargin);
+
+endfunction
+/*
+demo       
+order = 9;  //order of test filter
+//going to 10 or above leads to numerical instabilities and large errors
+fc = 1/2;   // sampling rate / 4
+n = 128;    // frequency grid size
+// butterworth filter of order 9 and fc=0.5
+B0 = [5.1819e-03   4.6637e-02   1.8655e-01   4.3528e-01   6.5292e-01   6.5292e-01   4.3528e-01   1.8655e-01  4.6637e-02   5.1819e-03];
+A0 = [  1.0000e+00  -8.6736e-16   1.2010e+00  -7.7041e-16   4.0850e-01  -1.7013e-16   4.2661e-02  -9.0155e-18 9.6666e-04  -5.3661e-20];
+[H0, w] = freqz(B0, A0, n);
+Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2);
+[Bh, Ah, Sig0] = invfreqz(H0, w, order, order);
+[Hh, wh] = freqz(Bh, Ah, n);
+[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS");
+[HLS _ ] = freqz(BLS, ALS, n);
+[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS");
+[HTLS _ ]= freqz(BTLS, ATLS, n);
+[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR");
+[HMLS _ ] = freqz(BMLS, AMLS, n);
+plot(w,[abs(H0) abs(Hh)])
+xlabel("Frequency (rad/sample)");
+ylabel("Magnitude");
+legend('Original','Measured');
+err = norm(H0-Hh);
+disp(sprintf('L2 norm of frequency response error = %f',err));
+        
+*/
+/*
+order = 9; 
+fc = 1/2; 
+n = 128;
+B0 = [5.1819e-03   4.6637e-02   1.8655e-01   4.3528e-01   6.5292e-01   6.5292e-01   4.3528e-01   1.8655e-01  4.6637e-02   5.1819e-03];
+A0 = [  1.0000e+00  -8.6736e-16   1.2010e+00  -7.7041e-16   4.0850e-01  -1.7013e-16   4.2661e-02  -9.0155e-18 9.6666e-04  -5.3661e-20];
+[H0, w] = freqz(B0, A0, n);
+Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2);
+[Bh, Ah, Sig0] = invfreqz(H0, w, order, order);
+[Hh, wh] = freqz(Bh, Ah, n);
+[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS");
+[HLS _ ] = freqz(BLS, ALS, n);
+[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS");
+[HTLS _ ]= freqz(BTLS, ATLS, n);
+[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR");
+[HMLS _ ] = freqz(BMLS, AMLS, n);
+plot(w,[abs(H0) abs(Hh)])
+xlabel("Frequency (rad/sample)");
+ylabel("Magnitude");
+legend('Original','Measured');
+err = norm(H0-Hh);
+disp(sprintf('L2 norm of frequency response error = %f',err));
+*/
diff --git a/macros/invimpinvar.sci b/macros/invimpinvar.sci
index 4e8bd7b..b94a8b0 100644
--- a/macros/invimpinvar.sci
+++ b/macros/invimpinvar.sci
@@ -1,42 +1,125 @@
-function [b_out, a_out] = invimpinvar (b, a, fs, tol)
-//This function converts digital filter with coefficients b and a to analog, 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 digital filter with coefficients b and a to analog, 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] = invimpinvar(b,a,10)
-//ay =
-//  -1.6940e-16   4.6223e+00  -4.5210e+00   7.2880e+02
-//by =
-// Columns 1 through 4:
-//   1.0000e+00   3.0900e+01   9.6532e+02   1.2232e+04
-// Column 5:
-//   1.1038e+05
-funcprot(0);
-rhs = argn(2)
-if(rhs<2)
-error("Wrong number of input arguments.")
-end
-
-
-	select(rhs)
-	case 2 then
-	[b_out,a_out] = callOctave("invimpinvar",b,a)
-	case 3 then
-	[b_out,a_out] = callOctave("invimpinvar",b,a,fs)
-	case 4 then
-	[b_out,a_out] = callOctave("invimpinvar",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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+
+// Impulse invariant conversion from s to z domain
+/*
+    [b_out, a_out] = invimpinvar (b, a, fs, tol) 
+    [b_out, a_out] = invimpinvar (b, a, fs) 
+    [b_out, a_out] = invimpinvar (b, a) 
+
+    Converts digital filter with coefficients b and a to analog, conserving impulse response.
+
+    This function does the inverse of impinvar so that the following example should restore the original values of a and b.
+
+    [b, a] = impinvar (b, a);
+    [b, a] = invimpinvar (b, a);
+
+    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%) 
+  Dependencies
+    residue
+    inv_residue
+    */
+function [b_out, a_out] = invimpinvar (b_in, a_in, fs, tol)
+
+  if (nargin <2)
+    error("invimpinvar: Insufficient input arguments");
+  end
+
+  if nargin < 3 then fs = 1; end
+  if nargin < 4 then 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
+
+  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
+  
+  // clean r_in for zero values
+  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");
+  end
+
+  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=i+1; // Next residue
+      m=m+1; // Next multiplicity
+    end
+    [r, sm, k]= inv_z_res(r_in(i-m+1:i), first_pole, ts); // Find s-domain residues
+    k_in            = 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=i+1; // Next z-domain residue/pole
+  end
+  [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);
+
+  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
+
+  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_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=j-1;
+  end
+
+  // 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
+
+/*
+// tests passed
+[b_out,a_out]=invimpinvar([1],[1 -0.5],0.01) 
+[b_out,a_out]=invimpinvar([1],[1 -1 0.25],0.01)
+[b_out,a_out]=invimpinvar([1 1],[1 -1 0.25],0.01)
+[b_out,a_out]=invimpinvar([1],[1 -1.5 0.75 -0.125],0.01)
+[b_out,a_out]=invimpinvar([1 1],[1 -1.5 0.75 -0.125],0.01)
+
+
+
+// FIXME : built in filter doesn't support complex parameters
+// Because of this thsese test cases are failing
+//[b_out,a_out]=invimpinvar([1],[1 0 0.25],0.01)
+// [b_out,a_out]=invimpinvar([1 1],[1 0 0.25],0.01)
+// [b_out,a_out]=invimpinvar([1],[1 0 0.5 0 0.0625],0.01)
+// [b_out,a_out]=invimpinvar([1 1],[1 0 0.5 0 0.0625],0.01)
+// [b_out,a_out]=invimpinvar([1 1 1],[1 0 0.5 0 0.0625],0.01
+// [b_out,a_out]=invimpinvar([1 1 1 1],[1 0 0.5 0 0.0625],0.01)
+
+*/
diff --git a/macros/marcumq.sci b/macros/marcumq.sci
index 630ec97..29f4416 100644
--- a/macros/marcumq.sci
+++ b/macros/marcumq.sci
@@ -1,38 +1,389 @@
+// 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/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+/*
+Calling Sequence 
+ q = marcumq (a, b)
+ q = marcumq (a, b, m)
+ q = marcumq (a, b, m, tol)
+ 
+Input and Output parameters
+a — Noncentrality parameter --- nonnegative scalar | array of nonnegative numbers 
+b — Argument of Marcum Q-function --- nonnegative scalar | array of nonnegative numbers
+m — Order of generalized Marcum Q-function --- positive integer | array of positive integers
+
+Compute the generalized Marcum Q function of order `m` with noncentrality parameter `a` and argument `b`. 
+If the order `m` is omitted, it defaults to 1. An optional relative tolerance `tol` may be included, 
+and the default value is `eps`.
+
+If the input arguments are commensurate vectors, this function will produce a table of values.
+
+This function computes Marcum’s Q function using the infinite Bessel series, 
+which is truncated when the relative error is less than the specified tolerance. 
+The accuracy is limited by that of the Bessel functions, so reducing the tolerance is probably not useful.
+
+References:
+- Marcum, "Tables of Q Functions", Rand Corporation.
+- R.T. Short, "Computation of Noncentral Chi-squared and Rice Random Variables", 
+  www.phaselockedsystems.com/publications
+*/
+
 function q = marcumq (a, b, m, tol)
-//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b. 
-//Calling Sequence
-//q = marcumq (a, b)
-//q = marcumq (a, b, m)
-//q = marcumq (a, b, m, tol)
-//Parameters 
-//a:
-//b:
-//m: default value 1
-//tol: default value eps
-//Description
-//This is an Octave function.
-//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b. 
-//The third argument m is the order, which by default is 1.
-//The fourth argument tol is the tolerance, which by default is eps.
-//If input arguments are vectors which correspond in size and degree, the output is a table of values.
-//This function calculates Marcum’s Q function using the infinite Bessel series, which is truncated when the relative error is less than the specified tolerance.
-//Examples
-//marcumq([1,2,3],4)
-//ans  =
-//    0.0028895    0.0341348    0.1965122 
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<2 | rhs>4)
-error("Wrong number of input arguments.")
-end
-	select(rhs)
-	case 2 then
-	q = callOctave("marcumq",a,b)
-	case 3 then
-	q = callOctave("marcumq",a,b,m)
-	case 4 then
-	q = callOctave("marcumq",a,b,m,tol)
-	end
+
+  if ((nargin < 2) || (nargin > 5))
+    error(" marcumq : wrong numbers of input arguments ");
+  end
+
+  if nargin < 3 then m = 1 end
+  if nargin < 4 then tol = %eps end
+  if nargin < 5 then max_iter = 100 end 
+
+  if (or (a < 0))
+    error ("marcumq: A must be a non-negative value");
+  end
+  if (or (b < 0))
+    error ("marcumq: B must be a non-negative value");
+  end
+  if (or (m < 1) || or (fix (m) ~= m))
+    error ("marcumq: M must be a positive integer");
+  end
+
+  [a, b] = tablify (a, b);
+  q = []
+  for i=1:size(a,1)
+    for j=1:size(a,2)
+      q(i,j) = mq(a(i,j),b(i,j),m,tol)
+    end
+  end
+
 endfunction
-	
+
+// Subfunction to compute the actual Marcum Q function.
+function q = mq (a, b, m, tol)
+
+  // Special cases.
+  if (b == 0)
+    q = 1;
+    N = 0;
+    return;
+  end
+  if (a == 0)
+    k = 0:(m - 1);
+    q = exp (-b^2 / 2) * sum (b.^(2 * k) ./ (2.^k .* factorial (k)));
+    N = 0;
+    return;
+  end
+
+  // The basic iteration.  If a<b compute Q_M, otherwise
+  // compute 1-Q_M.
+  k = m;
+  z = a * b;
+  t = 1;
+  k = 0;
+  if (a < b)
+    s = 1;
+    c = 0;
+    x = a / b;
+    d = x;
+    S = besseli (0, z, 1);
+    if (m > 1)
+      for k = 1:m - 1
+        t = (d + 1 / d) * besseli (k, z, 1);
+        S = S + t;
+        d = d * x;
+      end
+    end
+    N = k;
+    k = k + 1
+  else
+    s = -1;
+    c =  1;
+    x = b / a;
+    k = m;
+    d = x^m;
+    S = 0;
+    N = 0;
+  end
+
+  while  (abs (t / S) > tol)
+    t = d * besseli (abs (k), z, 1);
+    S = S + t;
+    d = d * x;
+    N = k;
+    k = k + 1 ;
+  end
+  q = c + s * exp (-(a - b)^2 / 2) * S;
+
+endfunction
+
+// Internal helper function to create a table of like dimensions from arguments.
+function [ta , tb] = tablify(a,b)
+  rows = size(a,1)
+  cols = size(b,2)
+  ta=[]
+  for i=1:cols
+    ta = [ta a ]
+  end
+  tb=[]
+  for i=1:rows
+    tb = [ tb ; b]
+  end
+endfunction
+/*
+test
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [0.000000, 0.100000, 1.100000, 2.100000, 3.100000, 4.100000];
+ Q = [1.000000, 0.995012, 0.546074, 0.110251, 0.008189, 0.000224;
+      1.000000, 0.995019, 0.546487, 0.110554, 0.008238, 0.000226;
+      1.000000, 0.996971, 0.685377, 0.233113, 0.034727, 0.002092;
+      1.000000, 0.999322, 0.898073, 0.561704, 0.185328, 0.027068;
+      1.000000, 0.999944, 0.985457, 0.865241, 0.526735, 0.169515;
+      1.000000, 0.999998, 0.999136, 0.980933, 0.851679, 0.509876;
+      1.000000, 1.000000, 0.999979, 0.998864, 0.978683, 0.844038;
+      1.000000, 1.000000, 1.000000, 0.999973, 0.998715, 0.977300;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999969, 0.998618;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+q = marcumq (a, b);
+assert_checkalmostequal (q, Q, %eps,1e-4);
+
+test 
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [5.100000, 6.100000, 7.100000, 8.100000, 9.100000, 10.10000];
+ Q = [0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000049, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.001606, 0.000037, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.024285, 0.001420, 0.000032, 0.000000, 0.000000, 0.000000;
+      0.161412, 0.022812, 0.001319, 0.000030, 0.000000, 0.000000;
+      0.499869, 0.156458, 0.021893, 0.001256, 0.000028, 0.000000;
+      0.839108, 0.493229, 0.153110, 0.021264, 0.001212, 0.000027;
+      0.976358, 0.835657, 0.488497, 0.150693, 0.020806, 0.001180;
+      0.998549, 0.975673, 0.833104, 0.484953, 0.148867, 0.020458;
+      0.999965, 0.998498, 0.975152, 0.831138, 0.482198, 0.147437;
+      1.000000, 0.999963, 0.998458, 0.974742, 0.829576, 0.479995;
+      1.000000, 1.000000, 0.999962, 0.998426, 0.974411, 0.828307;
+      1.000000, 1.000000, 1.000000, 0.999961, 0.998400, 0.974138;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999960, 0.998378;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999960;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b);
+ assert_checkalmostequal(q, Q, %eps,1e-4);
+
+
+test 
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [11.10000, 12.10000, 13.10000, 14.10000, 15.10000, 16.10000];
+ Q = [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.000026, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.001155, 0.000026, 0.000000, 0.000000, 0.000000, 0.000000;
+      0.020183, 0.001136, 0.000025, 0.000000, 0.000000, 0.000000;
+      0.146287, 0.019961, 0.001120, 0.000025, 0.000000, 0.000000;
+      0.478193, 0.145342, 0.019778, 0.001107, 0.000024, 0.000000;
+      0.827253, 0.476692, 0.144551, 0.019625, 0.001096, 0.000024;
+      0.973909, 0.826366, 0.475422, 0.143881, 0.019494, 0.001087;
+      0.998359, 0.973714, 0.825607, 0.474333, 0.143304, 0.019381;
+      0.999959, 0.998343, 0.973546, 0.824952, 0.473389, 0.142803;
+      1.000000, 0.999959, 0.998330, 0.973400, 0.824380, 0.472564;
+      1.000000, 1.000000, 0.999958, 0.998318, 0.973271, 0.823876;
+      1.000000, 1.000000, 1.000000, 0.999958, 0.998307, 0.973158;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999957, 0.998297;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999957;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b);
+ assert_checkalmostequal (q, Q,%eps,1e-4);
+
+
+test 
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [17.10000, 18.10000, 19.10000];
+ Q = [0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000000, 0.000000, 0.000000;
+      0.000024, 0.000000, 0.000000;
+      0.001078, 0.000024, 0.000000;
+      0.019283, 0.001071, 0.000023;
+      0.142364, 0.019197, 0.001065;
+      0.471835, 0.141976, 0.019121;
+      0.823429, 0.471188, 0.141630;
+      0.973056, 0.823030, 0.470608;
+      0.998289, 0.972965, 0.822671;
+      0.999957, 0.998281, 0.972883;
+      1.000000, 0.999957, 0.998274;
+      1.000000, 1.000000, 0.999956;
+      1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b);
+ assert_checkalmostequal(q, Q, %eps,1e-4);
+
+// The tests for M>1 were generating from Marcum's tables by
+// using the formula
+//   Q_M(a,b) = Q(a,b) + exp(-(a-b)^2/2)*sum_{k=1}^{M-1}(b/a)^k*exp(-ab)*I_k(ab)
+
+test 
+ M = 2;
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
+ Q = [1.000000, 0.999987, 0.353353, 0.000000, 0.000000, 0.000000;
+      1.000000, 0.999988, 0.353687, 0.000000, 0.000000, 0.000000;
+      1.000000, 0.999992, 0.478229, 0.000000, 0.000000, 0.000000;
+      1.000000, 0.999999, 0.745094, 0.000001, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.934771, 0.000077, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.992266, 0.002393, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999607, 0.032264, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999992, 0.192257, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.545174, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.864230, 0.000040, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.981589, 0.001555, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.998957, 0.024784, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.999976, 0.166055, 0.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.509823, 0.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.846066, 0.000032;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.978062, 0.001335;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.998699, 0.022409;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999970, 0.156421;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.495223;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.837820;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.976328;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.998564;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b, M);
+ assert_checkalmostequal (q, Q, %eps,1e-4);
+
+test 
+
+ M = 5;
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
+ Q = [1.000000, 1.000000, 0.926962, 0.000000, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.927021, 0.000000, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.947475, 0.000001, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.980857, 0.000033, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.996633, 0.000800, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999729, 0.011720, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999990, 0.088999, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.341096, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.705475, 0.000002, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.933009, 0.000134, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.993118, 0.003793, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.999702, 0.045408, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.999995, 0.238953, 0.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.607903, 0.000001;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.896007, 0.000073;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.987642, 0.002480;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999389, 0.034450;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999988, 0.203879;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.565165;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.876284;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.984209;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999165;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999983;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b, M);
+ assert_checkalmostequal (q, Q, %eps,1e-4);
+
+test passed
+ M = 10;
+ a = [0.00; 0.05; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00; 7.00; 8.00; 9.00; 10.00;
+      11.00; 12.00; 13.00; 14.00; 15.00; 16.00; 17.00; 18.00; 19.00; 20.00;
+      21.00; 22.00; 23.00; 24.00];
+ b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
+ Q = [1.000000, 1.000000, 0.999898, 0.000193, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999897, 0.000194, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999931, 0.000416, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999980, 0.002377, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999997, 0.016409, 0.000000, 0.000000;
+      1.000000, 1.000000, 0.999999, 0.088005, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.302521, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.638401, 0.000000, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.894322, 0.000022, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.984732, 0.000840, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.998997, 0.014160, 0.000000;
+      1.000000, 1.000000, 1.000000, 0.999972, 0.107999, 0.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.391181, 0.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.754631, 0.000004;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.951354, 0.000266;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.995732, 0.006444;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999843, 0.065902;
+      1.000000, 1.000000, 1.000000, 1.000000, 0.999998, 0.299616;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.676336;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.925312;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.992390;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999679;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999995;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+ q = marcumq (a, b, M);
+ assert_checkalmostequal (q, Q, %eps,1e-4);
+
+*/
\ No newline at end of file
diff --git a/macros/mpoles.sci b/macros/mpoles.sci
index bba9997..4556120 100644
--- a/macros/mpoles.sci
+++ b/macros/mpoles.sci
@@ -1,124 +1,128 @@
 // 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
-// Author:[insert name]
+// Original Source : https://octave.sourceforge.io/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
 // Organization: FOSSEE, IIT Bombay
 // Email: toolbox@scilab.in
+function [multp, idxp] = mpoles (p, tol, reorder)
 
+  if (nargin < 1)
+    error("mpoles: Invalid number of input arguments");
+  end
 
-function [multp, indx] = mpoles (p, tol, reorder)
-
-//Calling sequence:
-// [multp, idxp] = mpoles (p)
-// [multp, idxp] = mpoles (p, tol)
-//  [multp, idxp] = mpoles (p, tol, reorder)
-// Identify unique poles in p and their associated multiplicity.
-//
-// The output is ordered from largest pole to smallest pole.
-//
-// If the relative difference of two poles is less than tol then they are
-// considered to be multiples.  The default value for tol is 0.001.
-//
-// If the optional parameter reorder is zero, poles are not sorted.
-//
-// The output multp is a vector specifying the multiplicity of the poles.
-// multp(n) refers to the multiplicity of the Nth pole
-// p(idxp(n)).
-//
-//
-//
-
-// test case:
-
-// p = [2 3 1 1 2];
-// [m, n] = mpoles (p)
-// n  =
-//    2.
-//    5.
-//    1.
-//    4.
-//    3.
-// m  =
-//    1.
-//    1.
-//    2.
-//    1.
-//    2.
-//
-
-
-[nargout,nargin]=argn();
-
-  if (nargin < 1 | nargin > 3)
-    error("wrong number of input arguments");
+  if ~( type(p)== 1)
+    error ("mpoles: P must be a single or double floating point vector");
   end
 
-   if (nargin < 2 | isempty (tol))
-     tol = 0.001;
-   end
+  if (nargin < 2 || isempty (tol))
+    tol = 0.001;
+  elseif (~(isscalar (tol) && isreal (tol) && tol > 0))
+    error ("mpoles: TOL must be a real scalar greater than 0");
+  end
 
-   if (nargin < 3 | isempty (reorder))
-     reorder = %t;
-   end
+  if (nargin < 3 || isempty (reorder))
+    reorder = %t;
+  elseif (~(isscalar (reorder) && isreal (reorder)))
+    error ("mpoles: REORDER must be a numeric or logical scalar");
+  end
 
   Np = length (p);
-
-  // Force the poles to be a column vector.
-
-  p = p(:);
-
-  // Sort the poles according to their magnitidues, largest first.
+  p = p(:);  // force poles to be a column vector
 
   if (reorder)
-    // Sort with smallest magnitude first.
-    [p, ordr] = gsort (p,"r","i");
-    // Reverse order, largest maginitude first.
-    n = Np:-1:1;
-    p = p(n);
-    ordr = ordr(n);
+    //// sort with largest magnitude first
+    [_, order] = gsort (abs(p));
+    p = p(order);
   else
-    ordr = 1:Np;
+    order = (1:Np).';
   end
 
-  // Find pole multiplicty by comparing the relative differnce in the
-  // poles.
+  //// Create vector of tolerances for use in algorithm.
+  vtol = zeros (Np, 1, typeof(p));
+  p_nz = (p ~= 0);     // non-zero poles
+  vtol(~p_nz) = tol;  // use absolute tolerance for zero poles
 
-  multp = zeros (Np, 1);
-  indx = [];
+  //// Find pole multiplicity by comparing relative difference of poles.
+  multp = zeros (Np, 1, typeof(p));
+  idxp = [];
   n = find (multp == 0, 1);
-
-
-
-
-
-
   while (n)
-    dp = abs (p-p(n));
-    if (p(n) == 0.0)
-      if (or (abs (p) > 0 & (abs(p)<%inf)))
-        p0 = mean (abs (p(abs (p) > 0 &  abs(p)<%inf)));
-      else
-        p0 = 1;
-      end
-    else
-      p0 = abs (p(n));
+    dp = abs (p - p(n));
+    vtol(p_nz) = tol * abs (p(n));
+    k = find (dp < vtol);
+    //// Poles can only be members of one multiplicity group.
+    if (length (idxp))
+      k = k(~ismember (k, idxp));
     end
-    k = find (dp < tol * p0)';
-    // Poles can only be members of one multiplicity group.
-//    if (length(indx))
-//        mk=members(k,indx);
-//      k = (~ bool2s(mk~=0));
-//    end
     m = 1:length (k);
-    multp(k) = m';
-    indx = [indx; k];
+    multp(k) = m;
+    // disp("k")
+    // disp(k)
+    // disp("idxp")
+    // disp(idxp)
+    idxp = [idxp; k(:)];
     n = find (multp == 0, 1);
   end
-  multp = multp(indx);
-  indx = ordr(indx);
+  multp = multp(idxp);
+  idxp = order(idxp);
+endfunction
 
+function  y = ismember(a, b) // FIXME : do this in a more efficient .
+  y = zeros(size(a,1),size(a,2));
+  for i = 1:length(a)
+    for j = 1:length(b)
+      if a(i) == b(j)
+        y(i) = 1;
+        break;
+      end
+    end
+  end
 endfunction
+
+
+/*
+
+ // test
+ [mp, ip] = mpoles ([0 0], 0.01);
+ assert_checkequal (mp, [1; 2]);
+
+ // test
+ [mp, ip] = mpoles ([-1e4, -0.1, 0]);
+ assert_checkequal (mp, [1; 1; 1]);
+ assert_checkequal (ip, [1; 2; 3]);
+
+// Test single inputs
+// test
+ [mp, ip] = mpoles ([-1e4, -0.1, 0]);
+ assert_checkequal (mp,[1; 1; 1]);
+ assert_checkequal (ip, [1; 2; 3]);
+
+// Test relative tolerance criteria
+// test
+ [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1);
+ assert_checkequal (mp, [1; 1; 1]);
+ [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1 + %eps);
+ assert_checkequal (mp, [1; 1; 2]);
+
+// Test absolute tolerance criteria with a zero pole
+// test
+ [mp, ip] = mpoles ([0, -0.1, 0.3], .1);
+ assert_checkequal (mp, [1; 1; 1]);
+ [mp, ip] = mpoles ([0, -0.1, 0.3], .1 + %eps);
+ assert_checkequal (mp, [1; 1; 2]);
+
+//// Test input validation
+error <Invalid call> mpoles ()
+error <P must be a single or double floating point vector> mpoles (uint8 (1))
+error <TOL must be a real scalar greater than 0> mpoles (1, [1, 2])
+error <TOL must be a real scalar greater than 0> mpoles (1, 1i)
+error <TOL must be a real scalar greater than 0> mpoles (1, 0)
+error <REORDER must be a numeric or logical scalar> mpoles (1, 1, [1, 2])
+error <REORDER must be a numeric or logical scalar> mpoles (1, 1, {1})
+
+*/
\ No newline at end of file
diff --git a/macros/ncauer.sci b/macros/ncauer.sci
index d9a0684..dad3025 100644
--- a/macros/ncauer.sci
+++ b/macros/ncauer.sci
@@ -1,53 +1,85 @@
 // 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
-// Author:Sonu Sharma, RGIT Mumbai
+// Original Source : https://octave.sourceforge.io/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
 // Organization: FOSSEE, IIT Bombay
 // Email: toolbox@scilab.in
-
 function [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
-//Analog prototype for Cauer filter (Cauer filter and elliptic filters are same).
-
-//Calling Sequence
-//[Zz, Zp, Zg] = ncauer(Rp, Rs, n)
-
-//Parameters
-//n: Filter Order
-//Rp: Peak-to-peak passband ripple in dB
-//Rs: Stopband attenuation in dB
-
-//Description
-//It gives an analog prototype for Cauer filter of nth order, with a Peak-to-peak passband ripple of Rp dB and a stopband attenuation of Rs dB.
+    //Analog prototype for Cauer filter (Cauer filter and elliptic filters are same).
+    
+    //Calling Sequence
+    //[Zz, Zp, Zg] = ncauer(Rp, Rs, n)
+    
+    //Parameters
+    //n: Filter Order
+    //Rp: Peak-to-peak passband ripple in dB
+    //Rs: Stopband attenuation in dB
+    
+    //Description
+    //It gives an analog prototype for Cauer filter of nth order, with a Peak-to-peak passband ripple of Rp dB and a stopband attenuation of Rs dB.
+    
+    
+    //Examples
+    //n = 5;
+    //Rp = 5;
+    //Rs = 5;
+    //[Zz, Zp, Zg] = ncauer(Rp, Rs, n)
+    
+    //Zz =
+    //
+    //   0.0000 + 2.5546i   0.0000 + 1.6835i  -0.0000 - 2.5546i  -0.0000 - 1.6835i
+    //
+    //Zp =
+    //
+    //  -0.10199 + 0.64039i  -0.03168 + 0.96777i  -0.10199 - 0.64039i  -0.03168 - 0.96777i  -0.14368 + 0.00000i
+    //
+    //Zg =  0.0030628
+    // Dependencies
+    // ellipap
+    
+    funcprot(0);
+    lhs = argn(1)
+    rhs = argn(2)
+    if (rhs < 3 | rhs > 3)
+    error("ncauer : Wrong number of input arguments.")
+    end
+    
+    [Zz, Zp, Zg] = ellipap(n, Rp, Rs) ;
+    // temp fix to permanently fix this change ellipap
+    Zz = Zz';
+    Zp = Zp';
 
+    endfunction
 
-//Examples
-//n = 5;
-//Rp = 5;
-//Rs = 5;
-//[Zz, Zp, Zg] = ncauer(Rp, Rs, n)
+/* 
 
-//Zz =
-//
-//   0.0000 + 2.5546i   0.0000 + 1.6835i  -0.0000 - 2.5546i  -0.0000 - 1.6835i
-//
-//Zp =
-//
-//  -0.10199 + 0.64039i  -0.03168 + 0.96777i  -0.10199 - 0.64039i  -0.03168 - 0.96777i  -0.14368 + 0.00000i
-//
-//Zg =  0.0030628
+// Test Case 1 (ncauer 0.1, 60, 7)
+[z, p, g] = ncauer(0.1, 60, 7);
+assert_checkalmostequal(z, [2.5574 1.5522 1.3295 -2.5574 -1.5522 -1.3295]*%i, 1e-2);
+assert_checkalmostequal(p, [-0.3664+0.5837*%i -0.1802+0.9088*%i -0.0499+1.0285*%i -0.3664-0.5837*%i -0.1802-0.9088*%i -0.0499-1.0285*%i -0.4796], 1e-2);
+assert_checkalmostequal(g, 7.4425e-03, 1e-2);
 
+// Test Case 2 (ncauer 1.0, 30, 3)
+[z, p, g] = ncauer(1.0, 30, 3);
+assert_checkalmostequal(z, [1.9536 -1.9536]*%i, 1e-2);
+assert_checkalmostequal(p, [-0.2053+0.9870*%i -0.2053-0.9870*%i -0.5597], 1e-2);
+assert_checkalmostequal(g, 0.1490, 1e-2);
 
-funcprot(0);
-lhs = argn(1)
-rhs = argn(2)
-if (rhs < 3 | rhs > 3)
-error("ncauer : Wrong number of input arguments.")
-end
+// Test Case 3 (ncauer 0.25, 50, 6)
+[z, p, g] = ncauer(0.25, 50, 6);
+assert_checkalmostequal(z, [4.0596 1.6414 1.3142 -4.0596 -1.6414 -1.3142]*%i, 1e-2);
+assert_checkalmostequal(p, [-0.4210+0.3665*%i -0.2117+0.8503*%i -0.0550+1.0198*%i -0.4210-0.3665*%i -0.2117-0.8503*%i -0.0550-1.0198*%i], 1e-2);
+assert_checkalmostequal(g, 3.1618e-03, 1e-2);
 
-[Zz, Zp, Zg] = ellipap(n, Rp, Rs) ;
+// Test Case 4 (ncauer 0.8, 45, 4)
+[z, p, g] = ncauer(0.8, 45, 4);
+assert_checkalmostequal(z, [4.1768 1.8543 -4.1768 -1.8543]*%i, 1e-2);
+assert_checkalmostequal(p, [-0.3861+0.4640*%i -0.1234+1.0000*%i -0.3861-0.4640*%i -0.1234-1.0000*%i], 1e-2);
+assert_checkalmostequal(g, 5.6237e-03, 1e-2);
 
-endfunction
+*/
\ No newline at end of file
diff --git a/macros/residue.sci b/macros/residue.sci
index 70dc91f..58303d8 100644
--- a/macros/residue.sci
+++ b/macros/residue.sci
@@ -1,286 +1,318 @@
 // 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
-// Author:[insert name]
+// Original Source : https://octave.sourceforge.io/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
 // Organization: FOSSEE, IIT Bombay
 // Email: toolbox@scilab.in
-
-
+// Dependencies
+// prepad deconv mpoles
 function [r, p, k, e] = residue (b, a, varargin)
 
-
-//  [r, p, k, e] = residue (b, a)
-//  [b, a] = residue (r, p, k)
-// [b, a] = residue (r, p, k, e)
-// The first calling form computes the partial fraction expansion for the
-// quotient of the polynomials, b and a.
-//
-// The quotient is defined as
-
-// B(s)    M       r(m)        N
-// ---- = SUM ------------- + SUM k(i)*s^(N-i)
-// A(s)   m=1 (s-p(m))^e(m)   i=1
-
-// where M is the number of poles (the length of the r, p,
-// and e), the k vector is a polynomial of order N-1
-// representing the direct contribution, and the e vector specifies the
-// multiplicity of the m-th residue's pole.
-//
-//NOTE that the polynomials 'b' and 'a' should have real coefficients(because of the function 'filter' used in polyval)
-//
-//Test case
-//1.
-// b = [1, 1, 1];
-// a = [1, -5, 8, -4];
-// [r, p, k, e] = residue (b, a)
-//    result r = [-2; 7; 3]
-//    result p = [2; 2; 1]
-//    result k = [](0x0)
-//    result e = [1; 2; 1]
-//
-
-//2.
-//[r,p,k,e]=residue([1 2 1],[1 -5 8 -4])
-//OUTPUT
-//r =
-//  -3.0000
-//   9.0000
-//   4.0000
-//
-//p =
-//   2.0000
-//   2.0000
-//   1.0000
-//
-//f = [](0x0)
-//e =
-//   1
-//   2
-//   1
-//
-
-
- [nargout,nargin]=argn();
-
-  if (nargin < 2 | nargin > 4)
-    error ("wrong umber of input arguments");
-  end
-
-  toler = .001;
-
-  if (nargin >= 3)
-    if (nargin >= 4)
-      e = varargin(2);
-    else
-      e = [];
+    if (nargin < 2 || nargin > 4)
+      error("residue: Invalid number of arguments");
     end
-    // The inputs are the residue, pole, and direct part. Solve for the
-    // corresponding numerator and denominator polynomials
-    [r, p] = rresidue (b, a, varargin(1), toler, e);
-    return;
-  end
-
-  // Make sure both polynomials are in reduced form.
-
-  a = polyreduce (a);
-  b = polyreduce (b);
-
-  b = b / a(1);
-  a = a / a(1);
-
-  la = length (a);
-  lb = length (b);
-
-  // Handle special cases here.
-
-  if (la == 0 | lb == 0)
-    k =[];
-    r = [];
-    p = [];
-    e = [];
-    return;
-  elseif (la == 1)
-    k = b / a;
-    r = [];
-    p = [];
-    e = [];
-    return;
-  end
-
-  // Find the poles.
-
-  p = roots (a);
-  lp = length (p);
-
-  // Sort poles so that multiplicity loop will work.
-
-  [e, indx] = mpoles (p, toler, 0);
-  p = p(indx);
-
-  // For each group of pole multiplicity, set the value of each
-  // pole to the average of the group. This reduces the error in
-  // the resulting poles.
-
-  p_group = cumsum (e == 1);
-  for ng = 1:p_group($)
-    m = find (p_group == ng);
-    p(m) = mean (p(m));
-  end
-
-  // Find the direct term if there is one.
-
-  if (lb >= la)
-    // Also return the reduced numerator.
-    [k, b] = deconv (b, a);
+  
+    tol = .001;
+  
+    if (nargin >= 3)
+      if (nargin >= 4)
+        e = varargin(2);
+      else
+        e = [];
+      end
+      // The inputs are the residue2, pole, and direct part.
+      // Solve for the corresponding numerator and denominator polynomials.
+      [r, p] = rresidue (b, a, varargin(1), tol, e);
+      return;
+    end
+  
+    // Make sure both polynomials are in reduced form, and scaled.
+    a = polyreduce (a);
+    b = polyreduce (b);
+  
+    b = b / a(1);
+    a = a/ a(1);
+  
+    la = length (a);
     lb = length (b);
-  else
-    k = [];
-  end
-
-  // Determine if the poles are (effectively) zero.
-
-  small = max (abs (p));
-  if (type(a)==1 | type(b)==1)
-    small = max ([small, 1]) * 1.1921e-07 * 1e4 * (1 + length (p))^2;
-  else
-    small = max ([small, 1]) * %eps * 1e4 * (1 + length (p))^2;
-  end
-  p(abs (p) < small) = 0;
-
-  // Determine if the poles are (effectively) real, or imaginary.
-
-  index = (abs (imag (p)) < small);
-  p(index) = real (p(index));
-  index = (abs (real (p)) < small);
-  p(index) = 1*%i * imag (p(index));
-
-  // The remainder determines the residues.  The case of one pole
-  // is trivial.
-
-  if (lp == 1)
-    r = polyval (b, p);
-    return;
-  end
-
-  // Determine the order of the denominator and remaining numerator.
-  // With the direct term removed the potential order of the numerator
-  // is one less than the order of the denominator.
-
-  aorder = length (a) - 1;
-  border = aorder - 1;
-
-  // Construct a system of equations relating the individual
-  // contributions from each residue to the complete numerator.
-
-  A = zeros (border+1, border+1);
-  B = prepad (matrix (b, [length(b), 1]), border+1, 0,2);
-  for ip = 1:length (p)
-    ri = zeros (size (p,1),size(p,2));
-    ri(ip) = 1;
-    A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0,2).';
-  end
-
-  // Solve for the residues.
-if(size(A,1)~=size(B,1))
-    if(size(A,1)<size(B,1))
-        A=[A;zeros((size(B,1)-size(A,1)),(size(A,2)))];
+  
+    // Handle special cases here.
+    if (la == 0 || lb == 0)
+      k = [];
+      r = [];
+      p = [];
+      e = [];
+      return;
+    elseif (la == 1)
+      k = b / a;
+      r = []
+      p = []
+      e = [];
+      return;
+    end
+  
+    // Find the poles.
+    p = roots (a);
+    lp = length (p);
+  
+    // Sort poles so that multiplicity loop will work.
+    [e, idx] = mpoles (p, tol, 1);
+    p = p(idx);
+  
+    // For each group of pole multiplicity, set the value of each
+    // pole to the average of the group.  This reduces the error in
+    // the resulting poles.
+    p_group = cumsum (e == 1);
+    for ng = 1:p_group($)
+      m = find (p_group == ng);
+      p(m) = mean (p(m));
+    end
+  
+    // Find the direct term if there is one.
+    if (lb >= la)
+      // Also return the reduced numerator.
+      [k, b] = deconv (b, a);
+      lb = length (b);
     else
-        B=[zeros((size(A,1)-size(B,1)),(size(B,2)));B];
+      k = [];
     end
-
+  
+    // Determine if the poles are (effectively) zero.
+    small = max (abs (p));
+    if ( (type(a)==1) || (type(b)==1))
+      small = max ([small, 1]) *  %eps * 1e4 * (1 + length (p))^2;
+    else
+      small = max ([small, 1]) * %eps * 1e4 * (1 + length (p))^2;
     end
-  r = A \ B;
-  r=r(:,$);
-
-endfunction
-
-function [pnum, pden, e] = rresidue (rm, p, k, toler, e)
-
-  // Reconstitute the numerator and denominator polynomials from the
-  // residues, poles, and direct term.
-[nargout,nargin]=argn();
-  if (nargin < 2 | nargin > 5)
-    error ("wrong number of input arguments");
-  end
-
-  if (nargin < 5)
-    e = [];
-  end
-
-  if (nargin < 4)
-    toler = [];
-  end
-
-  if (nargin < 3)
-    k = [];
-  end
-
-  if (length (e))
-    indx = 1:length (p);
-  else
-    [e, indx] = mpoles (p, toler, 0);
-    p = p(indx);
-    rm = rm(indx);
-  end
-
-  indx = 1:length (p);
-
-  for n = indx
-    pn = [1, -p(n)];
-    if (n == 1)
-      pden = pn;
+    p(abs (p) < small) = 0;
+  
+    // Determine if the poles are (effectively) real, or imaginary.
+    idx = (abs (imag (p)) < small);
+    p(idx) = real (p(idx));
+    idx = (abs (real (p)) < small);
+    p(idx) = 1*%i * imag (p(idx));
+  
+    // The remainder determines the residue2s.  The case of one pole is trivial.
+    if (lp == 1)
+      r = polyval (b, p);
+      return;
+    end
+  
+    // Determine the order of the denominator and remaining numerator.
+    // With the direct term removed, the potential order of the numerator
+    // is one less than the order of the denominator.
+    aorder = length (a) - 1;
+    border = aorder - 1;
+  
+    // Construct a system of equations relating the individual
+    // contributions from each residue2 to the complete numerator.
+    A = zeros (border+1, border+1);
+  
+    B = prepad (matrix (b, [length(b), 1]), border+1, 0);
+    B = B(:); // incase b have only 1 element
+    for ip = 1:length (p)
+      ri = zeros (size (p,1),size(p,2));
+      ri(ip) = 1;
+     // A(:,ip) = prepad (rresidue (ri, p, [], tol), border+1, 0).';// invalid index error
+     temprr=rresidue (ri, p, [], tol)
+     ppr=prepad(temprr,border+1,0)
+     A(:,ip) = ppr 
+    end
+  
+    // Solve for the residue2s.
+    // FIXME: Use a pre-conditioner d to make A \ B work better (bug #53869).
+    //        It would be better to construct A and B so they are not close to
+    //        singular in the first place.
+    d = max (abs (A),'c');
+
+    r = (diag (d) \ A) \ (B ./ d);
+  
+  endfunction
+  
+  // Reconstitute the numerator and denominator polynomials
+  // from the residue2s, poles, and direct term.
+  function [pnum, pden, e] = rresidue (r, p, k, tol, e )
+    
+    if nargin < 3 then k = [] ; end
+    if nargin < 4 then tol = [] ; end
+    if nargin < 5 then e = [] ; end
+    if (~isempty (e))
+      idx = 1:length (p);
     else
-      pden = conv (pden, pn);
+      [e, idx] = mpoles (p, tol, 0);
+      p = p(idx);
+      r = r(idx);
     end
-  end
-
-  // D is the order of the denominator
-  // K is the order of the direct polynomial
-  // N is the order of the resulting numerator
-  // pnum(1:(N+1)) is the numerator's polynomial
-  // pden(1:(D+1)) is the denominator's polynomial
-  // pm is the multible pole for the nth residue
-  // pn is the numerator contribution for the nth residue
-
-  D = length (pden) - 1;
-  K = length (k) - 1;
-  N = K + D;
-  pnum = zeros (1, N+1);
-  for n = indx(abs (rm) > 0)
-    p1 = [1, -p(n)];
-    for m = 1:e(n)
-      if (m == 1)
-        pm = p1;
+  
+    idx = 1:length (p);
+    for n = idx
+      pn = [1, -p(n)];
+      if (n == 1)
+        pden = pn;
       else
-        pm = conv (pm, p1);
+        pden = conv (pden, pn);
       end
     end
-    pn = deconv (real(pden),real(pm));
-    pn = rm(n) * pn;
-    pnum = pnum + prepad (real(pn), N+1, 0, 2);
-  end
-
-  // Add the direct term.
-
-  if (length (k))
-    pnum = pnum + conv (pden, k);
-  end
-
-  // Check for leading zeros and trim the polynomial coefficients.
-  if (type(rm)==1 | type(p)==1 | type(k)==1)
-    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * 1.1921e-07;
-  else
-    small = max ([max(abs(pden)), max(abs(pnum)), 1]) *%eps;
-  end
-
-  pnum(abs (pnum) < small) = 0;
-  pden(abs (pden) < small) = 0;
-
-  pnum = polyreduce (pnum);
-  pden = polyreduce (pden);
-
+  
+    // D is the order of the denominator
+    // K is the order of the direct polynomial
+    // N is the order of the resulting numerator
+    // pnum(1:(N+1)) is the numerator's polynomial
+    // pden(1:(D+1)) is the denominator's polynomial
+    // pm is the multiple pole for the nth residue2
+    // pn is the numerator contribution for the nth residue2
+  
+    D = length (pden) - 1;
+    K = length (k) - 1;
+    N = K + D;
+    pnum = zeros (1, N+1);
+    for n = idx(abs (r) > 0)
+      p1 = [1, -p(n)];
+      pn = 1;
+      for j = 1:n - 1
+        pn = conv (pn, [1, -p(j)]);
+      end
+      for j = n + 1:length (p)
+        pn = conv (pn, [1, -p(j)]);
+      end
+      for j = 1:e(n) - 1
+        pn = deconv (pn, p1);
+      end
+      pn = r(n) * pn;
+      pnum = pnum + prepad (pn, N+1, 0, 2);
+    end
+  
+    // Add the direct term.
+    if (length (k))
+      pnum = pnum + conv (pden, k);
+    end
+  
+    pnum = polyreduce (pnum);
+    pden = polyreduce (pden);
+  
+  endfunction
+  
+function y = polyreduce(p)
+    // Remove leading zeros from the polynomial
+    idx = find(p, 1)
+    if isempty(idx) then
+        y = 0;
+    else
+        y = p(idx(1):$);
+    end
 endfunction
+/*
+//test passed
+ b = [1, 1, 1];
+ a = [1, -5, 8, -4];
+ [r, p, k, e] = residue (b, a);
+ assert_checkalmostequal (r, [-2; 7; 3], 1e-12);
+ assert_checkalmostequal (p, [2; 2; 1], 1e-12);
+ assert_checkalmostequal (k,[]);
+ assert_checkalmostequal (e, [1; 2; 1]);
+ k = [1 0];
+ b = conv (k, a) + prepad (b, length (k) + length (a) - 1, 0);
+ a = a;
+ [br, ar] = residue (r, p, k);
+ assert_checkalmostequal (br, b,1e-12);
+ assert_checkalmostequal (ar, a,1e-12);
+ [br, ar] = residue (r, p, k, e);
+ assert_checkalmostequal (br, b,1e-12);
+ assert_checkalmostequal (ar, a,1e-12);
+
+//test passed
+ b = [1, 0, 1];
+ a = [1, 0, 18, 0, 81];
+ [r, p, k, e] = residue (b, a); // error filter: Wrong type for input argument #1: Real matrix or polynomial expect
+ FIXME :  builtin filter doesn't support complex paramters 
+ r1 = [-5*%i; 12; +5*%i; 12]/54;
+ p1 = [+3*%i; +3*%i; -3*%i; -3*%i];
+ assert_checkalmostequal (r, r1, 1e-12);
+ assert_checkalmostequal (p, p1, 1e-12);
+ assert_checkalmostequal (k,[]);
+ assert_checkalmostequal (e, [1; 2; 1; 2]);
+ [br, ar] = residue (r, p, k);
+ assert_checkalmostequal (br, b, 1e-12);
+ assert_checkalmostequal (ar, a, 1e-12);
+
+//test passed
+ r = [7; 3; -2];
+ p = [2; 1; 2];
+ k = [1 0];
+ e = [2; 1; 1];
+ [b, a] = residue (r, p, k, e);
+ assert_checkalmostequal (b, [1, -5, 9, -3, 1], 1e-12);
+ assert_checkalmostequal (a, [1, -5, 8, -4], 1e-12);
+
+ 
+ [rr, pr, kr, er] = residue (b, a);
+ [_, m] = mpoles (rr);
+ [_, n] = mpoles (r);
+ assert_checkalmostequal (rr(m), r(n), 1e-12);
+ assert_checkalmostequal (pr(m), p(n), 1e-12);
+ assert_checkalmostequal (kr, k, 1e-12);
+ assert_checkalmostequal (er(m), e(n), 1e-12);
+
+//test passed
+ b = [1];
+ a = [1, 10, 25];
+ [r, p, k, e] = residue (b, a);
+ r1 = [0; 1];
+ p1 = [-5; -5];
+ assert_checkalmostequal (r, r1, 1e-12);
+ assert_checkalmostequal (p, p1, 1e-12);
+ assert_checkalmostequal (k,[]);
+ assert_checkalmostequal (e, [1; 2]);
+ [br, ar] = residue (r, p, k);
+ assert_checkalmostequal (br, b, 1e-12);
+ assert_checkalmostequal (ar, a, 1e-12);
+
+// The following //test is due to Bernard Grung
+//test <*34266> passed
+ z1 =  7.0372976777e6;
+ p1 = -3.1415926536e9;
+ p2 = -4.9964813512e8;
+ r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1;
+ r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1;
+ r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1;
+ r4 = z1/p2/p1;
+ r = [r1; r2; r3; r4];
+ p = [p1; p2; 0; 0];
+ k = [];
+ e = [1; 1; 1; 2];
+ b = [1, z1];
+ a = [1, -(p1 + p2), p1*p2, 0, 0];
+ [br, ar] = residue (r, p, k, e);
+ assert_checkalmostequal (br, [0,0,b],%eps,1e-8);
+ assert_checkalmostequal (ar, a, %eps,1e-8);
+
+//test <*49291> passed
+ rf = [1e3, 2e3, 1e3, 2e3];
+ cf = [316.2e-9, 50e-9, 31.6e-9, 5e-9];
+ [num, den] = residue (1./cf,-1./(rf.*cf),0);
+ assert_checkalmostequal (length (num), 4);
+ assert_checkalmostequal (length (den), 5);
+ assert_checkalmostequal (den(1), 1);
+
+//test <*51148> passed
+ r = [1.0000e+18, 3.5714e+12, 2.2222e+11, 2.1739e+10];
+ pin = [-1.9231e+15, -1.6234e+09, -4.1152e+07, -1.8116e+06];
+ k = 0;
+ [p, q] = residue (r, pin, k);
+ assert_checkalmostequal (p(4), 4.6828e+42, 1e-5);
+
+//test <*60384> passed
+ B = [1315.789473684211];
+ A = [1, 1.100000536842105e+04, 1.703789473684211e+03, 0];
+ poles1 = roots (A);
+ [r, p, k, e] = residue (B, A);
+ [B1, A1] = residue (r, p, k, e);
+ assert_checkalmostequal (B, B1);
+ assert_checkalmostequal (A, A1);
+
+%/