Added LP module and tests

2 files changed, 140 insertions, 0 deletions

diff --git a/lp.py b/lp.py
new file mode 100644
index 0000000..740ca97
--- /dev/null
+++ b/lp.py
@@ -0,0 +1,86 @@
+# Copyright (c) 2009 John Glover, National University of Ireland, Maynooth
+# 
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of 
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+import simpl
+import numpy as np
+
+def burg(signal, order):
+    """Using the burg method, calculate order coefficients.
+    Returns a numpy array."""
+    coefs = np.array([1.0])
+    # initialise f and b - the forward and backwards predictors
+    f = np.zeros(len(signal))
+    b = np.zeros(len(signal))
+    for i in range(len(signal)):
+        f[i] = b[i] = signal[i]
+    # burg algorithm
+    for k in range(order):
+        # fk is f without the first element
+        fk = f[1:]
+        # bk is b without the last element
+        bk = b[0:b.size-1]
+        # calculate mu
+        if sum((fk*fk)+(bk*bk)):
+            # check for division by zero
+            mu = -2.0 * sum(fk*bk) / sum((fk*fk)+(bk*bk))
+        else:
+            mu = 0.0
+        # update coefs
+        # coefs[::-1] just reverses the coefs array
+        coefs = np.hstack((coefs,0)) + (mu * np.hstack((0, coefs[::-1])))
+        # update f and b
+        f = fk + (mu*bk)
+        b = bk + (mu*fk)
+    return coefs[1:]
+
+def predict(signal, coefs, num_predictions):
+    """Using Linear Prediction, return the estimated next num_predictions
+    values of signal, using the given coefficients.
+    Returns a numpy array."""
+    predictions = np.zeros(num_predictions)
+    past_samples = np.zeros(len(coefs))
+    for i in range(len(coefs)):
+        past_samples[i] = signal[-1-i]
+    sample_pos = 0
+    for i in range(num_predictions):
+        # each sample in past_samples is multiplied by a coefficient
+        # results are summed
+        for j in range(len(coefs)):
+            predictions[i] -= coefs[j] * past_samples[(j+sample_pos) % len(coefs)]
+        sample_pos -= 1
+        if sample_pos < 0:
+            sample_pos = len(coefs) - 1
+        past_samples[sample_pos] = predictions[i] 
+    return predictions
+
+
+class LPPartialTracking(simpl.PartialTracking):
+    "Partial tracking, based on the McAulay and Quatieri (MQ) algorithm"    
+    def __init__(self):
+        simpl.PartialTracking.__init__(self)
+        self._matching_interval = 100  # peak matching interval, in Hz
+
+        
+    def update_partials(self, frame, frame_number):
+        """Streamable (real-time) LP peak-tracking.
+        """
+        frame_partials = []
+                
+        for peak in frame:
+            pass
+
+        return frame_partials
+
diff --git a/tests/lp.py b/tests/lp.py
new file mode 100644
index 0000000..c6051dd
--- /dev/null
+++ b/tests/lp.py
@@ -0,0 +1,54 @@
+# Copyright (c) 2010 John Glover, National University of Ireland, Maynooth
+# 
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of 
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+import simpl
+from simpl import lp
+import numpy as np
+import matplotlib.pyplot as plt
+import unittest
+
+audio, sampling_rate = simpl.read_wav("audio/flute.wav")
+
+# take just the first few frames
+audio = audio_in[0:4096]
+# Peak detection and partial tracking using SMS
+pd = simpl.SndObjPeakDetection()
+pd.max_peaks = 60
+peaks = pd.find_peaks(audio)
+pt = simpl.MQPartialTracking()
+pt.max_partials = 60
+partials = pt.find_partials(peaks)
+simpl.plot.plot_partials(partials)
+plt.show()
+
+exit()
+
+
+class TestLP(unittest.TestCase):
+    FLOAT_PRECISION = 3 # number of decimal places to check for accuracy
+
+    def test_predict(self):
+        """test_predict"""
+        coefs = np.array([1,2,3,4,5])
+        test_signal = np.ones(5)
+        predictions = lp.predict(test_signal, coefs, 2)
+        self.assertEquals(predictions[0], -sum(coefs))
+        self.assertEquals(predictions[1], -sum(coefs[1:])-predictions[0])
+    
+suite = unittest.TestSuite()
+suite.addTest(TestLP('test_predict'))
+unittest.TextTestRunner().run(suite)
+