[base] Add C++ implementation of Synthesis and Residual.

4 files changed, 368 insertions, 27 deletions

diff --git a/simpl/base.pyx b/simpl/base.pyx
index bb34e14..b126be1 100644
--- a/simpl/base.pyx
+++ b/simpl/base.pyx
@@ -92,6 +92,38 @@ cdef extern from "../src/simpl/base.h" namespace "simpl":
         vector[c_Peak*] update_partials(c_Frame* frame)
         vector[c_Frame*] find_partials(vector[c_Frame*] frames)
 
+    cdef cppclass c_Synthesis "simpl::Synthesis":
+        c_Synthesis()
+        int frame_size()
+        void frame_size(int new_frame_size)
+        int next_frame_size()
+        int hop_size()
+        void hop_size(int new_hop_size)
+        int sampling_rate()
+        void sampling_rate(int new_sampling_rate)
+        int max_partials()
+        void max_partials(int new_max_partials)
+        void synth_frame(c_Frame* frame)
+        vector[c_Frame*] synth(vector[c_Frame*] frames)
+
+    cdef cppclass c_Residual "simpl::Residual":
+        c_Synthesis()
+        int frame_size()
+        void frame_size(int new_frame_size)
+        int next_frame_size()
+        int hop_size()
+        void hop_size(int new_hop_size)
+        int sampling_rate()
+        void sampling_rate(int new_sampling_rate)
+        void residual_frame(int synth_size, double* synth,
+                            int original_size, double* original,
+                            int residual_size, double* residual)
+        void find_residual(int synth_size, double* synth,
+                           int original_size, double* original,
+                           int residual_size, double* residual)
+        void synth_frame(c_Frame* frame)
+        vector[c_Frame*] synth(vector[c_Frame*] frames)
+
 
 cdef class Peak:
     cdef c_Peak* thisptr
@@ -355,3 +387,99 @@ cdef class PartialTracking:
             f.set_frame(output_frames[i])
             partial_frames.append(f)
         return partial_frames
