Applying an suitable butterworth filter on raw signal using Python

Question

I have acquired a 10 seconds raw PPG (Photoplethysmogram) signal from my TI AFE4490. My hardware is calibrated and I'm using 250 samples per second to record those signal. I acquired 2500 points at the end.

I used Butterworth bandpass filter with lowcut=0.5 , highcut=15 and order=2 . You can see my raw and filtered signals bellow:

I also tried to filter this using a Butterworth lowpass filter with lowcut=15 and order=2 . As you can see, my raw and filtered signals are bellow:

I read at some articles that 0.5Hz and 15Hz are the good lowcut and highcut frequencies to this type of signal.

Before I apply the filters I used an Scipy Butterworth (from scipy docs ) algorithm to show me the filter response, and it was good.

My filtered signal seems to be good after that "start" (elevation at the beginning), but I don't know why that beginning. Anyone can tell me if that "start" it's normal at Butterworth filters? If yes, there is some method to fix it?

I appreciate your help.

My code:

RED, IR, nSamples, sRate = getAFESignal()

period = 1/sRate # Signal period.
# Desired cutoff frequency (in Hz) and filter order.
lowcut = 0.5
highcut = 15
orders = 2

plt.figure(1)

x = np.linspace(0, nSamples*period, nSamples, endpoint=True)

plt.subplot(2,1,1)
y = IR
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x,y, label='Noisy signal')


plt.subplot(2,1,2)
yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x, yf, label='Filtered signal')

plt.grid()
plt.show()

The function getAFEsignal() is just a function to read a .txt file and put all into two numpy arrays.

SleuthEye SleuthEye · Accepted Answer · 2018-08-28T03:33:40

The initial transient that you see in your graphs is the filter's step response as a sudden input is applied to the filter in its resting state. If you had just connected a physical instrument including such a bandpass filter, the instrument's sensor might have picked up input data samples jumping from 0 (while the probe is disconnected) to the first sample's value of ~0.126V. The response of the instrument's filter would have then shown a similar transient.

However, you are probably more interested in the steady-state response of the instrument after it is no longer disturbed by these external factors (such as the probe being connected), and has had time to settle to the properties of the signal of interest.

One way to achieve this is to use a data sample that is sufficiently long and discard the initial transient. Another approach is to force the filter's initial internal state to something close to what might be expected had a signal of similar amplitude been applied for some time before your first input sample. This can be done for example by setting the initial condition with scipy.signal.lfilter_zi.

Now, I assume you've used butter_bandpass_filter from SciPy Cookbook, which doesn't take care of the filter's initial conditions. Fortunately, it can be easily modified to this end:

from scipy.signal import butter, lfilter, lfilter_zi

def butter_bandpass_filter_zi(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    zi = lfilter_zi(b, a)
    y,zo = lfilter(b, a, data, zi=zi*data[0])
    return y

Another thing to note at this point is that you are calling butter_bandpass_filter as:

yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)

passing nSamples (the total number of samples, in your case 2500) as the 4th argument, whereas the function is expecting the sampling rate (in your case 250). The factor of 10 between the two quantities has an effect equivalent to reducing the filtered range from [0.5,15] Hz to [0.05,1.5] Hz. To get the intended bandpass frequency range, you should pass sRate as the fourth argument:

yf = butter_bandpass_filter_zi(IR, lowcut, highcut, sRate, order=orders)

Finally, you may notice that this last output is a little less triangular than the input. This is caused by some of the low-frequency content near 0.5Hz being filtered out. If that is what you were expecting then great. Otherwise, you can still play around with the cutoff frequencies to get whatever you feel yields the best results. For example (and I do not mean to say that this is a better frequency range), if you were to set lowcut=0.25 you would get a more triangular graph such as:

Applying an suitable butterworth filter on raw signal using Python

1 Answers