Homework 5

Due: Mar 4

Submit your code on OWL

Signal Processing Exercises

Here is some Python code to generate a complex signal composed of sinusoids at several frequencies. The signal is sampled at 1000 Hz.

import numpy as np
t = np.arange(0, 0.5, 1/1000)  # 500 ms of time at 1000 Hz
y1 = np.sin(2*np.pi*30 * t)    # pure sinusoid at  30 Hz
y2 = np.sin(2*np.pi*125 * t)   # pure sinusoid at 125 Hz
y3 = np.sin(2*np.pi*300 * t)   # pure sinusoid at 300 Hz
y = y1 + (0.8*y2) + (0.9*y3)

Plot the signal y in the time domain. Label the axes and include a title.
Plot the magnitude spectrum of y using the scipy.signal.welch() function. Label the axes and include a title.
Highpass filter the signal to remove the 30 Hz component. Plot the filtered signal in the time domain and also plot the magnitude spectrum of the filtered signal. Label the axes and include a title.
Lowpass the filtered signal from 3. to remove the 300 Hz component. Plot the filtered signal in the time domain and also plot the magnitude spectrum of the filtered signal. Label the axes and include a title.

Estimate your heart rate using your phone’s camera

Estimate your heart rate using a video taken with your webcam or your smartphone (or if you don’t have access to one, you can analyse my video instead, linked below).

Here’s what you will do to generate the data. First, make sure you’re in a brightly lit space (i.e. turn on the lights if you’re in a dark room). Better yet, find a light source like a desk lamp or a window, and stand in front of it with your webcam or smartphone camera. Take a minute or two to sit down and relax, to get your heart rate near its resting rate (usually around 60 beats per minute). Now gently rest the tip of your index finger against the lens of your webcam or smartphone camera. Now record a video. Record for somewhere between 15 and 30 seconds. Try not to move your finger, keep it as still as you can for the duration of the video.

If you used your smartphone, transfer the video file to your computer. Open it up and look at it on the screen. Here’s a sample frame from mine, which was approximately 15 seconds in duration (sample rate 24.0 frames per second):

Here is a link to the complete video (you can right-click and save as… to download it):

The idea here is that as your heart pumps blood through your arteries and veins, the amount of light that the camera detects, as it comes through your fingertip, will be modulated by the differences in blood flow. The variation in brightness will be quite subtle, perhaps even not visible with the naked eye, and so your task is to use signal processing techniques to try to detect it.

Here are some suggestions as to how to proceed. First load in the video file. Each frame of the video will have three channels: red, green and blue. Each frame will contain the brightness of each pixel for each of the three channels red, green and blue. You might as well average over all the pixels in the frame to get an average brightness measure. The signal that you ultimately want is probably the brightness of the red channel over time (i.e. for each frame of the video).

Here are some other things to consider:

make sure to use the correct sampling rate (frame rate) of your video file when doing signal analysis or filtering
in what frequency range do you expect your heart rate to reside?

Your tasks

Compute the average pixel brightness, averaged over all pixels, for the red channel, for each frame. This gives you a time series, one value for each frame of the video. Subtract the mean of the time series from each value. This will center the time series around zero. Plot the resulting time series on the y-axis, and time in seconds on the x-axis.
Band-pass filter the zero-centered time series from task 1 using cutoff frequencies to bracket what you think is likely to be your own heart rate during the recording. Remember, cutoff frequencies are defined in Hz (cycles per second) whereas heart rate is typically characterized as beats per minute. So a 60 beats per minute resting heart rate corresponds to 1.0 Hz. Plot the filtered time series.
Compute and plot a spectrum of the filtered time series. Identify the frequency corresponding to the peak of the spectrum and convert the value to beats per minute.

Hints

you will probably need to install the packages imageio and imageio-ffmpeg either using conda install or using pip or whatever Python package manager you use
then in your Python code you can use the imageio package to load in the video file like so:

import imageio
filename = "images/IMG_1383.MOV"
vid = imageio.get_reader(filename,  'ffmpeg')
fps = vid.get_meta_data()['fps']
print(f"read video at {fps} fps")

read video at 24.0 fps

read in frames from the video like so:

n = 0
for image in vid.iter_data():
    n += 1 # count frames
red = np.zeros(n)
i = 0
for image in vid.iter_data():
    red[i] = np.mean(image[:,:,0]) # the third dimension codes red/green/blue channels
    i += 1

Sample plot

When I do this for my video (linked above), band-pass filtering between 40 and 80 beats per minute, I get the following:

You don’t need to get the exact same thing. The important thing is that you can identify the frequency corresponding to the peak of the spectrum, and convert that to beats per minute.