1 files changed, 104 insertions, 0 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)
+*/