+
+
+cdef class Synthesis:
+    cdef c_Synthesis* thisptr
+
+    def __cinit__(self): self.thisptr = new c_Synthesis()
+    def __dealloc__(self): del self.thisptr
+
+    property sampling_rate:
+        def __get__(self): return self.thisptr.sampling_rate()
+        def __set__(self, int i): self.thisptr.sampling_rate(i)
+
+    property frame_size:
+        def __get__(self): return self.thisptr.frame_size()
+        def __set__(self, int i): self.thisptr.frame_size(i)
+
+    property hop_size:
+        def __get__(self): return self.thisptr.hop_size()
+        def __set__(self, int i): self.thisptr.hop_size(i)
+
+    property max_partials:
+        def __get__(self): return self.thisptr.max_partials()
+        def __set__(self, int i): self.thisptr.max_partials(i)
+
+    def synth_frame(self, Frame frame not None):
+        self.thisptr.synth_frame(frame.thisptr)
+        return frame.audio
+
+    def synth(self, frames):
+        cdef vector[c_Frame*] c_frames
+        for frame in frames:
+            c_frames.push_back((<Frame>frame).thisptr)
+        cdef vector[c_Frame*] output_frames = self.thisptr.synth(c_frames)
+        cdef np.ndarray[dtype_t, ndim=1] output = np.zeros(
+            output_frames.size() * self.thisptr.hop_size()
+        )
+        cdef np.npy_intp shape[1]
+        shape[0] = <np.npy_intp> self.thisptr.hop_size()
+        for i in range(output_frames.size()):
+            frame_audio = np.PyArray_SimpleNewFromData(1, shape, np.NPY_DOUBLE, output_frames[i].synth())
+            output[i * self.thisptr.hop_size():(i + 1) * self.thisptr.hop_size()] = frame_audio
+        return output
+
+
+cdef class Residual:
+    cdef c_Residual* thisptr
+
+    def __cinit__(self): self.thisptr = new c_Residual()
+    def __dealloc__(self): del self.thisptr
+
+    property sampling_rate:
+        def __get__(self): return self.thisptr.sampling_rate()
+        def __set__(self, int i): self.thisptr.sampling_rate(i)
+
+    property frame_size:
+        def __get__(self): return self.thisptr.frame_size()
+        def __set__(self, int i): self.thisptr.frame_size(i)
+
+    property hop_size:
+        def __get__(self): return self.thisptr.hop_size()
+        def __set__(self, int i): self.thisptr.hop_size(i)
+
+    def residual_frame(self, np.ndarray[dtype_t, ndim=1] synth,
+                       np.ndarray[dtype_t, ndim=1] original):
+        cdef np.ndarray[dtype_t, ndim=1] residual = np.zeros(len(synth))
+        self.thisptr.residual_frame(len(synth), <double*> synth.data,
+                                    len(original), <double*> original.data,
+                                    len(residual), <double*> residual.data)
+        return residual
+
+    def find_residual(self, np.ndarray[dtype_t, ndim=1] synth,
+                      np.ndarray[dtype_t, ndim=1] original):
+        cdef np.ndarray[dtype_t, ndim=1] residual = np.zeros(len(synth))
+        self.thisptr.find_residual(len(synth), <double*> synth.data,
+                                   len(original), <double*> original.data,
+                                   len(residual), <double*> residual.data)
+        return residual
+
+    def synth_frame(self, Frame frame not None):
+        self.thisptr.synth_frame(frame.thisptr)
+        return frame.audio
+
+    def synth(self, frames):
+        cdef vector[c_Frame*] c_frames
+        for frame in frames:
+            c_frames.push_back((<Frame>frame).thisptr)
+        cdef vector[c_Frame*] output_frames = self.thisptr.synth(c_frames)
+        cdef np.ndarray[dtype_t, ndim=1] output = np.zeros(
+            output_frames.size() * self.thisptr.hop_size()
+        )
+        cdef np.npy_intp shape[1]
+        shape[0] = <np.npy_intp> self.thisptr.hop_size()
+        for i in range(output_frames.size()):
+            frame_audio = np.PyArray_SimpleNewFromData(1, shape, np.NPY_DOUBLE, output_frames[i].synth())
+            output[i * self.thisptr.hop_size():(i + 1) * self.thisptr.hop_size()] = frame_audio
+        return output
diff --git a/src/simpl/base.cpp b/src/simpl/base.cpp
index 9c70403..9d5f7ad 100644
--- a/src/simpl/base.cpp
+++ b/src/simpl/base.cpp
@@ -428,3 +428,116 @@ Frames PartialTracking::find_partials(Frames frames) {
     _frames = frames;
     return _frames;
 }
+
+
+// ---------------------------------------------------------------------------
+// Synthesis
+// ---------------------------------------------------------------------------
+Synthesis::Synthesis() {
+    _frame_size = 512;
+    _hop_size = 512;
+    _max_partials = 100;
+    _sampling_rate = 44100;
+}
+
+int Synthesis::frame_size() {
+    return _frame_size;
+}
+
+void Synthesis::frame_size(int new_frame_size) {
+    _frame_size = new_frame_size;
+}
+
+int Synthesis::hop_size() {
+    return _hop_size;
+}
+
+void Synthesis::hop_size(int new_hop_size) {
+    _hop_size = new_hop_size;
+}
+
+int Synthesis::max_partials() {
+    return _max_partials;
+}
+
+void Synthesis::max_partials(int new_max_partials) {
+    _max_partials = new_max_partials;
+}
+
+int Synthesis::sampling_rate() {
+    return _sampling_rate;
+}
+
+void Synthesis::sampling_rate(int new_sampling_rate) {
+    _sampling_rate = new_sampling_rate;
+}
+
+void Synthesis::synth_frame(Frame* frame) {
+}
+
+Frames Synthesis::synth(Frames frames) {
+    for(int i = 0; i < frames.size(); i++) {
+        synth_frame(frames[i]);
+    }
+    return frames;
+}
+
+
+// ---------------------------------------------------------------------------
+// Residual
+// ---------------------------------------------------------------------------
+Residual::Residual() {
+    _frame_size = 512;
+    _hop_size = 512;
+    _sampling_rate = 44100;
+}
+
+int Residual::frame_size() {
+    return _frame_size;
+}
+
+void Residual::frame_size(int new_frame_size) {
+    _frame_size = new_frame_size;
+}
+
+int Residual::hop_size() {
+    return _hop_size;
+}
+
+void Residual::hop_size(int new_hop_size) {
+    _hop_size = new_hop_size;
+}
+
+int Residual::sampling_rate() {
+    return _sampling_rate;
+}
+
+void Residual::sampling_rate(int new_sampling_rate) {
+    _sampling_rate = new_sampling_rate;
+}
+
+void Residual::residual_frame(int synth_size, sample* synth,
+                              int original_size, sample* original,
+                              int residual_size, sample* residual) {
+}
+
+void Residual::find_residual(int synth_size, sample* synth,
+                             int original_size, sample* original,
+                             int residual_size, sample* residual) {
+    for(int i = 0; i < synth_size; i += _hop_size) {
+        residual_frame(_hop_size, &synth[i],
+                       _hop_size, &original[i],
+                       _hop_size, &residual[i]);
+    }
+}
+
+void Residual::synth_frame(Frame* frame) {
+}
+
+// Calculate and return a synthesised residual signal
+Frames Residual::synth(Frames frames) {
+    for(int i = 0; i < frames.size(); i++) {
+        synth_frame(frames[i]);
+    }
+    return frames;
+}
diff --git a/src/simpl/base.h b/src/simpl/base.h
index 385e2a9..bdb0d96 100644
--- a/src/simpl/base.h
+++ b/src/simpl/base.h
@@ -7,7 +7,7 @@
 
 using namespace std;
 
