Corrected the code

3 files changed, 127 insertions, 57 deletions

diff --git a/macros/filtfilt.sce b/macros/filtfilt.sce
new file mode 100644
index 0000000..c8b9aa2
--- /dev/null
+++ b/macros/filtfilt.sce
@@ -0,0 +1,104 @@
+//This is a Signal Processing toolbox function

+//Author: Rashmi Patankar, FOSSEE IIT Bombay

+// y = filtfilt (b, a, x)

+//Forward and reverse filter the signal. 

+//This corrects for phase distortion introduced by a one-pass filter, though it does square the magnitude response in the process. That’s the theory at least. In practice the phase correction is not perfect, and magnitude response is distorted, particularly in the stop band. 

+//Example 

+//[b, a]=butter(3, 0.1);                  # 5 Hz low-pass filter

+//t = 0:0.01:1.0;                         # 1 second sample

+//x=sin(2*pi*t*2.3)+0.25*randn(size(t));  # 2.3 Hz sinusoid+noise

+//y = filtfilt(b,a,x); z = filter(b,a,x); # apply filter

+//plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')

+ 

+

+function y = filtfilt(b, a, x)

+

+  // Check for correct number of input arguments

+  if nargin() ~= 3 then

+    error("filtfilt: Wrong number of input arguments");

+  end

+

+  // Handle row vectors

+  rotate = (size(x,1)==1);

+  if rotate then

+    x = x(:); // Make it a column vector

+  end

+

+  // Get dimensions and lengths

+  lx = size(x,1);

+  a = a(:).';

+  b = b(:).';

+  lb = length(b);

+  la = length(a);

+  n = max(lb, la);

+  lrefl = 3 * (n - 1);

+  if la < n then 

+      a(n) = 0; 

+  end

+  if lb < n then 

+      b(n) = 0; 

+  end

+

+  // Check input length

+  if (size(x, 1) <= lrefl) then

+    error("filtfilt: X must be a vector or matrix with length greater than " + string(lrefl));

+  end

+

+  // Compute initial state

+  kdc = sum(b) / sum(a);

+  if (abs(kdc) < %inf) then // Neither NaN nor +/- Inf

+    si = flipdim(cumsum(flipdim(b - kdc * a, 2)), 2);

+  else

+    si = zeros(size(a)); // Fall back to zero initialization

+  end

+  si(1) = [];

+  si = [si 0];  // Append a zero to si

+  si = si(1:(max(length(a), length(b)) - 1));

+  si = si';

+  disp(size(si));

+  // Filter each column

+  for c = 1:size(x,2)

+   // v = x(:);  // Convert single column to row vector

+    v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c); 2*x($,c)-x(($-1):-1:$-lrefl,c)];

+

+    // Forward and reverse filtering

+    v = filter(b, a, v, si*v(1));  // Forward filter

+    v = flipdim(filter(b, a, flipdim(v, 1), si*v($)), 1);  // Reverse filter

+    y(:,c) = v((lrefl+1):(lx+lrefl));

+    y(:,c) = y(:)';  // Transpose back to column vector

+  end

+

+  // Handle row vectors

+  if rotate then

+    y = flipdim(y, 2) // Rotate it back

+  end

+

+endfunction

+

+/*

+Test case 1:

+b = [1 2];  

+a = [1 -0.5];

+x = [1 2 3 4];

+y = filtfilt(b, a, x)

+Output: 38.396851  71.793701  105.08740  136.92480

+

+Test case 2:

+b = [0.2, 0.2];

+a = [1, -0.8];

+x = [1, 2, 3, 4, 5];

+y = filtfilt(b, a, x)

+Output:  5.0357884  7.2211355  9.3675393  11.382320  13.166217

+

+Test case 3:

+b = [1 2 5];

+a = [1 6 8];

+x = [1 2 3 4 5 6 7];

+y = filtfilt(b, a, x)

+

+Test case 4:

+b = [1, 2, 5, 7, 9];

+a = [1, -0.5, 8, 4, 3];

+x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 3, 5, 7];

+y = filtfilt(b, a, x)

+*/

diff --git a/macros/filtfilt.sci b/macros/filtfilt.sci
deleted file mode 100644
index 7f86113..0000000
--- a/macros/filtfilt.sci
+++ /dev/null
@@ -1,23 +0,0 @@
-function [y]=filtfilt(b,a,x)
-
-// Zero phase digital filtering
-// Calling Sequence
-//	[y]=filtfilt(b,a,x)
-// Parameters
-//	b: Real or complex valued vector or matrix
-//	a: Real or complex valued vector or matrix
-//	x: Real or complex valued vector or matrix
-// Description
-//	This is an Octave function
-//	In theory, it forwards and reverse filters the signal and corrects phase distortion upto an extent by a one-pass filter but squares the magnitude response in the process. Practically though, the correction isn't perfect and magnitude response, particularly the stop band is distorted.
-// Examples
-// 1.	y=filtfilt (1,2*%i,[%i -4 0])        // Number of Output argument should be equal to 1
-//	y =   [-0.25i   1   0]
-
-funcprot(0);
-rhs=argn(2);
-if (rhs~=3) then
-	error ("Wrong number of input arguments.")
-else y=callOctave("filtfilt",b,a,x)
-end
-endfunction
diff --git a/macros/gaussian.sci b/macros/gaussian.sci
index 59ddb3c..234512e 100644
--- a/macros/gaussian.sci
+++ b/macros/gaussian.sci
@@ -1,34 +1,23 @@
-function w = gaussian (m, a)
-//This function returns a Gaussian convolution window.
-//Calling Sequence
-//w = gaussian (m)
-//w = gaussian (m, a)
-//Parameters 
-//m: positive integer value
-//a: 
-//w: output variable, vector of real numbers
-//Description
-//This is an Octave function.
-//This function returns a Gaussian convolution window of length m supplied as input, to the output vector w.
-//The second parameter is the width measured in sample rate/number of samples and should be f for time domain and 1/f for frequency domain. The width is inversely proportional to a.
-//Examples
-//gaussian(5,6)
-//ans  =
-//    5.380D-32  
-//    1.523D-08  
-//    1.         
-//    1.523D-08  
-//    5.380D-32  
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<1 | rhs>2)
-error("Wrong number of input arguments.")
-end
-if(rhs==1)
-w = callOctave("gaussian",m)
-elseif(rhs==2)
-w = callOctave("gaussian",m, a)
-end
-endfunction
-
+//Author: Rashmi Patankar

+//This function returns a Gaussian convolution window.

+//Calling Sequence

+//w = gaussian (m)

+//w = gaussian (m, a)

+//Parameters 

+//m: positive integer value

+//a: integer value

+//w: output variable, vector of real numbers

+//Description

+//This function returns a Gaussian convolution window of length m supplied as input, to the output vector w.

+//The second parameter is the width measured in sample rate/number of samples and should be f for time domain and 1/f for frequency domain. The width is inversely proportional to a.

+function w = gaussian(m, a)

+    if nargin < 1 | nargin > 2 then

+        error('gaussian: Incorrect number of input arguments');

+    elseif ~(isscalar(m) & m == fix(m) & m > 0) then

+        error('gaussian: M must be a positive integer');

+    elseif nargin == 1 then

+        a = 1;

+    end

+

+    w = exp(-0.5 * (([0:m-1]' - (m-1)/2) .* a).^2);

+endfunction