filtfilt Does Not Perform Zero-Phase Filtering. Does It Matter?

The doc page for filtfilt claims that the function "performs zero-phase digital filtering" and then goes on to describe the expected characteristics of a zero-phase output. It's pretty easy to show that this claim is false.

Generate some input data

rng(100);
x = rand(1,50);

Define a simple FIR filter

b = 1:4;
a = 1;

The zero-phase filtered output y and its associated discrete times can be computed as

h = b;         % filter impulse resonse
y2 = conv(fliplr(h),conv(h,x)); % zero-phase, filtered output

Nh = numel(h); % length of impulse response
Nx = numel(x); % length of input
ny2 = -(Nh-1) : Nh + Nx - 2;  % time of output

As must be the case, the output is longer than the input, and the output extends to the left into negative time.

We can show that y2 satisfies the properties for zero-phase filtering

Compute the DTFT of the input and the impulse response

[X,w]  = freqz(x,1,1000);
H      = freqz(b,a,w);

Compute the DTFT of the output

Y2     = freqz(y2,1,w).*exp(-1j*ny2(1)*w);

The phase of Y2 relative to X is (essentially) zero

figure
plot(w,angle(Y2./X)),xlabel('rad'),ylabel('rad')

The magnitude of Y2 relative to X is the square of the magitude of the transfer function

figure
semilogy(w,abs(Y2./X),'DisplayName','abs(Y2/X)')
hold on
semilogy(w,abs(H).^2,'ro','MarkerIndices',1:20:numel(w),'DisplayName','abs(H)^2')
xlabel('rad'),ylabel('Amplitude')
legend

Now compare the output of filtfilt() to the zero-phase output

y4 = filtfilt(b,a,x);
ny4 = 0:numel(y4) - 1;

figure
stem(ny2,y2,'DisplayName','zero-phase')
hold on
stem(ny4,y4,'x','DisplayName','filtfilt')
legend("Location",'Best')

The differences are apparent at the edges, undoubtedly because filtfilt() also takes actation that "minimizes start-up and ending transients." However, such action kills the "zero-phase" intent of the function. Even though it looks like the output of filtfilt() is close to the zero-phase output, the phase response really isn't even close to zero-phase and the magnitude response isn't abs(H)^2, even after taking into account that filtfilt() only returns the values of the output that correspond to the time samples of the input (it doesn't return the computed values that extend past the left and right edges of the input).

Edit 13 April 2023

Here is the frequency domain assesment of the output of filtfilt as seen by the user. I didn't include it in the original post because the output of filtfilt to the user doesn't include the extra points that extend off the edges of the result that filtfilt computes internally but chops off before returning to the user to make the output points correspond to the input points. However, I've looked at those when debugging and they don't really change the story and I thinkg these plots may be of interest.

Y4     = freqz(y4,1,w);
figure
semilogy(w,abs(Y4./X),'DisplayName','abs(Y4/X)')
hold on
semilogy(w,abs(H).^2,'ro','MarkerIndices',1:20:numel(w),'DisplayName','abs(H)^2')
xlabel('rad'),ylabel('Amplitude')
legend

figure
plot(w,angle(Y4./X)),xlabel('rad'),ylabel('rad')

Are there applications where the true, zero-phase response (or at least its portion that corresponds to the time samples of the input) is desired?

Should filtfilt() include an option to not minimize transients for applications that need a true, zero-phase filter?

At a minimum, shouldn't the filtfilt() doc page be updated?

Related to that last question, I think it's interesting (if not telling) that other functions, like lowpass, that I think operate nearly exactly, if not exactly, the same as filtfilt() only claim that the function "compensates for the delay introduced by the filter."

[Paul]

Here is an example of an IIR filter taken from this thread that shows how to use filter (instead of conv) to achieve an approximate zero-phase response (a true zero-phase response for an IIR filter is not numerically computable, in general, because it would be an infinite length response).

Design the filter of interest

bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
  'PassbandFrequency2', 3600, 'StopbandFrequency2', 3700, 'StopbandAttenuation1', 30,...
  'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');
% apply filtfilt to a random (white) long input signal; output signal shows an undesirable transient
rng(100);
x = randn(1,2^20);

Don't run filtfilt here, see linked thread for further discussion on this example