-namespace simpl 
+namespace simpl
 {
 
 typedef double sample;
@@ -68,10 +68,10 @@ typedef std::vector<Partial*> Partials;
 
 // ---------------------------------------------------------------------------
 // Frame
-// 
+//
 // Represents a frame of audio information.
-// This can be: - raw audio samples 
-//              - an unordered list of sinusoidal peaks 
+// This can be: - raw audio samples
+//              - an unordered list of sinusoidal peaks
 //              - an ordered list of partials
 //              - synthesised audio samples
 //              - residual samples
@@ -129,7 +129,7 @@ typedef std::vector<Frame*> Frames;
 
 // ---------------------------------------------------------------------------
 // PeakDetection
-// 
+//
 // Detect spectral peaks
 // ---------------------------------------------------------------------------
 
@@ -214,6 +214,67 @@ class PartialTracking {
         virtual Frames find_partials(Frames frames);
 };
 
+
+// ---------------------------------------------------------------------------
+// Synthesis
+//
+// Synthesise audio from spectral analysis data
+// ---------------------------------------------------------------------------
+
+class Synthesis {
+    private:
+        int _frame_size;
+        int _hop_size;
+        int _max_partials;
+        int _sampling_rate;
+
+    public:
+        Synthesis();
+        int frame_size();
+        void frame_size(int new_frame_size);
+        int hop_size();
+        void hop_size(int new_hop_size);
+        int max_partials();
+        void max_partials(int new_max_partials);
+        int sampling_rate();
+        void sampling_rate(int new_sampling_rate);
+
+        virtual void synth_frame(Frame* frame);
+        virtual Frames synth(Frames frames);
+};
+
+// ---------------------------------------------------------------------------
+// Residual
+//
+// Calculate a residual signal
+// ---------------------------------------------------------------------------
+
+class Residual {
+    private:
+        int _frame_size;
+        int _hop_size;
+        int _sampling_rate;
+
+    public:
+        Residual();
+        int frame_size();
+        void frame_size(int new_frame_size);
+        int hop_size();
+        void hop_size(int new_hop_size);
+        int sampling_rate();
+        void sampling_rate(int new_sampling_rate);
+
+        virtual void residual_frame(int synth_size, sample* synth,
+                                    int original_size, sample* original,
+                                    int residual_size, sample* residual);
+        virtual void find_residual(int synth_size, sample* synth,
+                                   int original_size, sample* original,
+                                   int residual_size, sample* residual);
+
+        virtual void synth_frame(Frame* frame);
+        virtual Frames synth(Frames frames);
+};
+
 } // end of namespace Simpl
 
 #endif
diff --git a/tests/test_base.py b/tests/test_base.py
index 738d530..ce4297e 100644
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -4,6 +4,13 @@ import scipy.io.wavfile as wavfile
 from nose.tools import assert_almost_equals
 import simpl.base as base
 
