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."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." 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." filtfilt, zero-phase filter MATLAB Answers — New Questions
