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

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.