Camera measurement of heart rate - frequency bias towards 50 bpm (0.83 Hz)?

Question

I am trying to develop a RGB camera-based non-contact measurement of heart rate using principles of remote photoplethysmography. I want to be able to measure a range of heart rates I have been looking into multiple methods, which are all quite similar though. Typically, they all detect/track the face, and then they spatially average the colors in the facial region. They either do a FFT of the green or red temporal signal, or they do the FFT of a component coming from a blind-source separation method (like ICA or PCA). Now, I have been working on this, however I have noticed that the frequency seems to always stay around 50bpm (my heart rate is usually around 70bpm), occasionally jumping to the correct bpm. I have even checked this with a fitness tracker. I have tried adding filters, detrending, etc. but with no success. What could be the problem here? Could there be some sort of other frequencies I am not considering here from the environment? I have tried this out in various environments, and illuminations and that does not seem to affect the problem. Is FFT not robust enough for this method? I have tried Welch's PSD, which seems like it gives better results, but still has some bias towards 50 bpm (0.83 Hz).

Here is some code based on this paper:

import cv2
import datetime
import numpy as np
from scipy import signal
from scipy.fftpack import fft, fftfreq, fftshift
from sklearn.decomposition import PCA, FastICA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
# Change these variables based on the location of your cascade classifier
PATH_TO_HAAR_CASCADES = "..." 
face_cascade = cv2.CascadeClassifier(PATH_TO_HAAR_CASCADES+'haarcascade_frontalface_default.xml') # Full pathway must be used
firstFrame = None
time = []
R = []
G = []
B = []
pca = FastICA(n_components=3)
cap = cv2.VideoCapture(0)
if cap.isOpened() == False:
    print("Failed to open webcam")
frame_num = 0
plt.ion()
while cap.isOpened():
    ret, frame = cap.read()
    if ret == True:
        frame_num += 1
        if firstFrame is None:
            start = datetime.datetime.now()
            time.append(0)
            # Take first frame and find face in it
            firstFrame = frame
            cv2.imshow("frame",firstFrame)
            old_gray = cv2.cvtColor(firstFrame, cv2.COLOR_BGR2GRAY)
            faces = face_cascade.detectMultiScale(old_gray, 1.3, 5) 
            if faces == ():
                firstFrame = None
            else:
                for (x,y,w,h) in faces: 
                    x2 = x+w
                    y2 = y+h
                    cv2.rectangle(firstFrame,(x,y),(x+w,y+h),(255,0,0),2)
                    cv2.imshow("frame",firstFrame)
                    VJ_mask = np.zeros_like(firstFrame)
                    VJ_mask = cv2.rectangle(VJ_mask,(x,y),(x+w,y+h),(255,0,0),-1)
                    VJ_mask = cv2.cvtColor(VJ_mask, cv2.COLOR_BGR2GRAY)
                ROI = VJ_mask
                ROI_color = cv2.bitwise_and(ROI,ROI,mask=VJ_mask)
                cv2.imshow('ROI',ROI_color)
                R_new,G_new,B_new,_ = cv2.mean(ROI_color,mask=ROI)
                R.append(R_new)
                G.append(G_new)
                B.append(B_new)

        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

        else:
            current = datetime.datetime.now()-start
            current = current.total_seconds()
            time.append(current)
            frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
            ROI_color = cv2.bitwise_and(frame,frame,mask=ROI)
            cv2.imshow('ROI',ROI_color)
            R_new,G_new,B_new,_ = cv2.mean(ROI_color, mask=ROI)
            R.append(R_new)
            G.append(G_new)
            B.append(B_new)
            if frame_num >= 900:
                N = 900
                G_std = StandardScaler().fit_transform(np.array(G[-(N-1):]).reshape(-1, 1))
                G_std = G_std.reshape(1, -1)[0]
                R_std = StandardScaler().fit_transform(np.array(R[-(N-1):]).reshape(-1, 1))
                R_std = R_std.reshape(1, -1)[0]
                B_std = StandardScaler().fit_transform(np.array(B[-(N-1):]).reshape(-1, 1))
                B_std = B_std.reshape(1, -1)[0]
                T = 1/(len(time[-(N-1):])/(time[-1]-time[-(N-1)]))
                X_f=pca.fit_transform(np.array([R_std,G_std,B_std]).transpose()).transpose()
                N = len(X_f[0])
                yf = fft(X_f[1])
                yf = yf/np.sqrt(N)
                xf = fftfreq(N, T)
                xf = fftshift(xf)
                yplot = fftshift(abs(yf))
                plt.figure(1)
                plt.gcf().clear()
                fft_plot = yplot
                fft_plot[xf<=0.75] = 0
                print(str(xf[fft_plot[xf<=4].argmax()]*60)+' bpm')
                plt.plot(xf[(xf>=0) & (xf<=4)], fft_plot[(xf>=0) & (xf<=4)])
                plt.pause(0.001)

Some notes about the code:

