pythonmathgraphing

Detecting Local Minima in a 1D array of points


Graph in 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.


Solution

  • 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()
    

    filtered output