Added MQSynthesis, fixed a few bugs

7 files changed, 283 insertions, 206 deletions

diff --git a/__init__.py b/__init__.py
index 44f2c37..cd46097 100644
--- a/__init__.py
+++ b/__init__.py
@@ -18,10 +18,10 @@ from basetypes import Peak, PeakDetection, Partial, PartialTracking, Synthesis,
 from basetypes import compare_peak_amps, compare_peak_freqs
 from sndobj import SndObjPeakDetection, SndObjPartialTracking, SndObjSynthesis
 from sms import SMSPeakDetection, SMSPartialTracking, SMSSynthesis, SMSResidual
-from mq import MQPeakDetection, MQPartialTracking
+from mq import MQPeakDetection, MQPartialTracking, MQSynthesis
+from plot import plot_peaks, plot_partials
 
 import numpy
-#float = numpy.float64
 
 def array (n, type=float):
     return(numpy.array(n, dtype=type))
diff --git a/basetypes.py b/basetypes.py
index b2f0ae8..28311ad 100644
--- a/basetypes.py
+++ b/basetypes.py
@@ -267,15 +267,29 @@ class PartialTracking(object):
 class Synthesis(object):
     "Synthesise audio from spectral analysis data"
     def __init__(self):
+        self._frame_size = 512
         self._hop_size = 512
         self._max_partials = 100
+        self._sampling_rate = 44100
         
     # properties
+    frame_size = property(lambda self: self.get_frame_size(),
+                          lambda self, x: self.set_frame_size(x))
     hop_size = property(lambda self: self.get_hop_size(),
                         lambda self, x: self.set_hop_size(x))
     max_partials = property(lambda self: self.get_max_partials(),
                             lambda self, x: self.set_max_partials(x))
+    max_partials = property(lambda self: self.get_max_partials(),
+                            lambda self, x: self.set_max_partials(x))
+    sampling_rate = property(lambda self: self.get_sampling_rate(),
+                             lambda self, x: self.set_sampling_rate(x))
         
+    def get_frame_size(self):
+        return self._frame_size
+    
+    def set_frame_size(self, frame_size):
+        self._frame_size = frame_size
+
     def get_hop_size(self):
         return self._hop_size
     
@@ -288,6 +302,12 @@ class Synthesis(object):
     def set_max_partials(self, num_partials):
         self._max_partials = num_partials
         
+    def get_sampling_rate(self):
+        return self._sampling_rate
+    
+    def set_sampling_rate(self, sampling_rate):
+        self._sampling_rate = sampling_rate
+
     def synth_frame(self, partials):
         "Synthesises a frame of audio, given a list of peaks from tracks"
         raise Exception("NotYetImplemented")
@@ -336,3 +356,4 @@ class Residual(object):
     def synth(self, synth, original):
         "Calculate and return a synthesised residual signal"
         raise Exception("NotYetImplemented")
+
diff --git a/mq.py b/mq.py
index b5835dc..7ef7b36 100644
--- a/mq.py
+++ b/mq.py
@@ -63,7 +63,7 @@ class MQPeakDetection(simpl.PeakDetection):
                 p = simpl.Peak()
                 p.amplitude = current_mag
                 p.frequency = (bin - 1) * self._fundamental
-                p.phase = np.angle(spectrum[bin-1])
+                p.phase = np.angle(f[bin-1])
                 current_peaks.append(p)
             prev_mag = current_mag
             current_mag = next_mag
@@ -147,7 +147,7 @@ class MQPartialTracking(simpl.PartialTracking):
         for peak in self._current_frame:
             if peak.is_start_of_partial():
                 partial = simpl.Partial()
-                partial.starting_frame = frame_number-1
+                partial.starting_frame = frame_number
                 partial.add_peak(peak)
                 self.partials.append(partial)
             match = self._find_closest_match(peak, frame)
@@ -179,7 +179,7 @@ class MQPartialTracking(simpl.PartialTracking):
                 # create a new track by adding a peak in the current frame, matched to p,
                 # with amplitude 0
                 partial = simpl.Partial()
