Butterworth filter vs Elliptic filter - artefacts and huge transient

Question

I'm working on a project involving digital filters and DFT. I want to use sampled ECG data to determine the heart rate using power spectral density. While filtering the ECG samples (samples are from this database), I noticed a significant transient when applying an IIR high-pass Butterworth filter. Here's the code and the output:

clear;
clc;

fc = 360;
T = 1/fc; 

ecg_samples = load('105m (0).mat').val;
N = length(ecg_samples);
times = (0:N-1)*T;
ecg_samples = (ecg_samples-ecg_samples(1))./(max(ecg_samples)-ecg_samples(1));


fstopband_hhp = 0.5;
fpass_hhp = 0.6;   

hhp_butter = designfilt("highpassiir", ...
    PassbandFrequency = fpass_hhp, ...
    StopbandFrequency = fstopband_hhp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc, ...
    DesignMethod = "butter");

hhp_elliptic = designfilt("highpassiir", ...
    PassbandFrequency = fpass_hhp, ...
    StopbandFrequency = fstopband_hhp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc, ...
    DesignMethod = "ellip");

ecg_samples_filtered1 = filtfilt(hhp_butter,ecg_samples);
ecg_samples_filtered2 = filtfilt(hhp_elliptic,ecg_samples);

subplot(3,1,1);
plot(times,ecg_samples);
xlabel("time");
ylabel("amplitude");
title("before filtering");

subplot(3,1,2);
plot(times,ecg_samples_filtered1);
xlabel("time");
ylabel("amplitude");
title("butterworth filter");

subplot(3,1,3);
plot(times,ecg_samples_filtered2);
xlabel("time");
ylabel("amplitude");
title("elliptic filter");

As you can see the transient is not present using an elliptic filter with the same specs. Why does this huge transient appear only when using a Butterworth filter with these specific specs? If I change the filter cutoff frequencies the Butterworth filter doesn't produce the transient, e.g. changing from fpass_hhp = 0.6; to fpass_hhp = 0.7;

Answer 1

Here is what's happening

Your filter specification is unsuitable for a Butterworth filter especially with such a low corner frequency. The dead give away here is the filter order is which is 42 (which filtfilt() then doubles to a whopping 84).

The issue with such a high IIR filter order is numerically stability. Since the filter is recursive any rounding errors get amplified so you need sufficient numerical precision to keep the error low enough. For your filter even 64-bit floating isn't good enough. The elliptic filter works a lot better here since it only got an order of 8.

In general, designing a Butterworth using stopband attenuation and passband ripple is not a good fit. You are much better off just specifying the corner frequency and varying the order until you see something you like. Butterworth filter are flat in the passband have monotonically increasing stopband attenuation, so the concept of "ripple" or "min attenuation" doesn't really apply.

We can give you a few guidelines on filter design, but it is a complicated topic with some really non-trivial math underneath it, so I always suggest taking a formal class on the topic.

1. If you just want to do an IIR high or lowpass avoid filter orders of 10 or more. Revise the filter spec or change the filter type as needed.
2. Inspect your pole locations. Make sure your you have enough margin between the poles and the unit circle given your application requirements, data types, and implementation details.
3. Always implement IIR filters as cascaded second order section with proper zero-pole pairing and section ordering. Avoid MATLAB's filter() function or Python's scipy.lfilter(). Use sosfilt() or related flavors instead
4. Be aware that filtfilt() doubles the filter order, the passband ripple, the stopband ripple and moves the corner frequency. Make sure revise your base filter specification accordingly.
5. Avoid using filtfilt() directly. It's not documented which specific filter algorithm the function uses internally and it may change over time. Instead, implement it manually using sosfilt()' and flip()`
6. If your signal has any DC bias, remove it first (unless you need it), OR initialize the filter state with a proper DC solution.

If we really wanted to implement your super high order Butterworth manually this would look like

%% High order Butterworth using second order sections
x0 = ecg_samples(1:3600); % single input
x0 = x0 - mean(x0);  % remove DC bias

% BAD filter design spec
hhp = designfilt("highpassiir", ...
    'PassbandFrequency', fpass_hhp, ...
    'StopbandFrequency', fstopband_hhp, ...
    'PassbandRipple', 1, ...
    'StopbandAttenuation', 60, ...
    'SampleRate', fc, ...
    'DesignMethod', "butter");
 % get second order sections from filter
 sos =  hhp.Coefficients;
 % forward filter
 x1 = sosfilt(sos,x0);
 % backwards filter
 x2 = flip(sosfilt(sos,flip(x1)));
 plot(x2);

Answer 2

Note that the designed Butterworth filter is of order 47 while the designed elliptic filter just has a lower order 7 given the same specification which is a known fact as mentioned here.

Keeping the same spec, you can incorporate the modification below to get a similar result after zero-phase filtering by the designed Butterworth filter

[z,p,k] = zpk(hhp_butter);
[sos,g] = zp2sos(z,p,k,"down");
ecg_samples_filtered1 = filtfilt(sos,g,[ecg_samples(end:-1:1),ecg_samples,ecg_samples(1:end)]);
ecg_samples_filtered1 = ecg_samples_filtered1(length(ecg_samples)+(1:length(ecg_samples)));

Here is the result