Let’s improve on the emotion recognition from a previous article about FisherFace Classifiers. We will be using facial landmarks and a machine learning algorithm, and see how well we can predict emotions in different individuals, rather than on a single individual like in another article about theemotion recognising music player.

Important: The code in this tutorial is licensed under the GNU 3.0 open source license and you are free to modify and redistribute the code, given that you give others you share the code with the same right, and cite my name (use citation format below). You are not free to redistribute or modify the tutorial itself in any way. By reading on you agree to these terms. If you disagree, please navigate away from this page.

Citation format
van Gent, P. (2016). Emotion Recognition Using Facial Landmarks, Python, DLib and OpenCV. A tech blog about fun things with Python and embedded electronics. Retrieved from: http://www.paulvangent.com/2016/08/05/emotion-recognition-using-facial-landmarks/

IE users: I’ve gotten several reports that sometimes the code blocks don’t display correctly or at all on Internet Explorer. Please refresh the page and they should display fine.

Introduction and getting started
Using Facial Landmarks is another approach to detecting emotions, more robust and powerful than the earlier used fisherface classifier, but also requiring some more code and modules. Nothing insurmountable though. We need to do a few things:

• Get images from a webcam
• Detect Facial Landmarks
• Train a machine learning algorithm (we will use a linear SVM)
• Predict emotions

Those who followed the two previous posts about emotion recognition will know that the first step is already done.

Also we will be using:

• Python (2.7 or higher is fine, anaconda + jupyter notebook is a nice combo-package)
• OpenCV (I still use 2.4.9……so lazy, grab here)
• SKLearn (if you installed anaconda, it is already there, otherwise get it with pip install sklearn)
• Dlib (a C++ library for extracting the facial landmarks, see below for instructions)
• Visual Studio 2015 (get the community edition here, also select the Python Tools in the installation dialog).
  • Note that VS is not strictly required, I just build the modules against it. However it is a very nice IDE that also has good Python bindings and allows you to quickly make GUI applications to wrap around your Python scripts. I would recommend you give it a go

Installing and building the required libraries

I am on Windows, and building libraries on Windows always gives many people a bad taste in their mouths. I can understand why, however it’s not all bad and often the problems people run into are either solved by correctly setting PATH variables, providing the right compiler or reading the error messages and installing the right dependencies. I will walk you through the process of compiling and installing Dlib.

First install CMake. This should be straightforward, download the windows installer and install. Make sure to select the option “Add CMake to the system PATH” during the install. Choose whether you want this for all users or just for your account.

Download Boost-Python and extract the package. I extracted it into C:\boost but it can be anything. Fire up a command prompt and navigate to the directory. Then do:

1
2
3
4
5
6
7
bootstrap.bat #First run the bootstrap.bat file supplied with boost-python
#Once it finished invoke the install process of boost-python like this:
b2 install #This can take a while, go get a coffee
#Once this finishes, build the python modules like this
b2 -a --with-python address-model=64 toolset=mscv runtime-link=static #Again, this takes a while, reward yourself and get another coffee.

Once all is done you will find a folder named bin, or bin.v2, or something like this in your boost folder. Now it’s time to build Dlib.

Download Dlib and extract it somewhere. I used C:\Dlib but you can do it anywhere. Go back to your command prompt, or open a new one if you closed it, and navigate to your Dlib folder. Do this sequentially:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Set two flags so that the CMake compiler knows where to find the boost-python libraries 
set BOOST_ROOT=C:\boost #Make sure to set this to the path you extracted boost-python to!
set BOOST_LIBRARYDIR=C:\boost\stage\lib #Same as above
# Create and navigate into a directory to build into
mkdir build
cd build
# Build the dlib tools
cmake ..
#Navigate up one level and run the python setup program
cd ..
python setup.py install #This takes some time as well. GO GET ANOTHER COFFEE TIGER!

Open your Python interpreter and type “import dlib”. If you receive no messages, you’re good to go! Nice.

Testing the landmark detector
Before diving into much of the coding (which probably won’t be much because we’ll be recycling), let’s test the DLib installation on your webcam. For this you can use the following snippet. If you want to learn how this works, be sure to also compare it with the first script under “Detecting your face on the webcam” in the previous post. Much of the same OpenCV code to talk to your webcam, process the image by converting to grayscale, optimising the contrast with an adaptive histogram equalisation and displaying it is something we did there.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#Import required modules
import cv2
import dlib
#Set up some required objects
video_capture = cv2.VideoCapture(0) #Webcam object
detector = dlib.get_frontal_face_detector() #Face detector
predictor = dlib.shape_predictor("shape_predictor_68_face_landmarks.dat") #Landmark identifier. Set the filename to whatever you named the downloaded file
while True:
    ret, frame = video_capture.read()
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
    clahe_image = clahe.apply(gray)
    detections = detector(clahe_image, 1) #Detect the faces in the image
    for k,d in enumerate(detections): #For each detected face
        
        shape = predictor(clahe_image, d) #Get coordinates
            for i in range(1,68): #There are 68 landmark points on each face
                cv2.circle(frame, (shape.part(i).x, shape.part(i).y), 1, (0,0,255), thickness=2) #For each point, draw a red circle with thickness2 on the original frame
    cv2.imshow("image", frame) #Display the frame
    if cv2.waitKey(1) & 0xFF == ord('q'): #Exit program when the user presses 'q'
        break

This will result in your face with a lot of dots outlining the shape and all the “moveable parts”. The latter is of course important because it is what makes emotional expressions possible.

Note if you have no webcam and/or would rather like to try this on a static image, replace line #11 with something like frame = cv2.imread(“filename”) and comment out line #6 where we define the video_capture object. You will get something like:

 

my face has dots
people tell me my face has nice dots
experts tell me these are the best dots
I bet I have the best dots

Extracting features from the faces
The first thing to do is find ways to transform these nice dots overlaid on your face into features to feed the classifer. Features are little bits of information that describe the object or object state that we are trying to divide into categories. Is this description a bit abstract? Imagine you are in a room without windows with only a speaker and a microphone. I am outside this room and I need to make you guess whether there is a cat, dog or a horse in front of me. The rule is that I can only use visual characteristics of the animal, no names or comparisons. What do I tell you? Probably if the animal is big or small, that it has fur, that the fur is long or short, that it has claws or hooves, whether it has a tail made of flesh or just from hair, etcetera. Each bit of information I pass you can be considered afeature, and based the same feature set for each animal, you would be pretty accurate if I chose the features well.

How you extract features from your source data is actually where a lot of research is, it’s not just about creating better classifying algorithms but also about finding better ways to collect and describe data. The same classifying algorithm might function tremendously well or not at all depending on how well the information we feed it is able to discriminate between different objects or object states. If, for example, we would extract eye colour and number of freckles on each face, feed it to the classifier, and then expect it to be able to predict what emotion is expressed, we would not get far. However, the facial landmarks from the same image material describe the position of all the “moving parts” of the depicted face, the things you use to express an emotion. This is certainly useful information!

