FiltFilt function initial and final conditions

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

Answer 1

[Jan]

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.

[Robert Bristow-Johnson]

What, exactly, do you mean by "... filtfilt() mirrors the initial and final values to reduce the transitional effects."?

It's probably the same as me asking what, exactly, does the documentation mean by " filtfilt() minimizes start-up and ending transients by matching initial conditions."?

We have an IIR filter. What do they set the initial states to for each forward and backward pass?

I can imagine setting the IIR states to zero, doing the forward filtering and extending the length to let the IIR filter "ring out". Then reversing the extended-length intermediate sequence, zeroing the IIR states, filtering again until the end sample (we don't need to compute and save the ringout samples here). Then truncating the beginning samples (which are the filtered output of the previous "ring out" samples of the forward filtering) and finally reversing again back to the original order. But I do not see how this prevents either end from having a jump discontinuity.

Really, I just don't understand what filtfilt() does regarding these transients or what they do about the initial IIR states. It seems to me that filtfilt() should, to be totally legit, extend the length on both ends.

[Robert Bristow-Johnson]

I want to be clear about exactly what y=filtfilt(b,a,x) does.

is it this:

y = flipud( filter(b, a, flipud(filter(b, a, x)) ) );

?

It seems to me that it is not exactly that regarding the edges. But neither "mirrors the the initial and final values to reduce the transitional effects" nor "minimizes start-up and ending transients by matching initial conditions" suffice to spell out what filtfilt() actually does mathematically at the edges. I hope this isn't TMW proprietary information because we users have to know what are the mathematical consequences of using function calls in MATLAB.

Is this "built in" or is there source to this function?

[Robert Bristow-Johnson]

okay, i am looking at the source and getting a hint, but i have questions about the rationale.

it appears that the length of mirrored segment is 3 times the order of the filter. (it's this nfact number.) why? it seems to me that the number of samples to extend would depend on the Q of the dominant pole(s). like the time-constant is about twice the Q for a resonant frequency of 1. so then you need to go out about 5 time constants w.r.t. the dominant pole. what does that have to do with the filter order?

it appears that the mirroring is even-symmetry mirroring overlapping the edge sample which effectively doubles it. why? why not odd-symmetry mirroring? or why not just extend the value of the edge sample out as a constant for about 3 times the order?

i just do not understand many "why?s" in MATLAB functions. some of the decisions about function and convention seem arbitrary and not canonical.

[Bruno Luong]

+1 with you Robert, TMW is in general sloppy in their decision and then also in the documentation, because actually they might be not proud to write it down mathematically and rigorously.

For example the documentation of Mathematica is entirely another level.

I have opportunity to take a look of some of those filter tools, and many of them has been written 20 years ago. Some of them looks too complicated and seems the intern algorithm is not well documented. I guess it works but they lost their hand on that because the programmer is non longer with them.

[Jan]

@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]

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.

Answer 2

[Bruno Luong]

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.

% 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

Answer 3

[Bruno Luong]

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.

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

Answer 4

[Brecht Vermeulen]

I expected filtfilt to maintain exactly the initial and final sample. When heavily filtering this deviation might become quite big.

What I exactly needed is filtering of a cyclic signal, where the start and endpoint were exactly the same. I manually have written an adjustment code that did this.

Thank you for the clarification.

[Bruno Luong]

If your data is periodic, than not only you need to match the value of the first and end point but also the derivative of those points.

I'm not sure but to me it looks like one can write down the recursive filter equations in a periodic way and solve it.

I guess it still give a linear relation that can non longer apply recursively but rather can be solved by linear algebra and matrix.

Note: For non-recursive filter (not your case), it boils down to circular convolution.

[Jan]

 I expected filtfilt to maintain exactly the initial and
 final sample.

This is not the nature of filtfilt, which expands the original signal by a mirrored signal of the initial and final phase to allow a more or less generally matching smoothing in these parts also.

Transitional effects are a problem at filtering, which can be solved heuristically only. Any algorithm must either apply a poorer filtering at the corners due to less available data, or it must guess a "valid" continuation of the signal, to set the internal state of the filter to the matching parameters. Both methods have drawbacks. But if you know, that the signal is periodic, you have all you need to expand the signal: Add the final phase to the start, the start phase to the end, run filter in both directions or try filtfilt. Afterwards crop the added signals. This does not guarantee, that the first and last sample are perfectly identical, but if you have noise in the signal, this would not be realistic at all.

[Bruno Luong]

@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]

@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.