FiltFilt function initial and final conditions
50 views (last 30 days)
Show older comments
Hi,
I would like to filter a signal using filtfilt in order to have:
- zero phase - to have the initial and final sample of the filtered data the same as the raw data
Following this link: http://www.mechanicalvibration.com/filtfilt_Casual_versus_non_.html#foot60 the initial and final point are maintained using filtfilt. However when I do this, they are not.
The commands I use are: b = [0.0675, 0.1349, 0.0675]; % Low pass Butterworth filter a = [1.0, -1.1430, 0.4128]; Sig_filtered=filtfilt(b, a, Sig_raw);
Can anyone explain this to me or point me towards a filter function that can do this?
Thank you in advance.
Kr, Brecht
0 Comments
Accepted Answer
Jan
on 29 Aug 2018
As written in the documentation, and as you can see inside the code, filtfilt mirrors the initial and final values to reduce the transitional effects. This cannot preserve the marginal values in every case.
It is not clear, what you exactly want to achieve. Please explain uniquely and clearly what you need.
6 Comments
Jan
on 16 Oct 2018
Edited: Jan
on 16 Oct 2018
@Robert: In the documentation of filtfilt.m you find the source of the applied method:
Frederik Gustafsson, Determining the initial state in forward-backward filtering, IEEE Transactions on Signal Porcessing, pp 998-992, April 1996, Volume 44, Issue 4
Therefore I disagree with the opinion, that MathWorks is sloppy with the exact description of the mathematical decisions.
From a physical point of view, the filter is a resonator. The filtered output depends on the input signal and the state of the internal "oscillators". Then the method of reducing the transitional effects in filtfilt can be interpreted such, that it is assumed, that the signal can be expanded by subtracting the flipped initial phase from the first element. This is a fair guess for a sine wave, for a static signal, for pure noise. The spectrum of the signal is conserved approximately and this does not create a DC offset. But of course e.g. for a square or triangle wave the results are not perfect. There cannot be a perfect method in general. As soon as you do have knowledge about the real signal in before and after the filtered section, it is wise to use it.
Is this "built in" or is there source to this function?
It seems like you have found out already, that the code is shipped with Matlab. Simply try
type filtfilt.m
You find the code there and the above mentioned reference for the algorithm. This reference is mentioned in the documentation also: doc filtfilt (link)
Robert Bristow-Johnson
on 16 Oct 2018
Edited: Robert Bristow-Johnson
on 16 Oct 2018
Bruno, TMW and Cleve know that i am an historic pain-in-arse regarding MATLAB and will continue to be a grumpy old man as long as the indexing origin is hard-wired to 1 (i still can't fathom that, such a basic and fundamental problem). but now i have at least a couple of audio friends that work for TMW and they gave me a cute little T shirt so now i am less critical than i was on comp.dsp and comp.soft-sys.matlab in the '90s and early 2000s.
but even today, it's disgusting. my MATLAB code is disgusting. ugliest friggin' code i have ever written. it's all ugly ugly ugly. and that 1 origin convention is the main reason. not being able to have negative indices is another reason. and the order of coefficients for polyval() and polyfit() is wrong.
so i guess i am still a grumpy old DSPer.
More Answers (3)
Bruno Luong
on 16 Oct 2018
Edited: Bruno Luong
on 16 Oct 2018
Here is a demo of using filtfilt on periodic data
- The first method use filtfilt on one period
- The second method stitches rawdata one period on the left and one on the right, call filtfilt on the three periods, then crop the middle period (Jan's suggestion)
- The third method does like the second but smoothly interpolate in a transition of 0.1 period.
As you can see the first one have a problem of stitching at the boundary, the second and third methods give almost the same result, but I know the third one must have a smooth transition.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/198062/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/198063/image.png)
% True Periodic signal
x = linspace(0,2*pi,100);
y = sin(x) - 0.5*sin(2*x+1) + 0.6*sin(3*x+2) - 0.3*sin(4*x+4);
% Add noise
sigma = 0.3;
y = y + sigma*randn(size(y));
% Filter design
[b,a] = butter(8,0.1);
% First method
yff = filtfilt(b,a,y);
% Second & Third methods
nperiods = 3;
[xd,yd] = duplicate(x,y,nperiods);
y123 = filtfilt(b,a,yd);
y1 = crop(y123,nperiods,1);
y2 = crop(y123,nperiods,2);
% Second method
yt = y2;
% Third method
w = max(0,interp1(x(end)*[0.9 1],[0,1],x,'linear','extrap'));
yi = (1-w).*y2 + w.*y1;
% Generate graphics
close all
hold on
xp = x/(2*pi);
h(1)=plotdup(xp,y,'color',0.7*[1 1 1]);
h(2)=plotdup(xp,yff,'b');
h(3)=plotdup(xp,yt,'g');
h(4)=plotdup(xp,yi,'r');
xline(1);
xline(2);
legend(h,'noisy data','fiiltfilt','ff-truncated','ff-interpolate')
function [xd,yd] = duplicate(x,y,nperiods)
if nargin < 3
nperiods = 3;
end
x = x(:);
dx = x(end)-x(1);
xd = x(1:end-1) + dx*(0:nperiods-1);
xd = [xd(:)', dx*(nperiods-1)+x(end)];
yd = [repmat(y(1:end-1),1,nperiods), y(end)];
end
function y = crop(yd,nperiods,i)
yend = yd(end);
yd = reshape(yd(1:end-1),[],nperiods);
if i < nperiods
yend = yd(1,i+1);
end
y = [yd(:,i); yend]';
end
function h = plotdup(x,y,varargin)
[xd,yd] = duplicate(x,y);
h = plot(xd,yd,varargin{:});
end
0 Comments
Bruno Luong
on 16 Oct 2018
Edited: Bruno Luong
on 16 Oct 2018
Now a rigorous approach of zero-phase FIR on circular data. One need to build a circular matrix from the coefficients A and B, the similar to filt-filt apply the circular-filter in both directions to undo the phase-shift.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/198081/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/198083/image.png)
I call it fourth method, a ff-circular. The code is following
% True Periodic signal
x = linspace(0,2*pi,100);
y = sin(x) - 0.5*sin(2*x+1) + 0.6*sin(3*x+2) - 0.3*sin(4*x+4);
% Add noise
sigma = 0.3;
y = y + sigma*randn(size(y));
% Filter design
[b,a] = butter(8,0.1);
% First method
yff = filtfilt(b,a,y);
% Fourth method
n = length(y)-1;
Sa = cmat(a,n);
Sb = cmat(b,n);
yc=Sa\(Sb*y(1:end-1)');
yc=Sa\(Sb*flip(yc));
yc = flip(yc([end 1:end]))';
% Generate graphics
close all
hold on
xp = x/(2*pi);
h(1)=plotdup(xp,y,'color',0.7*[1 1 1]);
h(2)=plotdup(xp,yff,'b');
h(3)=plotdup(xp,yc,'r');
xline(1);
xline(2);
legend(h,'noisy data','fiiltfilt','ff-circular')
function S = cmat(a, n)
A = repmat(flip(a),n);
d = 0:length(a)-1;
S = spdiags(A,d,n,n) + ...
spdiags(A,d-n,n,n);
end
function [xd,yd] = duplicate(x,y,nperiods)
if nargin < 3
nperiods = 3;
end
x = x(:);
dx = x(end)-x(1);
xd = x(1:end-1) + dx*(0:nperiods-1);
xd = [xd(:)', dx*(nperiods-1)+x(end)];
yd = [repmat(y(1:end-1),1,nperiods), y(end)];
end
function h = plotdup(x,y,varargin)
[xd,yd] = duplicate(x,y);
h = plot(xd,yd,varargin{:});
end
0 Comments
Brecht Vermeulen
on 4 Sep 2018
4 Comments
Bruno Luong
on 16 Oct 2018
Edited: Bruno Luong
on 16 Oct 2018
@Jan: A better way than crude cropping is interpolating the filtered signal over overlapping zone (imagine the signal wrapped on the circle) by a weighted of two signals that make a smooth transitions with continuity, the weight selected to can be infinity smooth if such smoothness is desired.
Jan
on 16 Oct 2018
Edited: Jan
on 16 Oct 2018
@Bruno: I agree. I'm only not satisfied with such photoshopping of measurement data. Smoothing a signal until it fits the expectations exactly is a construction of a result and not necessarily a measurement. Another problem could be a high frequency signal near to the Nyquist frequency. Then a minimal timeshift can cause a cancellation of the signal. So cropping is less smart, but perhaps this crudeness preserves the nature of the original signal.
As usual in signal processing, filtering is not trivial. You cannot expect that filter or filtfilt remove the noise magically and a deep knowledge of the effects of filtering is required. A pure noise can be filtered, until a clean periodic sine-wave is found as "result", but this has no scientific power anymore.
But maybe the OP is aware of such problems. It was only "heavy filtering" combined with "maintain exactly the initial and final sample" which let me mention the potential overdoing at filtering.
See Also
Categories
Find more on Digital Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!