Detecting Local Minima in a 1D array of points

Question

For context, this graph represents the EAR (Eye Aspect Ratio) of my eye throughout a video, and the clear steep drops represent blinks. When I blink, the EAR rapidly decreases and increases again, since EAR measures how open my eye is. Given a 1D array of points, how would I detect the local minima that are representative of my blinks? In this example, it would detect all 5 blinks.

I tried this:

import numpy as np
from scipy.signal import find_peaks

# Invert the data to find local minima
inverted_data = -np.array(seq)

# Use find_peaks from SciPy to find the peaks in the inverted data
# Adjust the 'distance' parameter as needed for your dataset
peaks, _ = find_peaks(inverted_data, distance=5, height=-0.2)

print(peaks)

# The number of local minima is the length of the peaks array
number_of_minima = len(peaks)

print("Number of local minima:", number_of_minima)

But it didn't work very well because although the distance parameter is helpful, and potentially a good start, the height parameter is representative of the minimum EAR it takes to be considered a blink, and this varies for different eyes.

Answer 1

Demonstrating a very simple band-pass filter with some bogus data; comments inline:

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal

# https://en.wikipedia.org/wiki/Blinking#Adults
# Assuming adult, object focus. Delete all of these when you have the real sample rate.
blink_rate = 4/60
blink_sample_period = 35  # Eyeballed (!) from OP plot
samples = 220             # Eyeballed from OP plot

# Bogus; replace me with the real value
sample_rate = blink_rate * blink_sample_period

# Assume uniform time. Replace this if you know the actual sample times.
time = np.arange(0, samples/sample_rate, 1/sample_rate)

# Bogus, random "subsequence"
rand = np.random.default_rng(seed=0)
noise = rand.normal(scale=0.01, size=time.size)
blinks = 0.1*(0.5 + 0.5*np.cos(time*blink_rate * 2*np.pi + 0.5))**25
slope = 0.29 + (0.23 - 0.29)/time.size * time
seq = slope - blinks + noise

# Simple post-facto (not streamed) FFT bandpass filter
fft_input = np.fft.rfft(a=seq, norm='forward')
fft_freq = np.fft.rfftfreq(n=seq.size, d=1/sample_rate)
locut = 0.05  # Hz
hicut = 0.5   # Hz
fft_output = fft_input.copy()
fft_output[:round(locut*seq.size / sample_rate)] = 1e-5  # For log display. Just set these to 0
fft_output[round(hicut*seq.size / sample_rate):] = 1e-5
filtered = -np.fft.irfft(a=fft_output, norm='forward')

peaks, props = scipy.signal.find_peaks(
    x=filtered,
    prominence=0.5*filtered.max(),
    distance=5,  # samples
)

fig, (time_ax, freq_ax) = plt.subplots(ncols=2)
fig.suptitle('Eye aspect ratio blinking during video focus')

time_ax.plot(time, seq, label='input')
time_ax.plot(time, filtered, label='filtered')
time_ax.scatter(time[peaks], filtered[peaks], label='blink')
idx_ax = time_ax.secondary_xaxis('top', functions=(
    lambda t: t * sample_rate,
    lambda i: i / sample_rate,
))
idx_ax.set_xlabel('index')
time_ax.set_xlabel('time (s)')
time_ax.set_ylabel('Eye Aspect Ratio')
# time_ax.legend()

freq_ax.semilogy(fft_freq, np.abs(fft_input), label='input')
freq_ax.semilogy(fft_freq, np.abs(fft_output), label='output')
idx_ax = freq_ax.secondary_xaxis('top', functions=(
    lambda f: f*seq.size / sample_rate,
    lambda i: i/seq.size * sample_rate,
))
idx_ax.set_xlabel('index')
freq_ax.set_xlabel('freq (Hz)')
freq_ax.set_ylabel('FFT amplitude')
freq_ax.legend()

plt.show()