To get started, let’s take the code from the example above and change it so that it fits our current needs, like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import cv2
import dlib
detector = dlib.get_frontal_face_detector()
predictor = dlib.shape_predictor("shape_predictor_68_face_landmarks.dat")
def get_landmarks(image):
    detections = detector(image, 1)
    for k,d in enumerate(detections): #For all detected face instances individually
        shape = predictor(image, d) #Draw Facial Landmarks with the predictor class
        xlist = []
        ylist = []
        for i in range(1,68): #Store X and Y coordinates in two lists
            xlist.append(float(shape.part(i).x))
            ylist.append(float(shape.part(i).y))
            
        for x, y in zip(xlist, ylist): #Store all landmarks in one list in the format x1,y1,x2,y2,etc.
            landmarks.append(x)
            landmarks.append(y)
    if len(detections) > 0:
        return landmarks
    else: #If no faces are detected, return error message to other function to handle
        landmarks = "error"
        return landmarks

Here we extract the coordinates of all face landmarks. These coordinates are the first collection of features, and this might be the end of the road. You might also continue and try to derive other measures from this that will tell the classifier more about what is happening on the face. Whether this is necessary or not depends. For now let’s assume it is necessary, and look at ways to extract more information from what we have. Feature generation is always a good thing to try, if only because it brings you closer to the data and might give you ideas or alternative views at it because you’re getting your hands dirty. Later on we’ll see if it was really necessary at a classification level.

To start, look at the coordinates. They may change as my face moves to different parts of the frame. I could be expressing the same emotion in the top left of an image as in the bottom right of another image, but the resulting coordinate matrix would express different numerical ranges. However, the relationships between the coordinates will be similar in both matrices so some information is present in a location invariant form, meaning it is the same no matter where in the picture my face is.

Maybe the most straightforward way to remove numerical differences originating from faces in different places of the image would be normalising the coordinates between 0 and 1. This is easily done by: , or to put it in code:

1
2
xnorm = [(i-min(xlist))/(max(xlist)-min(xlist)) for i in xlist]
ynorm = [(i-min(ylist))/(max(ylist)-min(ylist)) for i in ylist]

However, there is a problem with this approach because it fits the entire face in a square with both axes ranging from 0 to 1. Imagine one face with its eyebrows up high and mouth open, the person could be surprised. Now imagine an angry face with eyebrows down and mouth closed. If we normalise the landmark points on both faces from 0-1 and put them next to each other we might see two very similar faces. Because both distinguishing features lie at the edges of the face, normalising will push both back into a very similar shape. The faces will end up looking very similar. Take a moment to appreciate what we have done; we have thrown away most of the variation that in the first place would have allowed us to tell the two emotions from each other! Probably this will not work. Of course some variation remains from the open mouth, but it would be better not to throw so much away.

A less destructive way could be to calculate the position of all points relative to each other. To do this we calculate the mean of both axes, which results in the point coordinates of the sort-of “centre of gravity” of all face landmarks. We can then get the position of all points relative to this central point. Let me show you what I mean. Here’s my face with landmarks overlaid:

First we add a “centre of gravity”, shown as a blue dot on the image below:

Lastly we draw a line between the centre point and each other facial landmark location:

Note that each line has both a magnitude (distance between both points) and a direction (angle relative to image where horizontal=0°), in other words, a vector.

But, you may ask, why don’t we take for example the tip of the nose as the central point? This would work as well, but would also throw extra variance in the mix due to short, long, high- or low-tipped noses. The “centre point method” also introduces extra variance; the centre of gravity shifts when the head turns away from the camera, but I think this is less than when using the nose-tip method because most faces more or less face the camera in our sets. There are techniques to estimate head pose and then correct for it, but that is beyond this article.

There is one last thing to note. Faces may be tilted, which might confuse the classifier. We can correct for this rotation by assuming that the bridge of the nose in most people is more or less straight, and offset all calculated angles by the angle of the nose bridge. This rotates the entire vector array so that tilted faces become similar to non-tilted faces with the same expression. Below are two images, the left one illustrates what happens in the code when the angles are calculated, the right one shows how we can calculate the face offset correction by taking the tip of the nose and finding the angle the nose makes relative to the image, and thus find the angular offset β we need to apply.

Now let’s look at how to implement what I described above in Python. It’s actually fairly straightforward. We just slightly modify the get_landmarks() function from above.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def get_landmarks(image):
    detections = detector(image, 1)
    for k,d in enumerate(detections): #For all detected face instances individually
        shape = predictor(image, d) #Draw Facial Landmarks with the predictor class
        xlist = []
        ylist = []
        for i in range(1,68): #Store X and Y coordinates in two lists
            xlist.append(float(shape.part(i).x))
            ylist.append(float(shape.part(i).y))
            
        xmean = np.mean(xlist) #Find both coordinates of centre of gravity
        ymean = np.mean(ylist)
        xcentral = [(x-xmean) for x in xlist] #Calculate distance centre <-> other points in both axes
        ycentral = [(y-ymean) for y in ylist]
        
        if xlist[26] == xlist[29]: #If x-coordinates of the set are the same, the angle is 0, catch to prevent 'divide by 0' error in function
            anglenose = 0
        else:
            anglenose = int(math.atan((ylist[26]-ylist[29])/(xlist[26]-xlist[29]))*180/math.pi) #point 29 is the tip of the nose, point 26 is the top of the nose brigde
        if anglenose < 0: #Get offset by finding how the nose brigde should be rotated to become perpendicular to the horizontal plane
            anglenose += 90
        else:
            anglenose -= 90
        landmarks_vectorised = []
        for x, y, w, z in zip(xcentral, ycentral, xlist, ylist):
            landmarks_vectorised.append(x) #Add the coordinates relative to the centre of gravity
            landmarks_vectorised.append(y)
            #Get the euclidean distance between each point and the centre point (the vector length)
            meannp = np.asarray((ymean,xmean))
            coornp = np.asarray((z,w))
            dist = np.linalg.norm(coornp-meannp)
            landmarks_vectorised.append(dist)
            #Get the angle the vector describes relative to the image, corrected for the offset that the nosebrigde has when the face is not perfectly horizontal
            anglerelative = (math.atan((z-ymean)/(w-xmean))*180/math.pi) - anglenose
            landmarks_vectorised.append(anglerelative)
            
    if len(detections) < 1: 
        landmarks_vestorised = "error" #If no face is detected set the data to value "error" to catch detection errors
    
    return landmarks_vectorised
    

That was actually quite manageable, no? Now it’s time to put all of the above together with some stuff from the first post. The goal is to read the existing dataset into a training and prediction set with corresponding labels, train the classifier (we use Support Vector Machines with linear kernel from SKLearn, but feel free to experiment with other available kernels such as polynomial or rbf, or other classifiers!), and evaluate the result. This evaluation will be done in two steps; first we get an overall accuracy after ten different data segmentation, training and prediction runs, second we will evaluate the predictive probabilities.

