Merge branch 'restructure' into threaded

* restructure:
  Changed make_filter back to just returning a function.
  Added tests for all the signals functions, and did some cleaning/formatting.
  Fixed bug in find_periodicities which resulted in incorrect frequencies
  Moved signal-processing and utility functions to separate files.

Conflicts:
	main.py

Author: Joseph Atkins-Turkish

facetracking.py

@@ -0,0 +1,21 @@
+
+def make_face_rects(rect):
+    """ Given a rectangle (covering a face), return two rectangles.
+        which cover the forehead and the area under the eyes """
+    x, y, w, h = rect
+
+    rect1_x = x + w / 4.0
+    rect1_w = w / 2.0
+    rect1_y = y + 0.05 * h
+    rect1_h = h * 0.9 * 0.2
+
+    rect2_x = rect1_x
+    rect2_w = rect1_w
+
+    rect2_y = y + 0.05 * h + (h * 0.9 * 0.55)
+    rect2_h = h * 0.9 * 0.45
+
+    return (
+        (int(rect1_x), int(rect1_y), int(rect1_w), int(rect1_h)),
+        (int(rect2_x), int(rect2_y), int(rect2_w), int(rect2_h))
+    )

main.py

@@ -7,7 +7,7 @@
 import itertools
 from scipy import interpolate, signal
 import sklearn.decomposition
-
+from signals import *
 import threading
 import Queue
 
@@ -46,7 +46,7 @@ def __init__(self,frameQueue,incrementalPCA,n_components,windowSize):
         self.daemon = True
         self.frameQueue = frameQueue
         self.windowSize = windowSize
-        if incremental:
+        if incrementalPCA:
             self.pca = sklearn.decomposition.IncrementalPCA(n_components=n_components)
         else:
             self.pca = sklearn.decomposition.PCA(n_components=n_components)
@@ -192,135 +192,6 @@ def get_moving_points(source_queue, do_draw=True, n_points=100):
         print "Heart rate by peak estimate: {} BPM".format(countbpm)
 ########################################################################################
 
-def make_face_rects(rect):
-    """ Given a rectangle (covering a face), return two rectangles.
-        which cover the forehead and the area under the eyes """
-    x, y, w, h = rect
-
-    rect1_x = x + w / 4.0
-    rect1_w = w / 2.0
-    rect1_y = y + 0.05 * h
-    rect1_h = h * 0.9 * 0.2
-
-    rect2_x = rect1_x
-    rect2_w = rect1_w
-
-    rect2_y = y + 0.05 * h + (h * 0.9 * 0.55)
-    rect2_h = h * 0.9 * 0.45
-
-    return (
-        (int(rect1_x), int(rect1_y), int(rect1_w), int(rect1_h)),
-        (int(rect2_x), int(rect2_y), int(rect2_w), int(rect2_h))
-    )
-
-
-
-
-def window(seq, n=2, skip=1):
-    "Returns a sliding window (of width n) over data from the iterable"
-    "   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   "
-    # http://stackoverflow.com/questions/6822725/rolling-or-sliding-window-iterator-in-python
-
-    it = iter(seq)
-    result = tuple(itertools.islice(it, n))
-    if len(result) == n:
-        yield result
-    i = 0
-    for elem in it:
-        result = result[1:] + (elem,)
-        i = (i + 1) % skip
-        if i == 0:
-            yield result
-
-
-def interpolate_points(points, ratio=30.0 / 250.0, axis=0):
-    """ Given an matrix of waveforms, interpolate them along an axis
-    such that the number of new is multiplied by (1/ratio) """
-    # Define the old time space, i.e. the index of each point
-    N = points.shape[axis]
-    indices = np.arange(0, N)
-    # Make an 'interpolation function' using scikit's interp1d
-    f = interpolate.interp1d(indices, points, kind='cubic', axis=axis)
-    # Define the new time axis,
-    xnew = np.arange(0, N - 1, ratio)
-    return f(xnew)
-
-
-def filter_unstable_movements(points):
-    """ Filter unstable movements, e.g. coughing """
-    """ In the paper, they removed points which had maximum movements
-    greater than the "mode" of rounded maximum movements. Or something.
-    This didn't really work or make sense when we tried it, so we tried
-    something else, then gave up... """
-    maximums = np.max(np.diff(points.T), axis=1)
-    median = np.median(maximums)
-    return points[:, maximums > median]
-
-
-def make_filter(order=5, low_freq=0.75, high_freq=5, sample_freq=250.0):
-    """ Make the butterworth filter function required by the pulse paper"""
-    nyq = 0.5 * sample_freq
-    low = low_freq / nyq
-    high = high_freq / nyq
-    b, a = signal.butter(order, [low, high], btype='bandpass')
-    return lambda x: signal.lfilter(b, a, x)
-
-
-def find_periodicities(X, sample_freq=250.0):
-    """ Find the periodicity of each signal in a matrix(along axis 0),
-    and the associated frequencies of the periods"""
-
-    # We're not sure if this is quite correct, but it's what the paper
-    # seemed to imply...
-    # This could also be made much neater, and optimised.
-
-    # Find the power spectrum of the signal (absolute fft squared)
-    power = np.abs(np.fft.rfft(X, axis=0))**2
-
-    # Build a list of the actual frequencies corresponding to each fft index, using numpy's rfftfreq
-    # n.b. This is where I'm having some trouble. I don't think I'm actually getting the right
-    # numbers out for the frequencies of these signals...
-
-    real_frequencies = np.fft.rfftfreq(
-        power.shape[0],  d=(1 / (sample_freq)))
-
-    # Find the most powerful non-zero frequency in each signal
-    max_indices = np.argmax(power[1:, :], axis=0) + 1
-
-    # The first haromic component of f = f*2
-    harmonic_indices = max_indices * 2
-
-    # Initialise arrays for return values
-    periodicities = []
-    frequencies = []
-    i = 0
-
-    # Loop over each signal
-    for i1, i2 in zip(max_indices, harmonic_indices):
-        # Get the real frequency highest power component
-        frequencies.append(real_frequencies[i1])
-
-        # Get the total power of the highest power component and its
-        # first harmonic, as a percentage of the total signal power
-        period_power = np.sum(power[[i1, i2], i])
-        total_power = np.sum(power[:, i])
-        percentage = period_power / total_power
-
-        # That number is (apparently) the signal periodicity
-        periodicities.append(percentage)
-        i += 1
-
-    return np.array(frequencies), np.array(periodicities)
-
-
-def getpeaks(x, winsize=151):
-    """ Return the indices of all points in a signal which are the largest poitns in
-    a window centered on themselves """
-    for i, win in enumerate(window(x, winsize, 1)):
-        ind = int(winsize / 2)
-        if np.argmax(win) == ind:
-            yield i + ind
-
 
 def main(incrementalPCA=False):
     """ Run the full algorithm on a video """