-                partial.starting_frame = frame_number-1
+                partial.starting_frame = frame_number
                 new_peak = simpl.Peak()
                 new_peak.amplitude = 0
                 new_peak.frequency = p.frequency
@@ -189,3 +189,57 @@ class MQPartialTracking(simpl.PartialTracking):
         self._current_frame = frame        
         return frame_partials
 
+
+class MQSynthesis(simpl.Synthesis):
+    def __init__(self):
+        simpl.Synthesis.__init__(self)
+        self._current_frame = simpl.zeros(self.frame_size)
+
+    def hz_to_radians(self, frequency):
+        if not frequency:
+            return 0.0
+        else:
+            return (frequency * 2.0 * np.pi) / self.sampling_rate
+    
+    def synth_frame(self, peaks):
+        "Synthesises a frame of audio, given a list of peaks from tracks"
+        self._current_frame *= 0.0
+
+        for p in peaks:
+            # get values for last amplitude, frequency and phase
+            # these are the initial values of the instantaneous amplitude/frequency/phase
+            current_freq = self.hz_to_radians(p.frequency)
+            if not p.previous_peak:
+                prev_amp = 0.0
+                prev_freq = current_freq 
+                prev_phase = p.phase - (current_freq * self.frame_size)
+                while prev_phase >= np.pi: prev_phase -= 2.0 * np.pi 
+                while prev_phase < -np.pi: prev_phase += 2.0 * np.pi
+            else:
+                prev_amp = p.previous_peak.amplitude
+                prev_freq = self.hz_to_radians(p.previous_peak.frequency)
+                prev_phase = p.previous_peak.phase
+
+            # amplitudes are linearly interpolated between frames
+            inst_amp = prev_amp
+            amp_inc = (p.amplitude - prev_amp) / self.frame_size
+
+            # freqs/phases are calculated by cubic interpolation
+            freq_diff = current_freq - prev_freq
+            x = ((prev_phase + (prev_freq * self.frame_size) - p.phase) +
+                 (freq_diff * (self.frame_size / 2.0)))
+            x /= (2.0 * np.pi)
+            m = int(np.round(x))
+            phase_diff = p.phase - prev_phase - (prev_freq * self.frame_size) + (2.0 * np.pi * m)
+            alpha = ((3.0 / (self.frame_size**2)) * phase_diff) - (freq_diff / self.frame_size)
+            beta = ((-2.0 / (self.frame_size**3)) * phase_diff) + (freq_diff / (self.frame_size**2))
+
+            # calculate output samples
+            for i in range(self.frame_size):
+                inst_amp += amp_inc
+                inst_phase = prev_phase + (prev_freq * i) + (alpha * (i**2)) + (beta * (i**3))
+                #print inst_phase
+                self._current_frame[i] += (2.0 * inst_amp) * np.cos(inst_phase)
+
+        return self._current_frame
+
diff --git a/readme.txt b/readme.txt
index 6fbfb52..864749a 100644
--- a/readme.txt
+++ b/readme.txt
@@ -116,4 +116,5 @@ sms:
 mq:
 - fix peak detection module
 - window size should change depending on detected pitch. No pitch detection algorithm described in the
-  original paper, so can use the SMS algorithm for now.
\ No newline at end of file
+  original paper, so can use the SMS algorithm for now.
+
diff --git a/sms.py b/sms.py
index 4c4584c..39c3c2a 100644
--- a/sms.py
+++ b/sms.py
@@ -282,6 +282,8 @@ class SMSSynthesis(simpl.Synthesis):
         pysms.sms_init()
         self._synth_params = pysms.SMS_SynthParams() 
         self._synth_params.iDetSynthType = pysms.SMS_DET_SIN
+        # use the default simpl hop size instead of the default SMS hop size
+        self._synth_params.sizeHop = self._hop_size 
         pysms.sms_initSynth(self._synth_params)
         self._current_frame = simpl.zeros(self.hop_size)
         self._analysis_frame = pysms.SMS_Data()
