1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
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)
*/
|