Déja-Vu All Over Again

The next thing we will be doing is returning to the two datasets from the original post. Let’s see how this approach stacks up.

First let’s write some code. The approach is to first extract facial landmark points from the images, randomly divide 80% of the data into a training set and 20% into a test set, then feed these into the classifier and train it on the training set. Finally we evaluate the resulting model by predicting what is in the test set to see how the model handles the unknown data. Basically a lot of the steps are the same as what we did earlier.

The quick and dirty (I will clean and ‘pythonify’ the code later, when there is time) solution based off of earlier code could be something like:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import cv2, glob, random, math, numpy as np, dlib, itertools
from sklearn.svm import SVC
__author__ = "Paul van Gent, 2016" #Please leave this line in
emotions = ["anger", "contempt", "disgust", "fear", "happiness", "neutral", "sadness", "surprise"] #Emotion list
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
detector = dlib.get_frontal_face_detector()
predictor = dlib.shape_predictor("shape_predictor_68_face_landmarks.dat") #Or set this to whatever you named the downloaded file
clf = SVC(kernel='linear', probability=True, tol=1e-3)#, verbose = True) #Set the classifier as a support vector machines with polynomial kernel
def get_files(emotion): #Define function to get file list, randomly shuffle it and split 80/20
    files = glob.glob("dataset\\%s\\*" %emotion)
    random.shuffle(files)
    training = files[:int(len(files)*0.9)] #get first 80% of file list
    prediction = files[-int(len(files)*0.1):] #get last 20% of file list
    return training, prediction
def get_landmarks(image):
    detections = detector(image, 1)
    for k,d in enumerate(detections): #For all detected face instances individually
        shape = predictor(image, d) #Draw Facial Landmarks with the predictor class
        xlist = []
        ylist = []
        for i in range(1,68): #Store X and Y coordinates in two lists
            xlist.append(float(shape.part(i).x))
            ylist.append(float(shape.part(i).y))
            
        xmean = np.mean(xlist) #Get the mean of both axes to determine centre of gravity
        ymean = np.mean(ylist)
        xcentral = [(x-xmean) for x in xlist] #get distance between each point and the central point in both axes
        ycentral = [(y-ymean) for y in ylist]
        if xlist[26] == xlist[29]: #If x-coordinates of the set are the same, the angle is 0, catch to prevent 'divide by 0' error in function
            anglenose = 0
        else:
            anglenose = int(math.atan((ylist[26]-ylist[29])/(xlist[26]-xlist[29]))*180/math.pi)
        if anglenose < 0:
            anglenose += 90
        else:
            anglenose -= 90
        landmarks_vectorised = []
        for x, y, w, z in zip(xcentral, ycentral, xlist, ylist):
            landmarks_vectorised.append(x)
            landmarks_vectorised.append(y)
            meannp = np.asarray((ymean,xmean))
            coornp = np.asarray((z,w))
            dist = np.linalg.norm(coornp-meannp)
            anglerelative = (math.atan((z-ymean)/(w-xmean))*180/math.pi) - anglenose
            landmarks_vectorised.append(dist)
            landmarks_vectorised.append(anglerelative)
    if len(detections) < 1: 
        landmarks_vectorised = "error"
    return landmarks_vectorised
def make_sets():
    training_data = []
    training_labels = []
    prediction_data = []
    prediction_labels = []
    for emotion in emotions:
        training, prediction = get_files(emotion)
        #Append data to training and prediction list, and generate labels 0-7
        for item in training:
            image = cv2.imread(item) #open image
            gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) #convert to grayscale
            clahe_image = clahe.apply(gray)
            get_landmarks(clahe_image)
            if landmarks_vectorised == "error":
                pass
            else:
                training_data.append(landmarks_vectorised) #append image array to training data list
                training_labels.append(emotions.index(emotion))
    
        for item in prediction:
            image = cv2.imread(item)
            gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
            clahe_image = clahe.apply(gray)
            get_landmarks(clahe_image)
            if landmarks_vectorised == "error":
                pass
            else:
                prediction_data.append(landmarks_vectorised)
                prediction_labels.append(emotions.index(emotion))
    return training_data, training_labels, prediction_data, prediction_labels   
accur_lin = []
for i in range(0,10):
    print("Making sets %s" %i) #Make sets by random sampling 80/20%
    training_data, training_labels, prediction_data, prediction_labels = make_sets()
    npar_train = np.array(training_data) #Turn the training set into a numpy array for the classifier
    npar_trainlabs = np.array(training_labels)
    print("training SVM linear %s" %i) #train SVM
    clf.fit(npar_train, training_labels)
    print("getting accuracies %s" %i) #Use score() function to get accuracy
    npar_pred = np.array(prediction_data)
    pred_lin = clf.score(npar_pred, prediction_labels)
    print "linear: ", pred_lin
    accur_lin.append(pred_lin) #Store accuracy in a list
print("Mean value lin svm: %.3f" %np.mean(accur_lin)) #Get mean accuracy of the 10 runs

Remember that in the previous post, for the standard set at 8 categories we managed to get 69.3% accuracy with the FisherFace classifier. This approach yields 84.1% on the same data, a lot better!

We then reduced the set to 5 emotions (leaving out contempt, fear and sadness), because the 3 categories had very few images, and got 82.5% correct. This approach gives 92.6%, also much improvement.

After adding the less standardised and more difficult images from google, we got 61.6% correct when predicting 7 emotions (the contempt category remained very small so we left that out). This is now 78.2%, also quite an improvement. This remains the lowest accuracy, showing that for a more diverse dataset the problem is also more difficult. Keep in mind that the dataset I use is still quite small in machine learning terms, containing about 1000 images spread over 8 categories.

Looking at features
So we derived different features from the data, but weren’t sure whether this was strictly necessary. So was this necessary? It depends! It depends on if doing so adds more unique variance related to what you’re trying to predict, it depends on what classifier you use, etc.

Let’s run different feature combinations as inputs through different classifiers and see what happens. I’ve run all iterations on the same slice of data with 4 emotion categories of comparable size (so that running the same settings again yields the same predictive value).

Using all of the features described so far leads to:
Linear SVM: 93.9%
Polynomial SVM: 83.7%
Random Forest Classifier: 87.8%

Now using just the vector length and angle:
Linear SVM: 87.8%
Polynomial SVM: 87.8%
Random Forest Classifier: 79.6%

Now using just the raw coordinates:
Linear SVM: 91.8%
Polynomial SVM: 89.8%
Random Forest Classifier: 59.2%

Now replacing all training data with zeros:
Linear SVM: 32.7%
Polynomial SVM: 32.7%
Random Forest Classifier: 32.7%