signals.py

@@ -0,0 +1,112 @@
+import numpy as np
+import itertools
+from scipy import interpolate, signal
+
+
+def window(seq, n=2, skip=1):
+    "Returns a sliding window (of width n) over data from the iterable"
+    "   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   "
+    # http://stackoverflow.com/questions/6822725/rolling-or-sliding-window-iterator-in-python
+
+    it = iter(seq)
+    result = tuple(itertools.islice(it, n))
+    if len(result) == n:
+        yield result
+    i = 0
+    for elem in it:
+        result = result[1:] + (elem,)
+        i = (i + 1) % skip
+        if i == 0:
+            yield result
+
+
+def interpolate_points(points, ratio=30.0 / 250.0, axis=0):
+    """ Given an matrix of waveforms, interpolate them along an axis
+    such that the number of new is multiplied by (1/ratio) """
+    # Define the old time space, i.e. the index of each point
+    N = points.shape[axis]
+    indices = np.arange(0, N)
+    # Make an 'interpolation function' using scikit's interp1d
+    f = interpolate.interp1d(indices, points, kind='cubic', axis=axis)
+    # Define the new time axis,
+    xnew = np.arange(0, N - 1, ratio)
+    return f(xnew)
+
+
+def filter_unstable_movements(points):
+    """ Filter unstable movements, e.g. coughing """
+    """ In the paper, they removed points which had maximum movements
+    greater than the "mode" of rounded maximum movements. Or something.
+    This didn't really work or make sense when we tried it, so we tried
+    something else, then gave up... """
+    maximums = np.max(np.diff(points.T), axis=1)
+    median = np.median(maximums)
+    return points[:, maximums > median]
+
+
+def make_filter(order=5, low_freq=0.75, high_freq=5, sample_freq=250.0):
+    """ Make the butterworth filter function required by the pulse paper"""
+    nyq = 0.5 * sample_freq
+    low = low_freq / nyq
+    high = high_freq / nyq
+    b, a = signal.butter(order, [low, high], btype='bandpass')
+    func = lambda x: signal.lfilter(b, a, x)
+    func.b = b
+    func.a = a
+    return func
+
+
+def find_periodicities(X, sample_freq=250.0):
+    """ Find the periodicity of each signal in a matrix(along axis 0),
+    and the associated frequencies of the periods"""
+
+    # We're not sure if this is quite correct, but it's what the paper
+    # seemed to imply...
+    # This could also be made much neater, and optimised.
+
+    # Find the power spectrum of the signal (absolute fft squared)
+    power = np.abs(np.fft.rfft(X, axis=0))**2
+
+    # Build a list of the actual frequencies corresponding to each fft index, using numpy's rfftfreq
+    # n.b. This is where I'm having some trouble. I don't think I'm actually getting the right
+    # numbers out for the frequencies of these signals...
+
+    real_frequencies = np.fft.rfftfreq(
+        X.shape[0],  d=(1 / (sample_freq)))
+
+    # Find the most powerful non-zero frequency in each signal
+    max_indices = np.argmax(power[1:, :], axis=0) + 1
+
+    # The first haromic component of f = f*2
+    harmonic_indices = max_indices * 2
+
+    # Initialise arrays for return values
+    periodicities = []
+    frequencies = []
+    i = 0
+
+    # Loop over each signal
+    for fundamental, harmonic in zip(max_indices, harmonic_indices):
+        # Get the real frequency highest power component
+        frequencies.append(real_frequencies[fundamental])
+
+        # Get the total power of the highest power component and its
+        # first harmonic, as a percentage of the total signal power
+        period_power = np.sum(power[[fundamental, harmonic], i])
+        total_power = np.sum(power[:, i])
+        percentage = period_power / total_power
+
+        # That number is (apparently) the signal periodicity
+        periodicities.append(percentage)
+        i += 1
+
+    return np.array(frequencies), np.array(periodicities)
+
+
+def getpeaks(x, winsize=151):
+    """ Return the indices of all points in a signal which are the largest poitns in
+    a window centered on themselves """
+    for i, win in enumerate(window(x, winsize, 1)):
+        ind = int(winsize / 2)
+        if np.argmax(win) == ind:
+            yield i + ind