@@ -295,8 +297,6 @@ class SMSSynthesis(simpl.Synthesis):
         SMSSynthesis._instances -= 1
         
     # properties
-    sampling_rate = property(lambda self: self.get_sampling_rate(),
-                             lambda self, x: self.set_sampling_rate(x))
     synthesis_type = property(lambda self: self.get_synthesis_type(),
                               lambda self, x: self.set_synthesis_type(x))
     num_stochastic_coeffs = property(lambda self: self.get_num_stochastic_coeffs(),
@@ -402,3 +402,4 @@ class SMSResidual(simpl.Residual):
         if pysms.sms_findResidual(synth, original, residual, self._analysis_params) == -1:
             raise Exception("Residual error: Synthesised audio and original audio have different lengths")
         return residual
+
diff --git a/sms/sineSynth.c b/sms/sineSynth.c
index aaca41d..58768bd 100644
--- a/sms/sineSynth.c
+++ b/sms/sineSynth.c
@@ -59,7 +59,7 @@ static void SinePhaSynth (sfloat fFreq, sfloat fMag, sfloat fPhase,
             (pLastFrame->pFSinFreq[iTrack] * sizeBuffer);
     fPhase = fTmp - floor(fTmp / TWO_PI) * TWO_PI;
   }
-  
+
   /* caculate the instantaneous amplitude */
   fMagIncr = (fMag - pLastFrame->pFSinAmp[iTrack]) / sizeBuffer;
   fInstMag = pLastFrame->pFSinAmp[iTrack];
diff --git a/sndobj/SinSyn.cpp b/sndobj/SinSyn.cpp
index 7452b95..7d31e15 100644
--- a/sndobj/SinSyn.cpp
+++ b/sndobj/SinSyn.cpp
@@ -1,4 +1,4 @@
- 
+
 ////////////////////////////////////////////////////////////////////////
 // This file is part of the SndObj library
 //
@@ -24,86 +24,86 @@
 
 SinSyn::SinSyn(){
 
-  m_factor = m_vecsize/m_sr;
-  m_facsqr = m_factor*m_factor;
-  m_ptable = 0;
-  m_size = 0;
-  m_LoTWOPI = 0.f;
-  m_maxtracks = 0;
-  m_freqs = 0;
-  m_amps = 0;
-  m_phases = 0;
-  m_trackID = 0;
-  m_scale =  0.f;
-  m_ratio = 0.f;
-  m_tracks = 0;
-  AddMsg("max tracks", 21);
-  AddMsg("scale", 23);
-  AddMsg("table", 24);
+    m_factor = m_vecsize/m_sr;
+    m_facsqr = m_factor*m_factor;
+    m_ptable = 0;
+    m_size = 0;
+    m_LoTWOPI = 0.f;
+    m_maxtracks = 0;
+    m_freqs = 0;
+    m_amps = 0;
+    m_phases = 0;
+    m_trackID = 0;
+    m_scale =  0.f;
+    m_ratio = 0.f;
+    m_tracks = 0;
+    AddMsg("max tracks", 21);
+    AddMsg("scale", 23);
+    AddMsg("table", 24);
 }
 
 SinSyn::SinSyn(SinAnal* input, int maxtracks, Table* table, 
-	       double scale, int vecsize, double sr)	  
-  :SndObj(input, vecsize, sr){
-
-  m_ptable = table;
-  m_size = m_ptable->GetLen();
-  m_LoTWOPI = m_size/TWOPI;
-
-  m_factor = m_vecsize/m_sr;
-  m_facsqr = m_factor*m_factor;
-  m_maxtracks = maxtracks;
-  m_tracks = 0;
-  m_scale = scale;
-  m_input = input;
-  m_freqs = new double[m_maxtracks];
-  m_amps = new double[m_maxtracks];
-  m_phases = new double[m_maxtracks];
-  m_trackID = new int[m_maxtracks];
-
-  memset(m_phases, 0, sizeof(double)*m_maxtracks);
-
-  m_incr = 0.f;
-  m_ratio = m_size/m_sr;
-  AddMsg("max tracks", 21);
-  AddMsg("scale", 23);
-  AddMsg("table", 24);
-  AddMsg("timescale", 24);
-  memset(m_trackID, 0, sizeof(int));
+        double scale, int vecsize, double sr)      
+:SndObj(input, vecsize, sr){
+
+    m_ptable = table;
+    m_size = m_ptable->GetLen();
+    m_LoTWOPI = m_size/TWOPI;
+
+    m_factor = m_vecsize/m_sr;
+    m_facsqr = m_factor*m_factor;
+    m_maxtracks = maxtracks;
+    m_tracks = 0;
+    m_scale = scale;
+    m_input = input;
+    m_freqs = new double[m_maxtracks];
+    m_amps = new double[m_maxtracks];
+    m_phases = new double[m_maxtracks];
+    m_trackID = new int[m_maxtracks];
+
+    memset(m_phases, 0, sizeof(double)*m_maxtracks);
+
+    m_incr = 0.f;
+    m_ratio = m_size/m_sr;
+    AddMsg("max tracks", 21);
+    AddMsg("scale", 23);
+    AddMsg("table", 24);
+    AddMsg("timescale", 24);
+    memset(m_trackID, 0, sizeof(int));
 
 }
 
 SinSyn::~SinSyn(){
 
-  delete[] m_freqs;  
-  delete[] m_amps;  
-  delete[] m_phases;  
-  delete[] m_trackID;
+    delete[] m_freqs;  
+    delete[] m_amps;  
+    delete[] m_phases;  
+    delete[] m_trackID;
 
 }
 
-void
+    void
 SinSyn::SetTable(Table *table)
 {
-  m_ptable = table;
-  m_size = m_ptable->GetLen();
-  m_LoTWOPI = m_size/TWOPI;
-  m_ratio = m_size/m_sr;
+    m_ptable = table;
+    m_size = m_ptable->GetLen();
+    m_LoTWOPI = m_size/TWOPI;
+    m_ratio = m_size/m_sr;
 }
 
 int
 SinSyn::Connect(char* mess, void* input){
 
-  switch (FindMsg(mess)){
+    switch (FindMsg(mess)){
+
+        case 24:
+            SetTable((Table *) input);
+            return 1;
 
-  case 24:
-    SetTable((Table *) input);
-    return 1;
+        default:
+            return SndObj::Connect(mess,input);
 
-  default:
-    return SndObj::Connect(mess,input);
-     
-  }
+    }
 
 }
 
@@ -111,158 +111,158 @@ SinSyn::Connect(char* mess, void* input){
 int
 SinSyn::Set(char* mess, double value){
 
-  switch(FindMsg(mess)){
+    switch(FindMsg(mess)){
 
-  case 21:
-    SetMaxTracks((int)value);
-    return 1;
+        case 21:
+            SetMaxTracks((int)value);
+            return 1;
 
-  case 23:
-    SetScale(value);
-    return 1;
-	
-  default:
-    return SndObj::Set(mess, value);
+        case 23:
+            SetScale(value);
+            return 1;
 
-  }
+        default:
+            return SndObj::Set(mess, value);
+
+    }
 }
 
 
 void
 SinSyn::SetMaxTracks(int maxtracks){
 
-  if(m_maxtracks){
+    if(m_maxtracks){
 
-    delete[] m_freqs;  
-    delete[] m_amps;  
-    delete[] m_phases;  
-    delete[] m_trackID;
+        delete[] m_freqs;  
+        delete[] m_amps;  
+        delete[] m_phases;  
+        delete[] m_trackID;
 
-  }
+    }
 
-  m_maxtracks = maxtracks;
-  m_freqs = new double[m_maxtracks];
-  m_amps = new double[m_maxtracks];
-  m_phases = new double[m_maxtracks];
-  m_trackID = new int[m_maxtracks];
+    m_maxtracks = maxtracks;
+    m_freqs = new double[m_maxtracks];
+    m_amps = new double[m_maxtracks];
+    m_phases = new double[m_maxtracks];
+    m_trackID = new int[m_maxtracks];
 
 }
 
 short
 SinSyn::DoProcess() {
-	
-  if(m_input){
-		
-    double ampnext,amp,freq, freqnext, phase,phasenext;
-    double a2, a3, phasediff, cph;
-    int i3, i, j, ID, track;
-    int notcontin = 0;
-    bool contin = false;
-    int oldtracks = m_tracks;
-    double* tab = m_ptable->GetTable(); 
-    if((m_tracks = ((SinAnal *)m_input)->GetTracks()) >
-       m_maxtracks) m_tracks = m_maxtracks;
-		
-    memset(m_output, 0, sizeof(double)*m_vecsize);
-   		
-    // for each track
-    i = j = 0;
-    while(i < m_tracks*3){			
-      i3 = i/3;
-      ampnext =  m_input->Output(i)*m_scale;
-      freqnext = m_input->Output(i+1)*TWOPI; 
-      phasenext = m_input->Output(i+2);
-      ID = ((SinAnal *)m_input)->GetTrackID(i3);
-
-      j = i3+notcontin;
-	
-      if(i3 < oldtracks-notcontin){
-          
-	if(m_trackID[j]==ID){	
-	  // if this is a continuing track  	
-	  track = j;
-	  contin = true;	
-	  freq = m_freqs[track];
-	  phase = m_phases[track];
-	  amp = m_amps[track];
-					
-	}
-	else {
-	  // if this is  a dead track
-	  contin = false;
-	  track = j;
-	  freqnext = freq = m_freqs[track];
-	  phase = m_phases[track];
-	  phasenext = phase + freq*m_factor;
-	  amp = m_amps[track]; 
-	  ampnext = 0.f;
-	}
-      }
-			
-      else{ 
-	// new tracks
-	contin = true;
-	track = -1;
-	freq = freqnext;
-	phase = phasenext - freq*m_factor;
-	amp = 0.f;
-      }
-			
-      // phasediff
-      phasediff = phasenext - phase;		
-      while(phasediff >= PI) phasediff -= TWOPI;
-      while(phasediff < -PI) phasediff += TWOPI;
-      // update phasediff to match the freq
-      cph = ((freq+freqnext)*m_factor/2. - phasediff)/TWOPI;
-      phasediff += TWOPI * Ftoi(cph + 0.5);
-      // interpolation coefs
-      a2 = 3./m_facsqr * (phasediff - m_factor/3.*(2*freq+freqnext));
-      a3 = 1./(3*m_facsqr)  * (freqnext - freq - 2*a2*m_factor);
-
-      // interpolation resynthesis loop	
-      double inc1, inc2, a, ph, cnt, frac;
-      int ndx;
-      a = amp;
-      ph = phase;
-      cnt = 0;
-      inc1 = (ampnext - amp)/m_vecsize;
-      inc2 = 1/m_sr;  
-      for(m_vecpos=0; m_vecpos < m_vecsize; m_vecpos++){
-		
-	if(m_enable) {    
- 	  // table lookup oscillator
-	  ph *= m_LoTWOPI;
-	  while(ph < 0) ph += m_size;  
-	  while(ph >= m_size) ph -= m_size;
-	  ndx = Ftoi(ph);
-	  frac = ph - ndx;
-	  m_output[m_vecpos] += a*(tab[ndx] + (tab[ndx+1] - tab[ndx])*frac); 
-	  a += inc1;
-	  cnt += inc2;
-	  ph = phase + cnt*(freq + cnt*(a2 + a3*cnt));
-					
-	}
-	else m_output[m_vecpos] = 0.f;		
-      }
-     
-      // keep amp, freq, and phase values for next time
-      if(contin){
-	m_amps[i3] = ampnext;
-	m_freqs[i3] = freqnext;
-	while(phasenext < 0) phasenext += TWOPI;
-	while(phasenext >= TWOPI) phasenext -= TWOPI;
-	m_phases[i3] = phasenext;
-	m_trackID[i3] = ID;    
-	i += 3;
-      } else notcontin++;
-      } 
-   
-    return 1;
-  }
-  else {
-    m_error  = 1;
-    return 0;
-  }
- 
+
+    if(m_input){
+
+        double ampnext,amp,freq, freqnext, phase,phasenext;
+        double a2, a3, phasediff, cph;
+        int i3, i, j, ID, track;
+        int notcontin = 0;
+        bool contin = false;
+        int oldtracks = m_tracks;
+        double* tab = m_ptable->GetTable(); 
+        if((m_tracks = ((SinAnal *)m_input)->GetTracks()) >
+                m_maxtracks) m_tracks = m_maxtracks;
+
+        memset(m_output, 0, sizeof(double)*m_vecsize);
+
+        // for each track
+        i = j = 0;
+        while(i < m_tracks*3){          
+            i3 = i/3;
+            ampnext =  m_input->Output(i)*m_scale;
+            freqnext = m_input->Output(i+1)*TWOPI; 
+            phasenext = m_input->Output(i+2);
+            ID = ((SinAnal *)m_input)->GetTrackID(i3);
+
+            j = i3+notcontin;
+
+            if(i3 < oldtracks-notcontin){
+
+                if(m_trackID[j]==ID){   
+                    // if this is a continuing track      
+                    track = j;
+                    contin = true;    
+                    freq = m_freqs[track];
+                    phase = m_phases[track];
+                    amp = m_amps[track];
+
+                }
+                else {
+                    // if this is  a dead track
+                    contin = false;
+                    track = j;
+                    freqnext = freq = m_freqs[track];
+                    phase = m_phases[track];
+                    phasenext = phase + freq*m_factor;
+                    amp = m_amps[track]; 
+                    ampnext = 0.f;
+                }
+            }
+
+            else{ 
+                // new tracks
+                contin = true;
+                track = -1;
+                freq = freqnext;
+                phase = phasenext - freq*m_factor;
+                amp = 0.f;
+            }
+
+            // phasediff
+            phasediff = phasenext - phase;        
+            while(phasediff >= PI) phasediff -= TWOPI;
+            while(phasediff < -PI) phasediff += TWOPI;
+            // update phasediff to match the freq
+            cph = ((freq+freqnext)*m_factor/2. - phasediff)/TWOPI;
+            phasediff += TWOPI * Ftoi(cph + 0.5);
+            // interpolation coefs
+            a2 = 3./m_facsqr * (phasediff - m_factor/3.*(2*freq+freqnext));
+            a3 = 1./(3*m_facsqr)  * (freqnext - freq - 2*a2*m_factor);
+
+            // interpolation resynthesis loop 
+            double inc1, inc2, a, ph, cnt, frac;
+            int ndx;
+            a = amp;
+            ph = phase;
+            cnt = 0;
+            inc1 = (ampnext - amp)/m_vecsize;
+            inc2 = 1/m_sr;  
+            for(m_vecpos=0; m_vecpos < m_vecsize; m_vecpos++){
+
+                if(m_enable) {    
+                    // table lookup oscillator
+                    ph *= m_LoTWOPI;
+                    while(ph < 0) ph += m_size;  
+                    while(ph >= m_size) ph -= m_size;
+                    ndx = Ftoi(ph);
+                    frac = ph - ndx;
+                    m_output[m_vecpos] += a*(tab[ndx] + (tab[ndx+1] - tab[ndx])*frac); 
+                    a += inc1;
+                    cnt += inc2;
+                    ph = phase + cnt*(freq + cnt*(a2 + a3*cnt));
+
+                }
+                else m_output[m_vecpos] = 0.f;      
+            }
+
+            // keep amp, freq, and phase values for next time
+            if(contin){
+                m_amps[i3] = ampnext;
+                m_freqs[i3] = freqnext;
+                while(phasenext < 0) phasenext += TWOPI;
+                while(phasenext >= TWOPI) phasenext -= TWOPI;
+                m_phases[i3] = phasenext;
+                m_trackID[i3] = ID;    
+                i += 3;
+            } else notcontin++;
+        } 
+
+        return 1;
+    }
+    else {
+        m_error  = 1;
+        return 0;
+    }
+
 }