Now this is interesting! First note that there isn’t much difference in the accuracy of the support vector machine classifiers when using the extra features we generate. This type of classifier already preprocesses the data quite extensively. The extra data we generate does not contain much if any extra information to this classifier, so it only marginally improves the performance of the linear kernel, and actually hurts the polynomial kernel because data with a lot of overlapping variance can also make a classification task more difficult. By the way, this is a nice 2D visualisation of what an SVC tries to achieve, complexity escalates when adding one dimension. Now remember that the SVC operates in an N-dimensional space and try to imagine what a set of hyperplanes in 4, 8, 12, 36 or more dimensions would look like. Don’t drive yourself crazy.

Random Forest Classifiers do things a lot differently. Essentially they are a forest of decision trees. Simplified, each tree is a long list of yes/no questions, and answering all questions leads to a conclusion. In the forest the correlation between each tree is kept as low as possible, which ensures every tree brings something unique to the table when explaining patterns in the data. Each tree then votes on what it thinks the answer is, and most votes win. This approach benefits extensively from the new features we generated, jumping from 59.2% to 87.8% accuracy as we combine all derived features with the raw coordinates.

So you see, the answer you likely get when you ask any scientist a direct question holds true here as well: it depends. Check your data, think twice and try a few things.

The last that may be noticed is that, when not adding any data at all and in stead presenting the classifiers with a matrix of zeros, they still perform slightly above the expected chance level of 25%. This is because the categories are not identically sized.

Looking at mistakes
Lastly, let’s take a look at where the model goes wrong. Often this is where you can learn a lot, for example this is where you might find that a single category just doesn’t work at all, which can lead you to look critically at the training material again.

One advantage of the SVM classifier we use is that it is probabilistic. This means that it assigns probabilities to each category it has been trained for (and you can get these probabilities if you set the ‘probability’ flag to True). So, for example, a single image might be “happy” with 85% probability, “angry” with “10% probability, etc.

To get the classifier to return these things you can use its predict_proba() function. You give this function either a single data row to predict or feed it your entire dataset. It will return a matrix where each row corresponds to one prediction, and each column represents a category. I wrote these probabilities to a table and included the source image and label. Looking at some mistakes, here are some notable things that were classified incorrectly (note there are only images from my google set, the CK+ set’s terms prohibit me from publishing images for privacy reasons):

anger: 0.03239878
contempt: 0.13635423
disgust: 0.0117559
fear: 0.00202098
neutral: 0.7560004
happy: 0.00382895
sadness: 0.04207027
surprise: 0.0155695

The correct answer is contempt. To be honest I would agree with the classifier, because the expression really is subtle. Note that contempt is the second most likely according to the classifier.

anger: 0.0726657
contempt: 0.24655082
disgust: 0.06427896
fear: 0.02427595
neutral: 0.20176133
happy: 0.03169822
sadness: 0.34911036
surprise: 0.00965867

The correct answer is disgust. Again I can definitely understand the mistake the classifier makes here (I might make the same mistake..). Disgust would be my second guess, but not the classifier’s. I have removed this image from the dataset because it can be ambiguous.

anger: 0.00304093
contempt: 0.01715202
disgust: 0.74954754
fear: 0.04916257
neutral: 0.00806644
happy: 0.13546932
sadness: 0.02680473
surprise: 0.01075646

The correct answer is obviously happy. This is a mistake that is less understandable but still the model is quite sure (~75%). There definitely is no hint of disgust in her face. Do note however, that happiness would be the classifier’s second guess. More training material might rectify this situation.

anger: 0.0372873
contempt: 0.08705531
disgust: 0.12282577
fear: 0.16857784
neutral: 0.09523397
happy: 0.26552763
sadness: 0.20521671
surprise: 0.01827547

The correct answer is sadness. Here the classifier is not sure at all (~27%)! Like in the previous image, the second guess (~20%) is the correct answer. This may very well be fixed by having more (and more diverse) training data.

anger: 0.01440529
contempt: 0.15626157
disgust: 0.01007962
fear: 0.00466321
neutral: 0.378776
happy: 0.00554828
sadness: 0.07485257
surprise: 0.35541345

The correct answer is surprise. Again a near miss (~38% vs ~36%)! Also note that this is particularly difficult because there are few baby faces in the dataset. When I said earlier that the extra google images are very challenging for a classifier, I meant it!

Upping the game – the ultimate challenge
Although the small google dataset I put together is more challenging than the lab-conditions of the CK/CK+ dataset, it is still somewhat controlled. For example I filtered out faces that were more sideways than frontal-facing, where the emotion was very mixed (happily surprised for example), and also where the emotion was so subtle that even I had trouble identifying it.

A far greater (and more realistic still) challenge is the SFEW/AFEW dataset, put together from a large collection of movie scenes. Read more about it here. The set is not publicly available but the author was generous enough to share the set with me so that I could evaluate the taken approach further.

Guess what, it fails miserably! It attained about 44.2% on the images when training on 90% and validating on 10% of the set. Although this is on par with what is mentioned in the paper, it shows there is still a long way to go before computers can recognize emotions with a high enough accuracy in real-life settings. There are also video clips included on which we will spend another post together with convolutional neural nets at a later time.

This set is particularly difficult because it contains different expressions and facial poses and rotations for similar emotions. This was the purpose of the authors; techniques by now are good enough to recognise emotions on controlled datasets with images taken in lab-like conditions, approaching upper 90% accuracy in many recent works (even our relatively simple approach reached early 90). However these sets do not represent real life settings very much, except maybe when using laptop webcams, because you always more or less face this device and sit at a comparable distance when using the laptop. This means for applications in marketing and similar fields the technology is already usable, albeit with much room for improvement still available and requiring some expertise to implement it correctly.

Final reflections
Before concluding I want you to take a moment, relax and sit back and think. Take for example the SFEW set with real-life examples, accurate classification of which quickly gets terribly difficult. We humans perform this recognition task remarkably well thanks to our highly complex visual system, which has zero problems with object rotation in all planes, different face sizes, different facial characteristics, extreme changes in lighting conditions or even partial occlusion of a face. Your first response might be “but that’s easy, I do it all the time!”, but it’s really, really, really not. Think for a moment about what an enormously complex problem this really is. I can show you a mouth and you would already be quite good at seeing an emotion. I can show you about 5% of a car and you could recognize it as a car easily, I can even warp and destroy the image and your brain would laugh at me and tell me “easy, that’s a car bro”. This is a task that you solve constantly and in real-time, without conscious effort, with virtually 100% accuracy and while only using the equivalent of ~20 watts for yourentire brain (not just the visual system). The average still-not-so-good-at-object-recognition CPU+GPU home computer uses 350-450 watts when computing. Then there’s supercomputers like theTaihuLight, which require about 15.300.000 watts (using in one hour what the average Dutch household uses in 5.1 years). At least at visual tasks, you still outperform these things by quite a large margin with only 0.00013% of their energy budget. Well done, brain!