tests.py

@@ -0,0 +1,80 @@
+import unittest
+import numpy as np
+import signals
+import matplotlib.pyplot as plt
+import math
+
+
+class TestSignalFunctions(unittest.TestCase):
+
+    def make_test_signal(self, sample_rate):
+        """ Make a 10s signal with three frequency components
+        f = [100, 10, 1]/(2pi) """
+        x = np.arange(0, 10, 1.0 / float(sample_rate))
+        y = np.sin(x * 10.0)
+        y += np.sin(x * 1.0) * 0.7
+        y += np.sin(x * 100.0) * 0.1
+        return (x, y)
+
+    def test_get_peaks(self):
+        """ Check that get_peaks finds approximately the right locations and
+        number of peaks in a sine wave """
+        x = np.arange(0, math.pi * 8, 0.01)
+        y = np.sin(x)
+        test_peaks = x[list(signals.getpeaks(y, 151))] / (math.pi / 2.0)
+        real_peaks = [1, 5, 9, 13]
+        np.testing.assert_allclose(test_peaks, real_peaks, rtol=1e-3)
+
+    def test_make_filter(self):
+        """ Check that the bandpass filter function filters the correct signals.
+        N.B graph must be checked by eye! """
+        filt = signals.make_filter(
+            order=5, low_freq=0.75, high_freq=5, sample_freq=250.0)
+        x, y = self.make_test_signal(250.0)
+
+        yf = filt(y)
+        f = plt.figure()
+        plt.plot(x, y)
+        plt.plot(x, yf)
+        plt.xlabel('Time (250 Hz sample rate)')
+        plt.ylabel('Signal')
+        plt.legend(['Original', 'Filtered'])
+        f.savefig('test_make_filter.png', bbox_inches='tight')
+
+    def test_window(self):
+        """ Test the moving window function """
+        windows = list(signals.window(xrange(10), 4, 2))
+        expected = [(0, 1, 2, 3), (2, 3, 4, 5), (4, 5, 6, 7), (6, 7, 8, 9)]
+        self.assertEquals(windows, expected)
+
+    def test_interpolate_points(self):
+        """ Test point interpolation
+        N.B. Graph must be checked by eye! """
+        x, y = self.make_test_signal(30.0)
+        y = np.vstack((y, -y))
+        y0 = signals.interpolate_points(y.T, 30.0 / 250.0, axis=0)
+        y1 = signals.interpolate_points(y, 30.0 / 250.0, axis=1)
+        self.assertTrue(y0.shape[0] >= 2490)
+        self.assertTrue(y1.shape[1] >= 2490)
+        f = plt.figure()
+        plt.subplot(2, 1, 1)
+        plt.plot(y.T)
+        plt.subplot(2, 1, 2)
+        plt.plot(y0)
+        f.savefig('test_interpolate_points.png', bbox_inches='tight')
+
+    def test_find_periodicities(self):
+        """ Test that find_periodicities picks the correct signal as the
+        most periodic one, and returns a sensible frequency """
+        sample_freq = 30.0
+        x, y_periodic = self.make_test_signal(sample_freq)
+        y_flat = x / 8.0
+        y = np.array([y_periodic, y_flat])
+        frequencies, periodicies = signals.find_periodicities(
+            y.T, sample_freq=sample_freq)
+        self.assertGreater(periodicies[0], periodicies[1])
+        self.assertAlmostEqual(frequencies[0], 10 / (2 * math.pi), places=1)
+
+
+if __name__ == '__main__':
+    unittest.main()