%y = filtfilt(bp,x);

Apply forward/backward filtering w/o padding.

% forward/backward filter w/o padding
y2 = flip(filter(bp,flip(filter(bp,x))));

Apply forward/backward filtering with padding going in both directions. For this example, a pad length of 2x the length of the impulse response seems to work well. A generalized approach would be needed for a general solution

Npad = 2*numel(impz(bp));
y3full = flip(filter(bp,[flip(filter(bp,[x,zeros(1,Npad)])),zeros(1,Npad)]));

The corresponding input sequence and discrete-time vector

x3full = [zeros(1,Npad),x,zeros(1,Npad)];
nfull  = -Npad:numel(x)-1+Npad;

Compare the magnitude of the input/output response to the square of the magnitude of the filter frequency response and plot the phase of the input/output response, which should be zero (not expected to be perfect, particularly where the magnitude is small)

[h3,w] = tfestimate(x3full,y3full);
hbp = freqz(bp,w);
figure
plot(subplot(211),w,db(abs(h3)),w,db(abs(hbp).^2)),grid,ylabel('Magnitude (db)')
plot(subplot(212),w,angle(h3)*180/pi),grid
ylim([-1,1]);ylabel('Phase (deg)'),xlabel('\omega (rad/sample)')

Extract the portion of the zero-phase response that corresponds to the original input signal

y3     = y3full((1:numel(x))+Npad);

Compare to the forward/backward output w/o padding, we see a difference at one end (I'm actually surprised there's not a difference at both ends).

e = abs(y2-y3);
figure
plot(e)

figure
plot(e)
xlim([numel(e)-500,numel(e)])

In this example, that difference over the small interval at the end doesn't really make any visible difference in the estimated frequency response based on x and y2 (not shown), probably because the input signal is so long.

What happens with a shorter input?

x = x(1:2^10);
y2 = flip(filter(bp,flip(filter(bp,x))));
y3full = flip(filter(bp,[flip(filter(bp,[x,zeros(1,Npad)])),zeros(1,Npad)]));
x3full = [zeros(1,Npad),x,zeros(1,Npad)];
[h3,w] = tfestimate(x3full,y3full);
hbp = freqz(bp,w);
figure
plot(subplot(211),w,db(abs(h3)),w,db(abs(hbp).^2)),grid,ylabel('Magnitude (db)')
plot(subplot(212),w,angle(h3)*180/pi),grid
ylim([-1,1]);ylabel('Phase (deg)'),xlabel('\omega (rad/sample)')

y3 = y3full((1:numel(x))+Npad);
e = abs(y2-y3);
figure
plot(e)

h2 = tfestimate(x,y2,w);
figure
plot(subplot(211),w,db(abs(h2)),w,db(abs(hbp).^2)),grid,ylabel('Magnitude (db)')
plot(subplot(212),w,angle(h2)*180/pi),grid

Now we see that not padding has a clear effect on the estimated frequency, though to be fair a similar result would be obtained if only considering the elements of y3full that correspond to the original input, x.

Answer 1

[Bruno Luong]

In the theory of linear filter, there is no such thing as start/stop, the signal extend to infinity in both ends.

The definition of "zero-phase" is that a filter is zero-phase if and only if the transfer function in z-domain is purely real.

filtfilt transfer function satisfies this requirement so it is a zero-phase filter.

[Paul]

I agree with the first statement, though not sure how it's relevant to this Question. It's also a basic point about signals that I feel tends to be get overlooked in this forum.

If by "transfer function in the z-domain" you mean its frequency response, I'd add that in addition to the frequency response being real, it must also be positive and an even function of frequency.

I'm not sure how to interpret the third statement in light of the fact that filtfilt() does more than just forward/reverse filter the input sequence.

filtfilt claims that

"The result has these characteristics:

• Zero phase distortion.
• A filter transfer function equal to the squared magnitude of the original filter transfer function"

The first claim is not true, as I've implicitly shown above. The second claim is, at best, misleading because the algorithm does more than just apply the input transfer function (twice) to the input data. That is, it reads as if the DT Fourier transform of the output, i.e., "the result," satisfies Y(w) = abs(H(w))^2*U(w), which isn't true. It's also worded poorly insofar as the result doesn't "have" a transfer function.