Anyway, to try and tackle this problem digitally we need another approach. In another post we will look at various forms of neural nets (modeled after your brain) and how these may or may not solve the problem, and also at some other feature extraction techniques.

The CK+ dataset was used for validating and training of the classifier in this article, references to the set are:

• Kanade, T., Cohn, J. F., & Tian, Y. (2000). Comprehensive database for facial expression analysis. Proceedings of the Fourth IEEE International Conference on Automatic Face and Gesture Recognition (FG’00), Grenoble, France, 46-53.
• Lucey, P., Cohn, J. F., Kanade, T., Saragih, J., Ambadar, Z., & Matthews, I. (2010). The Extended Cohn-Kanade Dataset (CK+): A complete expression dataset for action unit and emotion-specified expression. Proceedings of the Third International Workshop on CVPR for Human Communicative Behavior Analysis (CVPR4HB 2010), San Francisco, USA, 94-101.

The SFEW/AFEW dataset used for evaluation is authored by and described in:

• A. Dhall, R. Goecke, S. Lucey and T. Gedeon, “Collecting Large, Richly Annotated Facial- Expression Databases from Movies”, IEEE MultiMedia 19 (2012) 34-41.
• A. Dhall, R. Goecke, J. Joshi, K. Sikka and T. Gedeon, “ Emotion Recognition In The Wild Challenge 2014: Baseline, Data and Protocol”, ACM ICMI 2014.

This the third part in a four part series about how to use Python for heart rate analysis. In this part you will learn about how to improve peak detection using a dynamic threshold, signal filtering, and outlier detection.

Important: The code in this tutorial is licensed under the GNU 3.0 open source license and you are free to modify and redistribute the code, given that you give others you share the code with the same right, and cite my name (use citation format below). You are not free to redistribute or modify the tutorial itself in any way. By reading on you agree to these terms. If you disagree, please navigate away from this page.

Citation format
van Gent, P. (2016). Analyzing a Discrete Heart Rate Signal Using Python. A tech blog about fun things with Python and embedded electronics. Retrieved from: http://www.paulvangent.com/2016/03/30/analyzing-a-discrete-heart-rate-signal-using-python-part-3

IE users: I’ve gotten several reports that sometimes the code blocks don’t display correctly or at all on Internet Explorer. Please refresh the page and they should display fine.

• Part 1: Opening the data, detecting the first peaks and calculating the BPM;
• Part 2: Extracting complex measures from the heart rate signal;
• Part 3: Signal filtering, improving detection with a dynamic threshold;
• Part 4: Detecting and rejecting noisy signal parts.

The following tutorial assumes basic knowledge of the Python programming language, FIR-filters andfast fourier transform methods. The tutorial should be suitable for those with most levels of Python skill. It is divided into separate parts so that you can easily skip over those parts you understand anyway. The data was collected with a simple data logger based on an ATTiny45, as described in another post.

All data in this tutorial is recorded with the Pulse Sensor, which I’ve found to be a great little device to easily log heart rate data from the fingertips. I have used it / am using it for research with video games, driving simulators and on-road experiments with car drivers.

Some Theory and Background
So far we’ve been over how to analyze a heart rate signal and extract the most widely used time domain and frequency domain measures from it. However, the signal we used was ideal. Now consider this signal:

A challenge! This is the other extreme of the signal qualities you will run into. To be honest I cheated and generated this signal by measuring while I was attaching the sensor to my finger (between 0 and about 4000). Immediately after this the blood vessels in the finger need to adapt to being compressed by the sensor (at about 4000-5000), after which the signal becomes stable. At about 7500, 9000 and 12000 I forcefully moved the sensor over my finger to create extra noise. I was surprised at how well the sensor suppresses noise by itself, so well done guys at pulsesensor.com.

Although this signal is manually ‘destroyed’ at points, in practice it can will happen that parts of your data contain noise or artefacts. The sensor might move which creates noise, it might not be correctly attached or become detached during the measurement, the participant might sneeze, move, or any other noise inducing factor might interfere.

We will see how to deal with this in a few stages:

• Getting started, evaluate the result of passing this signal to our algorithm from part two;
• Careful: Sampling Frequency;
• Filter the signal to remove unwanted frequencies (noise);
• Improving detection accuracy with a dynamic threshold;
• Detecting incorrectly detected / missed peaks;
• Removing errors and reconstructing the R-R signal to be error-free.

Getting Started
First let’s see how our algorithm from part 2 handles this signal. Download the dataset here.

1
2
3
  (...) line 37, in detect_peaks
    maximum = max(window)
ValueError: max() arg is an empty sequence

Great, it doesn’t. What is happening?

Some may have noticed it before, our detect_peaks() function will break in one eventuality: when the heart rate signal goes from being smaller than the moving average, to becoming equal to it without moving above it for at least two data points in a row. The most likely case of this happening is if the signal drops to zero for a period of time. The function will then skip to the else statement and try to detect peaks in the ROI, but no ROI has been marked because the signal never rose above the moving average, so max(window) throws the error because an empty list has no maximum. This happensbecause I’m an idiot because in the loop the situation where both signals are of equal height is not accounted for, which happens to happen early on in the signal, where I’m attaching the sensor to my finger:

If you didn’t notice this in the code before, don’t worry, neither did I at first. Now to update thedetect_peaks() function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def detect_peaks(dataset):
    window = []
    peaklist = []
    listpos = 0 
    for datapoint in dataset.hart:
        rollingmean = dataset.hart_rollingmean[listpos]
        if (datapoint <= rollingmean) and (len(window) <= 1): #Here is the update in (datapoint <= rollingmean)
            listpos += 1
        elif (datapoint > rollingmean):
            window.append(datapoint)
            listpos += 1
        else:
            maximum = max(window)
            beatposition = listpos - len(window) + (window.index(max(window)))
            peaklist.append(beatposition)
            window = []
            listpos += 1
    measures['peaklist'] = peaklist
    measures['ybeat'] = [dataset.hart[x] for x in peaklist]

For now also comment out line 19, we will get to this later on.

1
2
3
4
5
6
def rolmean(dataset, hrw, fs):
    mov_avg = pd.rolling_mean(dataset.hart, window=(hrw*fs))
    avg_hr = (np.mean(dataset.hart)) 
    mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
    #mov_avg = [x*1.2 for x in mov_avg]
    dataset['hart_rollingmean'] = mov_avg

And plot again, this time with peak detection:

The module is not throwing errors anymore and is detecting peaks, but it is hardly the result we’re looking for. Let’s start by looking at noise reduction to see if we can improve the signal.

Sampling Frequency
Before getting to signal filtering let’s determine the sample rate of our file. The file from the previous two parts was 100Hz, what about this one? In practice the actual recording frequency may vary with different devices or systems. It can also be (slightly) different than what the devices are rated for. The accuracy of all the calculated measures is dependent on accurately knowing the sample rate, so this is important to check. Even a difference of 100Hz and 101Hz can lead to significantly different results if you combine data from various sources. I usually log either ‘world time’ or ‘time since start of recording’ with every data line for the purpose of calculating and verifying sample rate afterwards.