+float_precision = 5
+frame_size = 512
+hop_size = 512
+audio_path = os.path.join(
+    os.path.dirname(__file__), 'audio/flute.wav'
+)
+
 
 class TestFrame(object):
     def test_buffers(self):
@@ -39,8 +46,10 @@ class TestFrame(object):
         f.add_peak(p)
 
         assert f.num_peaks == 1
-        assert f.peak(0).amplitude == p.amplitude
-        assert f.peaks[0].amplitude == p.amplitude
+        assert_almost_equals(f.peak(0).amplitude, p.amplitude,
+                             float_precision)
+        assert_almost_equals(f.peaks[0].amplitude, p.amplitude,
+                             float_precision)
 
         f.clear()
         assert f.num_peaks == 0
@@ -56,21 +65,16 @@ class TestFrame(object):
         p.phase = 0.0
 
         f.partial(0, p)
-        assert f.partial(0).amplitude == p.amplitude
-        assert f.partial(0).frequency == p.frequency
+        assert_almost_equals(f.partial(0).amplitude, p.amplitude,
+                             float_precision)
+        assert_almost_equals(f.partial(0).frequency, p.frequency,
+                             float_precision)
 
 
 class TestPeakDetection(object):
-    float_precision = 5
-    frame_size = 512
-    hop_size = 512
-    audio_path = os.path.join(
-        os.path.dirname(__file__), 'audio/flute.wav'
-    )
-
     @classmethod
     def setup_class(cls):
-        cls.audio = wavfile.read(cls.audio_path)[1]
+        cls.audio = wavfile.read(audio_path)[1]
         cls.audio = np.asarray(cls.audio, dtype=np.double)
         cls.audio /= np.max(cls.audio)
 
@@ -78,21 +82,14 @@ class TestPeakDetection(object):
         pd = base.PeakDetection()
         pd.find_peaks(self.audio)
 
-        assert len(pd.frames) == len(self.audio) / self.hop_size
+        assert len(pd.frames) == len(self.audio) / hop_size
         assert len(pd.frames[0].peaks) == 0
 
 
 class TestPartialTracking(object):
-    float_precision = 5
-    frame_size = 512
-    hop_size = 512
-    audio_path = os.path.join(
-        os.path.dirname(__file__), 'audio/flute.wav'
-    )
-
     @classmethod
     def setup_class(cls):
-        cls.audio = wavfile.read(cls.audio_path)[1]
+        cls.audio = wavfile.read(audio_path)[1]
         cls.audio = np.asarray(cls.audio, dtype=np.double)
         cls.audio /= np.max(cls.audio)
 
@@ -103,5 +100,47 @@ class TestPartialTracking(object):
         pt = base.PartialTracking()
         frames = pt.find_partials(frames)
 
-        assert len(frames) == len(self.audio) / self.hop_size
+        assert len(frames) == len(self.audio) / hop_size
         assert len(frames[0].partials) == 100
+
+
+class TestSynthesis(object):
+    @classmethod
+    def setup_class(cls):
+        cls.audio = wavfile.read(audio_path)[1]
+        cls.audio = np.asarray(cls.audio, dtype=np.double)
+        cls.audio /= np.max(cls.audio)
+
+    def test_synthesis(self):
+        pd = base.PeakDetection()
+        frames = pd.find_peaks(self.audio)
+
+        pt = base.PartialTracking()
+        frames = pt.find_partials(frames)
+
+        s = base.Synthesis()
+        synth_audio = s.synth(frames)
+
+        assert len(synth_audio)
+
+
+class TestResidual(object):
+    @classmethod
+    def setup_class(cls):
+        cls.audio = wavfile.read(audio_path)[1]
+        cls.audio = np.asarray(cls.audio, dtype=np.double)
+        cls.audio /= np.max(cls.audio)
+
+    def test_synthesis(self):
+        pd = base.PeakDetection()
+        frames = pd.find_peaks(self.audio)
+
+        pt = base.PartialTracking()
+        frames = pt.find_partials(frames)
+
+        s = base.Synthesis()
+        synth_audio = s.synth(frames)
+
+        r = base.Residual()
+        residual_audio = r.find_residual(synth_audio, self.audio)
+        assert len(residual_audio)