Lowpass Filter for Heart Rate Signal in Python

General Discussions
, ,
1

I am working on filtering a noisy heart rate signal using Python. Since heart rates should not exceed around 220 beats per minute, I want to filter out all noise above 220 bpm. I first converted 220 bpm into 3.66666666 Hertz and then to rad/s, which gives me 23.0383461 rad/sec. The sampling frequency of my data is 30Hz, so I converted that to rad/s to get 188.495559 rad/s.

After researching online, I found some bandpass filter code and adapted it into a lowpass filter. Here’s the code I came up with:

from scipy.signal import butter, lfilter
from scipy.signal import freqs

def butter_lowpass(cutOff, fs, order=5):
    nyq = 0.5 * fs
    normalCutoff = cutOff / nyq
    b, a = butter(order, normalCutoff, btype='low', analog=True)
    return b, a

def butter_lowpass_filter(data, cutOff, fs, order=4):
    b, a = butter_lowpass(cutOff, fs, order=order)
    y = lfilter(b, a, data)
    return y

cutOff = 23.1  # cutoff frequency in rad/s
fs = 188.495559  # sampling frequency in rad/s
order = 20  # filter order

y = butter_lowpass_filter(data, cutOff, fs, order)
plt.plot(y)

I’m confused because I’m pretty sure the butter function expects the cutoff and sampling frequency in rad/s, but the output seems odd. Is it actually expecting these values in Hz?

Also, I don’t fully understand the purpose of these two lines:

nyq = 0.5 * fs
normalCutoff = cutOff / nyq

I know this is related to normalization, but I thought the Nyquist frequency was 2 times the sampling frequency, not 0.5. Why is the Nyquist used as a normalizer?

Can someone explain how to properly create a low pass filter in Python with these functions?

I also plotted the frequency response of the filter using the following code:

w, h = signal.freqs(b, a)
plt.plot(w, 20 * np.log10(abs(h)))
plt.xscale('log')
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green')  # cutoff frequency
plt.show()

However, the plot doesn’t show the expected cutoff at 23 rad/s. Can someone help clarify this?

2

Sure, I’ve worked with low pass filters in Python for biomedical signals like this! Your issue seems to arise from confusion between rad/s and Hz. SciPy’s butter function actually expects frequencies in Hz, not rad/s. Here’s how to adjust:

  1. Convert Frequencies from Rad/s to Hz:
# Convert rad/s to Hz
cutOff_Hz = cutOff / (2 * np.pi)  # 23.0383461 rad/s to Hz
fs_Hz = fs / (2 * np.pi)  # Convert fs from rad/s to Hz
  1. Pass the Values in Hz to the Filter:
# Now call the butter function with frequencies in Hz
b, a = butter(order, cutOff_Hz / (0.5 * fs_Hz), btype='low', analog=False)

This ensures you’re working in the correct units for digital signal processing. Adjusting this should resolve the discrepancy you’re seeing in the filter’s behavior.

3

Hey there! I’ve done some digital signal processing in the past, so let me clarify the part about the Nyquist frequency and normalization.

  1. Nyquist Frequency: The Nyquist frequency is half of the sampling frequency, not double. This is a fundamental concept in signal processing: the maximum frequency that can be sampled without aliasing is half the sampling rate. That’s why the formula nyq = 0.5 * fs is used.
  2. Normalization: To use the butter function properly, the cutoff frequency is normalized relative to the Nyquist frequency:
normalCutoff = cutOff / nyq

This normalization ensures that the cutoff frequency is expressed as a fraction of the Nyquist frequency (between 0 and 1), as required for digital filters.

To recap:

  • nyq = 0.5 * fs: The Nyquist frequency is half the sampling rate.
  • normalCutoff = cutOff / nyq: This normalizes the cutoff frequency.

Let me know if you need further clarification!

4

Hi, I’ve worked on several biomedical signal projects, so I totally get where your confusion might come from. Let’s also talk about the analog vs. digital filter settings in your code.

  1. Digital vs. Analog Filters: The butter function allows you to create both analog and digital filters. Since you’re working with sampled data (like your heart rate signal), you need a digital filter. Right now, you’re using analog=True, which might be the issue.
  2. Fixing for Digital Filtering: Switch to a digital filter by setting analog=False:
b, a = butter(order, normalCutoff, btype='low', analog=False)
  1. Frequency Response: Once you fix the units and ensure you’re using a digital filter, the frequency response plot should match your expectations. If it still doesn’t, double-check the normalization and conversion from rad/s to Hz.

To summarize:

  • Use analog=False for sampled signals.
  • Double-check your conversions and normalization.

Hope this clears up the remaining confusion!