With this it is straightforward to calculate the sampling rate:

1
2
3
4
5
6
7
8
9
10
#Simple way to get sample rate
sampletimer = [x for x in dataset.timer] #dataset.timer is a ms counter with start of recording at '0'
measures['fs'] = ((len(sampletimer) / sampletimer[-1])*1000) #Divide total length of dataset by last timer entry. This is in ms, so multiply by 1000 to get Hz value
#If your timer is a date time string, convert to UNIX timestamp to more easily calculate with, use something like this:
unix_time = []
for x in dataset.datetime:
    dt = datetime.datetime.strptime(Datum, "%Y-%m-%d %H:%M:%S.%f")
    unix_time.append(time.mktime(dt.timetuple()) + (dt.microsecond / 1000000.0))
measures['fs'] = (len(unix_time) / (unix_time[-1] - unix_time[0]))

The sample rate of the file provided here is actually 117Hz! The measures can now be calculated automatically using the determined sample rate.

Note that this is not the whole story, apart from sample rate you should also check your data for sampling spread; are all samples spaced evenly apart, e.g. are there no ‘gaps’ or ‘skips’ in the datastream? If your data contains gaps and skips, provided they are only a few samples long at maximum, consider interpolating the missing values before processing. If your sampling rate varies much over time you’re in trouble, use a different datalogging device.

Now that we know the sampling frequency more accurately, we can move on to signal filtering.

Signal Filtering
The first thing you should do when encountering artefacts or noisy signals (some would argue: always) is try to filter your signal. Why? Because in any real-life measuring situation your signal, whatever it may be, will consist of a signal part and an error part. This is because the perfect sensor does not exist, so it will always pick up interference from sources other than the one you are trying to measure. An example of common interference is power line noise. The AC power from the wall contacts in most countries has a frequency of 50Hz (some, such as the U.S., use 60Hz). This noise is present in many raw ECG-measurements as well.

An often used filter to reduce noise is the Butterworth Filter, which is characterized by a very even response to frequencies within the specified range. It acts as a ‘frequency gate’; suppressing frequencies beyond the specified cutoff range, more so as the frequencies move further away from it. This cutoff point is not a hard line. What I mean is that if you set the cutoff frequency at for example 5Hz, a signal of 5.1Hz will still be let through the filter, it will just be slightly reduced in amplitude (or ‘volume’, if that makes more sense). A signal of 10Hz on the other hand will only get through very weakly, if at all. This is also dependent on the filter order, as is illustrated nicely here.
Still difficult to understand? Think about it like this: you’re at a party and two people are talking to you simultaneously, leading to a situation where you cannot understand either of them. Now place a filter between you and the two others. The filter will reduce the speaking volume of person 1 without altering the volume of person 2. Now you can understand person 2 just fine. This is what the filter does except with frequencies.

Anyway, let’s get to coding and see if the signal can benefit from filtering. First download and open the dataset if you have not done it yet, define the filter using scipy.signal, filter and finally plot the signal. Assuming you’re working with the code from the previous part, define the filter and plot as such:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
from scipy.signal import butter, lfilter #Import the extra module required
#Define the filter
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs #Nyquist frequeny is half the sampling frequency
    normal_cutoff = cutoff / nyq 
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a
    
def butter_lowpass_filter(data, cutoff, fs, order):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y
dataset = get_data('data2.csv')
dataset = dataset[6000:12000].reset_index(drop=True) #For visibility take a subselection of the entire signal from samples 6000 - 12000 (00:01:00 - 00:02:00)
filtered = butter_lowpass_filter(dataset.hart, 2.5, 100.0, 5)#filter the signal with a cutoff at 2.5Hz and a 5th order Butterworth filter
#Plot it
plt.subplot(211)
plt.plot(dataset.hart, color='Blue', alpha=0.5, label='Original Signal')
plt.legend(loc=4)
plt.subplot(212)
plt.plot(filtered, color='Red', label='Filtered Signal')
plt.ylim(200,800) #limit filtered signal to have same y-axis as original (filter response starts at 0 so otherwise the plot will be scaled)
plt.legend(loc=4)
plt.show()

Now there doesn’t seem to be much improvement in this signal. If you look closely the lines are a little smoother, but there wasn’t a lot of high-amplitude, high-frequency noise there to begin with. It could even be argued that filtering slightly worsens the parts with the lower frequency noise, because there it suppresses the R-peaks a little as well. This is a good example of why you should always plot your signal when you decide to filter it. Filtering the signal changes it, and it is up to you to decide whether this change is for the better.

An example of where a Butterworth Filter was very useful, is this noisy signal I worked with in another project:

No question this improved the signal enough to process it further.

Improving Detection Accuracy With a Dynamic Threshold
The first and maybe most obvious way to reduce the incorrect labeling of the secondary peaks is to raise the moving average we use as a threshold. But raise to what level? This will be different for many different signals. We need measures to help determine which threshold level is probably the most accurate.

A few simple measures can help, we can:

• Look at signal length and count how many peaks there are vs how many we would expect;
• Determine and use the standard deviation of RR intervals (let’s call it RRSD).

The amount of detected peaks holds valuable information but only works to reject obvious incorrect thresholds. Depending on your application (most of mine are with people sitting still), we can reject unlikely bpm values. For example I reject thresholds resulting in average bpm of <30bpm and >130bpm. In other situations (physical excercise) the threshold can be different.

The RRSD tells us something about how spread out the differences between all RR-intervals are. Generally if there is not too much noise, the detection that has both the lowest RRSD that is not zeroand also a believable BPM value will be the best fit. RRSD must not be zero because that means there is no heart rate variability, which is a strong indication of mistakes in the detection of R-peaks.

This simple approach works because missing a beat will lead to an RR interval that is about twice as big as the average RR interval, while incorrectly labeling a beat will lead to an RR interval that is at mostabout half as big as the average RR interval, but generally smaller. Both situations lead to anincreased RRSD. In essence we exploit the fact that the heart rate signal contains a fairly stable, recurring pattern.

To illustrate let’s plot four peak detection rounds in a subselection of the dataset, with the moving average raised by 0%, 10%, 25% and 35% (top to bottom):

In the second-to-last plot all R-peaks are detected correctly and nothing has been marked as an R-peak incorrectly. Note that, although the BPM on its own could be valid for all four plots (none of them is an absolute mess), the RRSD strongly points to the plot which is most correct. I say ‘most correct’ because there are situations where some errors will remain no matter how you position the threshold, more on this later. Also note how the missing of a single peak in the last plot already causes the RRSD to increase quite a bit compared to the third one.