1. Make sure to put the location of your cascade classifier in the code

2. While measuring, keep still as it only detects the region of interest in the first frame

EDIT (4/9/2018): Is this post better suited for a better website, like Stack Overflow? I have seen several heart rate questions on StackOverflow, but since I think the problem is more with the signal processing of the video, I thought it would fit here more.

EDIT (4/11/2018): Here are some figures of the signals as requested by @A_A. This data was obtained with the same program but with a frame rate of 20 fps. However, the data and result is similar for 30fps signal as well.

Normalized RGB signals (color is the signal from that color channel):

After the RGB signals are passed into the ICA, the components (filtered with a 4th order Butterworth bandpass):

Finally, here is the FFT of the signal. The maximum peak is what is marked as the heart rate. Note that lower frequencies below 0.75 Hz (45 bpm) are already discarded.

EDIT (4/12/2018): I redid the plots, and normalized the FFT, which somehow made the FFT clearer. Also, I did a plot with the no face and just background. The peak in the RGB signals might be because I accidentally came into the field of view.

Answer 1

This is an interesting journey, which is detailed further here.

I was intrigued by the attempt to derive a Photoplethysmography (PPG) signal just from the visible spectrum as it is sensed by a common web camera. An indicative publication can be found here (link shared in chat).

I have to say, I am still a bit skeptical about this and I am thinking that any PPG related signal originating from a web cam is largely due to its sensor's sensitivity to the IR spectrum. This of course, can be easily checked and this has been part of my interest too.

If the PPG signal is not due to the IR then the statement that this work would confirm is "The web cam is capable of picking up the subtle pulsation of the human skin due to the circulation of blood".

A web cam produces an RGB image by combining three grayscale images formed (roughly) by sampling different regions (Red, Green, Blue) of the optical spectrum. Each grayscale image has a colour depth of 8 bits and produces a range of 256 intensity values. A common web cam is capable of capturing light flickering due to the mains AC (albeit aliased) that the human eye cannot observe. However, a pulsating redness around the face of a person is at low enough frequency to be observed but it goes unnoticeable. This means that it is not intense enough and would therefore occupy a limited range of those 256 values. This means that even if some signal is there, everything else within the web cam's Field of View (FOV) is noise. This includes AC lighting flicker which is so noticeable in digital cameras for manufacturers to have included a compensation feature in their products. That's a very difficult environment from the point of view of Signal To Noise Ratio (SNR).

In what follows, I will try to summarise the key points that have emerged from our discussion on the chat and provide a few comments in terms of answers and recommendations. So, in no particular order, here we go:

The signal pathway as it was shared by @TanMath in chat is this: . Based on the first few graphs that were posted here, I suspected that the system suffered from aliasing and that the bandpass filter was ringing.

Aliasing: The spectrum seems to have a "segmented" structure, with dips every few components. A similar structure appears when you modulate 100% by a square wave but much more pronounced. Here it was not as bad but still there.

Ringing: This has also been a question from @TanMath, "What do you mean by ringing?". Very briefly, the bandpass filter's pass band is between 0.6-1.5 Hz and the filter is of order 4. Narrow filters of high order tend to "ring". A Finite Impulse Response (FIR) filter's output is $y = x \circledast h$ where $\circledast$ denotes convolution and $x$ is the input signal. A high order bandpass means that $h$ is long and tends to look like a sinusoid. If $x$ contains sudden changes or impulses, this means that $y$ will contain a large proportion of $h$ as the impulse excites it and its values fade. A similar consideration applies in Infinite Impulse Response (IIR) with the extreme example of the resonator which is an extremely selective bandpass filter. Extremely selective it might be but at that region it is so close to the unit circle that it might even be dipping in and out of being an oscillator. Sustaining its oscillations without any input.

So, we proceeded to run quite a few experiments to understand the behaviour of the system and it slowly dawned that the way the data is acquired means that in addition to the challenging environment we are operating in, we make things difficult by not sampling regularly.

This pastebin (shared in chat) contains the timestamps of each video frame as obtained by the web cam and it tells a horrific story. Here is a simple plot of it:

The x-axis is the frame number (1,2,3,4,...) and the y-axis is the physical time it corresponds to. On first examination, there are two different rates present in it as indicated by two different slopes. But, if you differentiate the time vector, here is what you get:

Now, ideally, this should have been a flat line denoting the same time difference between two successive frames or a constant frame rate. But that is not what we see here. What we observe is about 7-8 different discrete timings all co-existing. Now, if you take the Discrete Fourier Transform of this, you get something like this:

This is the frequency spectrum of the sampling period. Ideally, if the frame-rate was constant, you would get just the DC component ($f = 0$ Hz) but obviously here the sampling period is not constant which means that anything in the input of your system is effectively convolved with this spectrum or in other words, it gets modulated by this spectrum.

This is why I recommended that a pulsating phantom is created, an animated GIF at a known frequency. This, along other phantoms, were created and imaged and they are available here. The sinusoidal one looks like this:

There are a few things visible here, which I will come back to, but one thing you can see here is the sinusoid "skipping phase", because of the jitter in the sampling. Thankfully, the timestamp VS RGB values signal was also available which meant that we could re-interpolate to a regular grid and obtain the spectrum of the regular waveform which looked like this:

Which you might think that it is expectable, given the earlier plots of the spectrum of the sampling frequency but if you look closely:

You can now see the "segmented nature" I am talking about. Parts of the spectrum are almost symmetrical and there are strong peaks not only at the fundamental frequency we are trying to measure.

Now, I have to say, what we see here is a portion of the sinusoid that does not include the transient, has its mean subtracted and is resampled at the equivalent of 30 Frames Per Second (FPS). This was achieved with plain simple linear interpolation like this:

Given $t,R$ are two vectors containing the timestamp and light intensity read off the camera at that light intensity:

U = 0:(1./30):(33.367-(1./30)); # Regular sampling at 30FPS and equivalent duration    
V = interp1(t,R,U); # V now contains the interpolated waveform

$V$ now contains an interpolated waveform which can be treated like a regularly sampled waveform without explicit reference to its vector of timestamps. (Other interpolation schemes available in interp1 might work better for you, for example cubic or spline handles that change in $Fs$ towards the end better ).

What is this transient? The transient is the camera adjusting exposure automatically, compensating for flicker and who knows what else is encoded in its driver and beyond our control. This is important for this application because automatic adjustment of exposure is a closed loop automatic control system (similar to Automatic Gain Control) and it is basically a "normalisation" step that happens in hardware. If you let the camera stabilise and then turn on a spot light in its face, it will adjust to this flood of light by reducing the exposure to make sure that the sensor doesn't clip. This is visible in the "step response" phantom we obtained, here and it looks like this:

Effectively, the exposure control of the camera is adjusting the SNR.

So, is this salvageable?

I cannot say with certainty if it will be possible to measure a reliable PPG related signal in general but it will be impossible to decide with this system. Therefore, my recommendations here are a first step towards something more reliable which could provide some indication that a signal is there.

I am happy to amend this post with stage 2 recommendations covering the use of Independent Component Analysis (which inserts a whole different can of worms here) and the system running on the Raspberry Pi after we have something more reliable to work with. Until then:

1. Remove all frame processing code from the acquisition routine. It is unlikely that the processing required for the estimation of the heart rate finishes before the acquisition of the next frame and in any case, it is impossible to guarantee that this routine takes exactly the same amount of time for each frame. To achieve this, you need to move towards a multithreaded architecture. Create a thread whose only job is to grab a frame and its timestamp and have these two being pushed to a deque. A deque is a First In Fist Out (FIFO) stack which is also thread safe. In other words, you don't have to add any other special control structures to coordinate the access of your data acquisition thread and your processing thread to the same data structure. You can have your data acquisition thread pushing values to the deque on one end and your processing thread pop values off of it from the other. There are various ways by which you can achieve the triggering of a function at a predefined rate in python including Timers, Events or a plain simple sleep within the separate thread. In anything other than a Real Time operating system, these methods are still unreliable. What we are trying to do here is move towards the lowest jitter possible.

2. Despite the use of more regular sampling, make sure that you make use of the frame timestamp as much as possible. (For example, to resample a block of acquired samples at a more regular rate).

3. Turn off all compensation features of the camera including automatic exposure control and flicker compensation. This means that you are going to have to control the lighting very carefully. If you notice, the authors of the cited paper are taking measuements of faces in a plain white background under natural light. I think that artificial lighting will creep in to the video feed which you probably have already experienced. It is impossible to control aliasing in that case.

4. You are dealing with a very difficult SNR environment. This is visible from the constant phantom results here. This is the response of the system for no input. It is still subject to the automatic exposure control of the camera but it gives you an indication of the kind of natural background "noise" you can expect. All components hover around the 100 level WITHOUT any input. The question then becomes, how much variation is expected from the PPG related signal and how can we maximise that variation and suppress the noise. Along this line of thinking, I would recommend that you work on the signal that is the result of the average of the three colours $\frac{R+G+B}{3}$.

Hope this helps.

Answer 2

Do not forget to consider parasite frequences from the environment, for instance the respiratory rate, the illumination changes due to the electricity network (50-60Hz). You can use a band-pass filtering that keep frequences between 40bpm and 200bpm for instance.

Also, you may not consider all color chanels. As the blood absorbs green radiation (the complementary color of red), it could be more relevant to observe only the green channel (or the total illumination i.e $1/3 * (R + G + B)$).

You can also try to use auto-correlation on the signal and then do a peak detection to find frequences.

One last point, as the blood flow is not identical on all the face, you have to avoid averaging the color. Consider taking measure points on the most relevants parts of the face (cheeks, forehead), then applying your method to each measure and finally to average the results.

You can also take into account other signals like small forehead movements (thesis here).