Now how to implement this? We cannot simply run 10.000 variations, each with a slightly more raised moving average. Apart from that this would cost us severely in overall performance of the algorithm, it would also be very redundant because many iterations would yield the same correct result (and many others the same incorrect result!). A possible solution is to check with intervals, and afterwards evaluate the most likely RRSD and BPM pair, like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def detect_peaks(dataset, ma_perc, fs): #Change the function to accept a moving average percentage 'ma_perc' argument
    rolmean = [(x+((x/100)*ma_perc)) for x in dataset.hart_rollingmean] #Raise moving average with passed ma_perc
    window = []
    peaklist = []
    listpos = 0 
    for datapoint in dataset.hart:
        rollingmean = rolmean[listpos]
        if (datapoint <= rollingmean) and (len(window) <= 1): #Here is the update in (datapoint <= rollingmean)
            listpos += 1
        elif (datapoint > rollingmean):
            window.append(datapoint)
            listpos += 1
        else:
            maximum = max(window)
            beatposition = listpos - len(window) + (window.index(max(window)))
            peaklist.append(beatposition)
            window = []
            listpos += 1
    measures['peaklist'] = peaklist
    measures['ybeat'] = [dataset.hart[x] for x in peaklist]
    measures['rolmean'] = rolmean
    calc_RR(dataset, fs)
    measures['rrsd'] = np.std(measures['RR_list'])
def fit_peaks(dataset, fs):
    ma_perc_list = [5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 150, 200, 300] #List with moving average raise percentages, make as detailed as you like but keep an eye on speed
    rrsd = []
    valid_ma = []
    for x in ma_perc_list: #Detect peaks with all percentages, append results to list 'rrsd'
        detect_peaks(dataset, x, fs)
        bpm = ((len(measures['peaklist'])/(len(dataset.hart)/fs))*60)
        rrsd.append([measures['rrsd'], bpm, x])
        
    for x,y,z in rrsd: #Test list entries and select valid measures
        if ((x > 1) and ((y > 30) and (y < 130))):
            valid_ma.append([x, z])
            
    measures['best'] = min(valid_ma, key = lambda t: t[0])[1] #Save the ma_perc for plotting purposes later on (not needed)
    detect_peaks(dataset, min(valid_ma, key = lambda t: t[0])[1], fs) #Detect peaks with 'ma_perc' that goes with lowest rrsd

Now run and plot on the dataset (timed the entire algorithm so far including loading and preprocessing at about 151ms, single core performance on an i7-4790, so it’s still still pretty quick. Multithreading will speed this up a lot):

That is already a lot better. It’s not dropping any correct R-peaks, but there are still a few incorrect detections remaining, and there is also the part from 0 to about 5000 which contains little to no heart rate signal. We will come back to this noisy segment and see how to detect and exclude segments of noise in part 4.

For now let’s take out the noisy part at the beginning and see how we can detect and reject outliers.

Finding Incorrectly Detected Peaks
The last thing to look at is how to detect and reject abnormal peak positions. One way to do this is to make use of the fact that the heart rate changes gradually rather than abruptly. Your bpm for example will not skip from 60bpm to 120bpm in a single beat or vice versa, so let’s make use of that.

Again this means returning to RR-intervals. Remember that if a peak is skipped in the detection, or two peaks are marked in stead of one, the resulting RR-interval will be a lot larger or smaller than the average interval. We can set a threshold and mark intervals that exceed it, similar to how we detected the peaks. For the data I collect this is enough.

There is, however, one potential complication. If you analyze a very long signal at once, wherein the heart rate changes a lot over time, this can lead to incorrect rejections. Imagine a signal with a gradually increasing heart rate, starting from 60 bpm and ending at 180bpm. This means a steady trend of decreasing RR-intervals, which is indicative of the speeding up of the heart rate rather than mistakes in the detection of R-peaks. By using just thresholds based on the mean RR-interval, this will result in a rejection of the first and last portion of the signal. If this happens in your data you could detrend RR_list first. Using scipy.signal, this is easy:

1
2
3
4
from scipy import signal
RR_list = measures['RR_list'] #First retrieve list from dictionary
RR_list_detrended = signal.detrend(RR_list, type='linear')

However, if your signal contains a period of large increases followed by similarly large decreases in heart rate, you will need to employ other methods. The solution is beyond the scope of this tutorial series, but if you have this problem you may want to use a high pass filter with a very low cutoff frequency. Another way could be to split the signal in to smaller portions (so that the increase and decrease trends are separated), detrend linearly and average the measures. If the smaller portions are not of equal length, be sure to weight the measures before averaging.
Naturally, do not calculate any measures with the detrended RR-signal, only use it for the detection of errors in peak marking.

Back to outlier rejection. For the thresholds, in practice I’ve found a threshold level of the mean of RR-differences with a band of 250-300ms on both ends works well. Using the previous code and limiting the dataset to [5000:15000] (to exclude the noisy beginning for now), implement like so:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
RR_list = measures['RR_list'] #Get measures
peaklist = measures['peaklist']
ybeat = measures['ybeat']
upper_threshold = (np.mean(RR_list) + 300) #Set thresholds
lower_threshold = (np.mean(RR_list) - 300)
#detect outliers
cnt = 0
removed_beats = []
removed_beats_y = []
RR2 = []
while cnt < len(RR_list):
    if (RR_list[cnt] < upper_threshold) and (RR_list[cnt] > lower_threshold):
        RR2.append(RR_list[cnt])
        cnt += 1
    else:   
        removed_beats.append(peaklist[cnt])
        removed_beats_y.append(ybeat[cnt])
        cnt += 1
measures['RR_list_cor'] = RR2 #Append corrected RR-list to dictionary
plt.subplot(211)
plt.title('Marked Uncertain Peaks')
plt.plot(dataset.hart, color='blue', alpha=0.6, label='heart rate signal')
plt.plot(measures['rolmean'], color='green')
plt.scatter(measures['peaklist'], measures['ybeat'], color='green')
plt.scatter(removed_beats, removed_beats_y, color='red', label='Detection uncertain')
plt.legend(framealpha=0.6, loc=4)
plt.subplot(212)
plt.title("RR-intervals with thresholds")
plt.plot(RR_list)
plt.axhline(y=upper_threshold, color='red')
plt.axhline(y=lower_threshold, color='red')
plt.show()

Resulting in:

It seems we got all the little buggers. The result is a plot with all correctly detected R-peaks marked green. The incorrect ones are marked red. The generated list measures[‘RR_list_cor’] has the RR-list without those belonging to incorrect peaks in it.

In part 4 we will look into how to mark and exclude noise segments and a few other optimizations.

Rounding up
Now tidy up all code, and also to update some functions to accept new variables and insert the peak rejection function.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.interpolate import interp1d
from scipy.signal import butter, lfilter, detrend
measures = {}
#Preprocessing
def get_data(filename):
    dataset = pd.read_csv(filename)
    return dataset
def get_samplerate(dataset):
    sampletimer = [x for x in dataset.datetime]
    measures['fs'] = ((len(sampletimer) / sampletimer[-1])*1000)
def rolmean(dataset, hrw, fs):
    mov_avg = pd.rolling_mean(dataset.hart, window=(hrw*fs), center=False)
    avg_hr = (np.mean(dataset.hart)) 
    mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
    dataset['hart_rollingmean'] = mov_avg
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a
def butter_lowpass_filter(data, cutoff, fs, order):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y    
def filtersignal(data, cutoff, fs, order):
    hart = [math.pow(x, 3) for x in dataset.hart]
    hartfiltered = butter_lowpass_filter(hart, cutoff, fs, order)
    dataset['hart'] = hartfiltered
#Peak detection
def detect_peaks(dataset, ma_perc, fs):
    rolmean = [(x+((x/100)*ma_perc)) for x in dataset.hart_rollingmean]
    window = []
    peaklist = []
    listpos = 0 
    for datapoint in dataset.hart:
        rollingmean = rolmean[listpos]
        if (datapoint <= rollingmean) and (len(window) <= 1): #Here is the update in (datapoint <= rollingmean)
            listpos += 1
        elif (datapoint > rollingmean):
            window.append(datapoint)
            listpos += 1
        else:
            maximum = max(window)
            beatposition = listpos - len(window) + (window.index(max(window)))
            peaklist.append(beatposition)
            window = []
            listpos += 1
    measures['peaklist'] = peaklist
    measures['ybeat'] = [dataset.hart[x] for x in peaklist]
    measures['rolmean'] = rolmean
    calc_RR(dataset, fs)
    measures['rrsd'] = np.std(measures['RR_list'])
def fit_peaks(dataset, fs):
    ma_perc_list = [5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 150, 200, 300] #List with moving average raise percentages, make as detailed as you like but keep an eye on speed
    rrsd = []
    valid_ma = []
    for x in ma_perc_list:
        detect_peaks(dataset, x, fs)
        bpm = ((len(measures['peaklist'])/(len(dataset.hart)/fs))*60)
        rrsd.append([measures['rrsd'], bpm, x])
    for x,y,z in rrsd:
        print x, y, z
        if ((x > 1) and ((y > 30) and (y < 130))):
            valid_ma.append([x, z]) 
    measures['best'] = min(valid_ma, key = lambda t: t[0])[1]
    detect_peaks(dataset, min(valid_ma, key = lambda t: t[0])[1], fs)
def check_peaks(dataset):
    RR_list = measures['RR_list']
    peaklist = measures['peaklist']
    ybeat = measures['ybeat']
    upper_threshold = np.mean(RR_list) + 300
    lower_threshold = np.mean(RR_list) - 300
    removed_beats = []
    removed_beats_y = []
    RR_list_cor = []
    peaklist_cor = [peaklist[0]]
    cnt = 0
    while cnt < len(RR_list):
        if (RR_list[cnt] < upper_threshold) and (RR_list[cnt] > lower_threshold):
            RR_list_cor.append(RR_list[cnt])
            peaklist_cor.append(peaklist[cnt+1]) 
            cnt += 1
        else:    
            removed_beats.append(peaklist[cnt+1])
            removed_beats_y.append(ybeat[cnt+1])
            cnt += 1
    measures['RR_list_cor'] = RR_list_cor
    measures['peaklist_cor'] = peaklist_cor
#Calculating all measures
def calc_RR(dataset, fs):
    peaklist = measures['peaklist']
    RR_list = []
    cnt = 0
    while (cnt < (len(peaklist)-1)):
        RR_interval = (peaklist[cnt+1] - peaklist[cnt])
        ms_dist = ((RR_interval / fs) * 1000.0)
        RR_list.append(ms_dist)
        cnt += 1
    RR_diff = []
    RR_sqdiff = []
    cnt = 0 
    while (cnt < (len(RR_list)-1)): 
        RR_diff.append(abs(RR_list[cnt] - RR_list[cnt+1])) 
        RR_sqdiff.append(math.pow(RR_list[cnt] - RR_list[cnt+1], 2))
        cnt += 1
    measures['RR_list'] = RR_list
    measures['RR_diff'] = RR_diff
    measures['RR_sqdiff'] = RR_sqdiff
    
def calc_ts_measures(dataset):
    RR_list = measures['RR_list_cor']
    RR_diff = measures['RR_diff']
    RR_sqdiff = measures['RR_sqdiff']
    measures['bpm'] = 60000 / np.mean(RR_list)
    measures['ibi'] = np.mean(RR_list)
    measures['sdnn'] = np.std(RR_list)
    measures['sdsd'] = np.std(RR_diff)
    measures['rmssd'] = np.sqrt(np.mean(RR_sqdiff))
    NN20 = [x for x in RR_diff if (x>20)]
    NN50 = [x for x in RR_diff if (x>50)]
    measures['nn20'] = NN20
    measures['nn50'] = NN50
    measures['pnn20'] = float(len(NN20)) / float(len(RR_diff))
    measures['pnn50'] = float(len(NN50)) / float(len(RR_diff))
def calc_fd_measures(dataset, fs):
    peaklist = measures['peaklist_cor']
    RR_list = measures['RR_list_cor']
    RR_x = peaklist[1:]
    RR_y = RR_list
    RR_x_new = np.linspace(RR_x[0],RR_x[-1],RR_x[-1])
    f = interp1d(RR_x, RR_y, kind='cubic')
    n = len(dataset.hart)
    frq = np.fft.fftfreq(len(dataset.hart), d=((1/fs)))
    frq = frq[range(n/2)]
    Y = np.fft.fft(f(RR_x_new))/n
    Y = Y[range(n/2)]
    measures['lf'] = np.trapz(abs(Y[(frq>=0.04) & (frq<=0.15)]))
    measures['hf'] = np.trapz(abs(Y[(frq>=0.16) & (frq<=0.5)]))
#Plotting it
def plotter(dataset):
    peaklist = measures['peaklist']
    ybeat = measures['ybeat']
    plt.title("Best fit: mov_avg %s percent raised" %measures['best'])
    plt.plot(dataset.hart, alpha=0.5, color='blue', label="heart rate signal")
    plt.plot(measures['rolmean'], color ='green', label="moving average")
    plt.scatter(peaklist, ybeat, color='red', label="RRSD:%.2f\nBPM:%.2f" %(np.std(measures['RR_list']), measures['bpm']))#, label="average: %.1f BPM" %measures['bpm'])
    plt.legend(loc=4, framealpha=0.6) 
    plt.show() 
#Wrapper function
def process(dataset, hrw, fs):
    filtersignal(dataset, 2.5, fs, 5)
    rolmean(dataset, hrw, fs)
    fit_peaks(dataset, fs)
    calc_RR(dataset, fs)
    check_peaks(dataset)
    calc_ts_measures(dataset)
    calc_fd_measures(dataset, fs)