Added multirate

src/sppysound/analysis/Analysis.py

@@ -133,11 +133,10 @@ class Analysis(object):
         # get all valid frames from current grain
         frames = frames[selection & valid_inds]
 
+        return formatter(frames)
         #if less than half the frames are valid then the grain is not valid.
         if frames.size < valid_inds[selection].nonzero()[0].size/2:
             return np.nan
-        else:
-            return formatter(frames)
 
     def analysis_formatter(self, frames, selection, format):
         """Calculate the average analysis value of the grain using the match format specified."""

src/sppysound/analysis/F0Analysis.py

@@ -182,25 +182,20 @@ class F0Analysis(Analysis):
                 except Warning:
                     pass
 
-            Z = feature_zcr(Gamma)
-            if Z > 0.15:
+            # compute T0 and harmonic ratio:
+            if np.isnan(Gamma).any():
                 HR = np.nan
                 f0 = np.nan
             else:
-                # compute T0 and harmonic ratio:
-                if np.isnan(Gamma).any():
-                    HR=np.nan
+                blag = np.argmax(Gamma)
+                HR = Gamma[blag]
+                interp, HR = parabolic(Gamma, blag)
+                if not interp:
                     f0 = np.nan
+                    HR = np.nan
                 else:
-                    blag = np.argmax(Gamma)
-                    HR = Gamma[blag]
-                    interp, HR = parabolic(Gamma, blag)
-                    if not interp:
-                        f0 = np.nan
-                        HR = np.nan
-                    else:
-                        # get fundamental frequency:
-                        f0 = samplerate / interp
+                    # get fundamental frequency:
+                    f0 = samplerate / interp
             if f0 > samplerate/2:
                 raise ValueError("F0 value ({0}) is above the nyquist rate "
                                  "({1}). This shouldn't happen...".format(f0,
@@ -224,8 +219,8 @@ class F0Analysis(Analysis):
         '''
         samplerate = self.AnalysedAudioFile.samplerate
         frames = args[0]
-        frames = multirate.interp(frames, 2)
-        samplerate *= 2
+        # frames = multirate.interp(frames, 4)
+        samplerate *= 1
         data = self.create_f0_analysis(frames, samplerate, **kwargs)
         f0 = data[:, 0]
         harmonic_ratio = data[:, 1]

src/sppysound/analysis/F0HarmRatioAnalysis.py

@@ -57,7 +57,7 @@ class F0HarmRatioAnalysis(Analysis):
     def calc_F0HarmRatio_frame_times(F0HarmRatioframes, sample_frames, samplerate):
 
         """Calculate times for frames using sample size and samplerate."""
-        samplerate *= 4
+        samplerate *= 1
 
         # Get number of frames for time and frequency
         timebins = F0HarmRatioframes.shape[0]

src/sppysound/config.py

@@ -1,30 +1,30 @@
 # Specify analysis parameters for root mean square analysis.
 rms = {
-    "window_size": 100,
+    "window_size": 70,
     "overlap": 8,
 }
 
 f0 = {
-    "window_size": 8192,
+    "window_size": 2048,
     "overlap": 8,
-    "ratio_threshold": 0.0
+    "ratio_threshold": 0.1
 }
 
 # Specify analysis parameters for variance analysis.
 variance = {
-    "window_size": 100,
+    "window_size": 70,
     "overlap": 8
 }
 
 # Specify analysis parameters for temporal kurtosis analysis.
 kurtosis = {
-    "window_size": 100,
+    "window_size": 70,
     "overlap": 8
 }
 
 # Specify analysis parameters for temporal skewness analysis.
 skewness = {
-    "window_size": 100,
+    "window_size": 70,
     "overlap": 8
 }
 
@@ -42,20 +42,20 @@ database = {
 # Sets the weighting for each analysis. a higher weighting gives an analysis
 # higher presendence when finding the best matches.
 matcher_weightings = {
-    "f0" : 2.,
-    "spccntr" : 1.,
-    "spcsprd" : 2.,
-    "spcflux" : 2.,
-    "spccf" : 2.,
-    "spcflatness": 3.,
+    "f0" : 1,
+    "spccntr" : 0.,
+    "spcsprd" : 0.,
+    "spcflux" : 0.,
+    "spccf" : 0.,
+    "spcflatness": 0.,
     "zerox" : 0.,
     "rms" : 0,
     "peak": 0.,
-    "centroid": 1.,
-    "kurtosis": 1.,
-    "skewness": 1.,
-    "variance": 2.,
-    "harm_ratio": 2.
+    "centroid": 0.,
+    "kurtosis": 0.,
+    "skewness": 0.,
+    "variance": 0.,
+    "harm_ratio": 0.
 }
 
 # Specifies the method for averaging analysis frames to create a single value
@@ -88,7 +88,7 @@ analysis = {
 matcher = {
     # Force the re-matching of analyses
     "rematch": False,
-    "grain_size": 100,
+    "grain_size": 70,
     "overlap": 8,
     # Defines the number of matches to keep for synthesis. Note that this must
     # also be specified in the synthesis config
@@ -103,13 +103,13 @@ synthesizer = {
     # between source and target.
     "enforce_intensity": True,
     # Specify the ratio limit that is the grain can be scaled by.
-    "enf_intensity_ratio_limit": 5.,
+    "enf_intensity_ratio_limit": 25.,
     # Artificially modify the pitch by the difference in f0 values between
     # source and target.
     "enforce_f0": True,
     # Specify the ratio limit that is the grain can be modified by.
-    "enf_f0_ratio_limit": 10.,
-    "grain_size": 100,
+    "enf_f0_ratio_limit": 100.,
+    "grain_size": 70,
     "overlap": 8,
     # Normalize output, avoid clipping of final output by scaling the final
     # frames.

src/sppysound/database.py

@@ -429,11 +429,11 @@ class Matcher:
                 all_target_analyses[i] = target_data
 
 
-            nan_columns = np.all(np.isnan(all_target_analyses), axis=0)
-            all_target_analyses[:, nan_columns] = 0.
+            # nan_columns = np.all(np.isnan(all_target_analyses), axis=0)
+            # all_target_analyses[:, nan_columns] = 0.
             # Impute values for Nans
-            all_target_analyses = imp.fit_transform(all_target_analyses)
-            # all_target_analyses[np.isnan(all_target_analyses)] = np.inf
+            # all_target_analyses = imp.fit_transform(all_target_analyses)
+            all_target_analyses[np.isnan(all_target_analyses)] = np.inf
 
             for sind, source_entry in enumerate(self.source_db.analysed_audio):
                 self.logger.info("K-d Tree Matching: {0} to {1}".format(source_entry.name, target_entry.name))
@@ -453,11 +453,11 @@ class Matcher:
 
 
                 # Impute values for Nans
-                nan_columns = np.all(np.isnan(all_source_analyses), axis=0)
-                all_source_analyses[:, nan_columns] = 0.
-                all_source_analyses = imp.fit_transform(all_source_analyses)
+                # nan_columns = np.all(np.isnan(all_source_analyses), axis=0)
+                # all_source_analyses[:, nan_columns] = 0.
+                # all_source_analyses = imp.fit_transform(all_source_analyses)
 
-                # all_source_analyses[np.isnan(all_source_analyses)] = np.inf
+                all_source_analyses[np.isnan(all_source_analyses)] = np.inf
 
                 source_tree = spatial.cKDTree(all_source_analyses.T, leafsize=100)
                 results_vals, results_inds = source_tree.query(all_target_analyses.T, k=self.match_quantity, p=2)

src/sppysound/multirate.py

@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+"""Module providing Multirate signal processing functionality.
+Largely based on MATLAB's Multirate signal processing toolbox with consultation
+of Octave m-file source code.
+
+Ref: https://github.com/mubeta06/python/blob/master/signal_processing/sp/multirate.py
+"""
+
+import sys
+import fractions
+import numpy
+from scipy import signal
+
+
+def downsample(s, n, phase=0):
+    """Decrease sampling rate by integer factor n with included offset phase.
+    """
+    return s[phase::n]
+
+
+def upsample(s, n, phase=0):
+    """Increase sampling rate by integer factor n  with included offset phase.
+    """
+    return numpy.roll(numpy.kron(s, numpy.r_[1, numpy.zeros(n-1)]), phase)
+
+
+def decimate(s, r, n=None, fir=False):
+    """Decimation - decrease sampling rate by r. The decimation process filters
+    the input data s with an order n lowpass filter and then resamples the
+    resulting smoothed signal at a lower rate. By default, decimate employs an
+    eighth-order lowpass Chebyshev Type I filter with a cutoff frequency of
+    0.8/r. It filters the input sequence in both the forward and reverse
+    directions to remove all phase distortion, effectively doubling the filter
+    order. If 'fir' is set to True decimate uses an order 30 FIR filter (by
+    default otherwise n), instead of the Chebyshev IIR filter. Here decimate
+    filters the input sequence in only one direction. This technique conserves
+    memory and is useful for working with long sequences.
+    """
+    if fir:
+        if n is None:
+            n = 30
+        b = signal.firwin(n, 1.0/r)
+        a = 1
+        f = signal.lfilter(b, a, s)
+    else: #iir
+        if n is None:
+            n = 8
+        b, a = signal.cheby1(n, 0.05, 0.8/r)
+        f = signal.filtfilt(b, a, s)
+    return downsample(f, r)
+
+
+def interp(s, r, l=4, alpha=0.5):
+    """Interpolation - increase sampling rate by integer factor r. Interpolation
+    increases the original sampling rate for a sequence to a higher rate. interp
+    performs lowpass interpolation by inserting zeros into the original sequence
+    and then applying a special lowpass filter. l specifies the filter length
+    and alpha the cut-off frequency. The length of the FIR lowpass interpolating
+    filter is 2*l*r+1. The number of original sample values used for
+    interpolation is 2*l. Ordinarily, l should be less than or equal to 10. The
+    original signal is assumed to be band limited with normalized cutoff
+    frequency 0=alpha=1, where 1 is half the original sampling frequency (the
+    Nyquist frequency). The default value for l is 4 and the default value for
+    alpha is 0.5.
+    """
+    b = signal.firwin(2*l*r+1, alpha/r);
+    a = 1
+    return r*signal.lfilter(b, a, upsample(s, r))[r*l+1:-1]
+
+
+def resample(s, p, q, h=None):
+    """Change sampling rate by rational factor. This implementation is based on
+    the Octave implementation of the resample function. It designs the
+    anti-aliasing filter using the window approach applying a Kaiser window with
+    the beta term calculated as specified by [2].
+
+    Ref [1] J. G. Proakis and D. G. Manolakis,
+    Digital Signal Processing: Principles, Algorithms, and Applications,
+    4th ed., Prentice Hall, 2007. Chap. 6
+    Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
+    Discrete-time signal processing, Signal processing series,
+    Prentice-Hall, 1999
+    """
+    gcd = fractions.gcd(p,q)
+    if gcd>1:
+        p=p/gcd
+        q=q/gcd
+
+    if h is None: #design filter
+        #properties of the antialiasing filter
+        log10_rejection = -3.0
+        stopband_cutoff_f = 1.0/(2.0 * max(p,q))
+        roll_off_width = stopband_cutoff_f / 10.0
+
+        #determine filter length
+        #use empirical formula from [2] Chap 7, Eq. (7.63) p 476
+        rejection_db = -20.0*log10_rejection;
+        l = numpy.ceil((rejection_db-8.0) / (28.714 * roll_off_width))
+
+        #ideal sinc filter
+        t = numpy.arange(-l, l + 1)
+        ideal_filter=2*p*stopband_cutoff_f*numpy.sinc(2*stopband_cutoff_f*t)
+
+        #determine parameter of Kaiser window
+        #use empirical formula from [2] Chap 7, Eq. (7.62) p 474
+        beta = signal.kaiser_beta(rejection_db)
+
+        #apodize ideal filter response
+        h = numpy.kaiser(2*l+1, beta)*ideal_filter
+
+    ls = len(s)
+    lh = len(h)
+
+    l = (lh - 1)/2.0
+    ly = numpy.ceil(ls*p/float(q))
+
+    #pre and postpad filter response
+    nz_pre = numpy.floor(q - numpy.mod(l,q))
+    hpad = h[-lh+nz_pre:]
+
+    offset = numpy.floor((l+nz_pre)/q)
+    nz_post = 0;
+    while numpy.ceil(((ls-1)*p + nz_pre + lh + nz_post )/q ) - offset < ly:
+        nz_post += 1
+    hpad = hpad[:lh + nz_pre + nz_post]
+
+    #filtering
+    xfilt = upfirdn(s, hpad, p, q)
+
+    return xfilt[offset-1:offset-1+ly]
+
+
+def upfirdn(s, h, p, q):
+    """Upsample signal s by p, apply FIR filter as specified by h, and
+    downsample by q. Using fftconvolve as opposed to lfilter as it does not seem
+    to do a full convolution operation (and its much faster than convolve).
+    """
+    return downsample(signal.fftconvolve(h, upsample(s, p)), q)
+
+def main():
+    """Show simple use cases for functionality provided by this module. Each
+    example below attempts to mimic the examples provided by mathworks MATLAB
+    documentation, http://www.mathworks.com/help/toolbox/signal/
+    """
+    import pylab
+    argv = sys.argv
+    if len(argv) != 1:
+        print >>sys.stderr, 'usage: python -m pim.sp.multirate'
+        sys.exit(2)
+
+    #Downsample
+    x = numpy.arange(1, 11)
+    print 'Down Sampling %s by 3' % x
+    print  downsample(x, 3)
+    print 'Down Sampling %s by 3 with phase offset 2' % x
+    print  downsample(x, 3, phase=2)
+
+    #Upsample
+    x = numpy.arange(1, 5)
+    print 'Up Sampling %s by 3' % x
+    print upsample(x, 3)
+    print 'Up Sampling %s by 3 with phase offset 2' % x
+    print upsample(x, 3, 2)
+
+    #Decimate
+    t = numpy.arange(0, 1, 0.00025)
+    x = numpy.sin(2*numpy.pi*30*t) + numpy.sin(2*numpy.pi*60*t)
+    y = decimate(x,4)
+    pylab.figure()
+    pylab.subplot(2, 1, 1)
+    pylab.title('Original Signal')
+    pylab.stem(numpy.arange(len(x[0:120])), x[0:120])
+    pylab.subplot(2, 1, 2)
+    pylab.title('Decimated Signal')
+    pylab.stem(numpy.arange(len(y[0:30])), y[0:30])
+
+    #Interp
+    t = numpy.arange(0, 1, 0.001)
+    x = numpy.sin(2*numpy.pi*30*t) + numpy.sin(2*numpy.pi*60*t)
+    y = interp(x,4)
+    pylab.figure()
+    pylab.subplot(2, 1, 1)
+    pylab.title('Original Signal')
+    pylab.stem(numpy.arange(len(x[0:30])), x[0:30])
+    pylab.subplot(2, 1, 2)
+    pylab.title('Interpolated Signal')
+    pylab.stem(numpy.arange(len(y[0:120])), y[0:120])
+
+    #upfirdn
+    L = 147.0
+    M = 160.0
+    N = 24.0*L
+    h = signal.firwin(N-1, 1/M, window=('kaiser', 7.8562))
+    h = L*h
+    Fs = 48000.0
+    n = numpy.arange(0, 10239)
+    x  = numpy.sin(2*numpy.pi*1000/Fs*n)
+    y = upfirdn(x, h, L, M)
+    pylab.figure()
+    pylab.stem(n[1:49]/Fs, x[1:49])
+    pylab.stem(n[1:45]/(Fs*L/M), y[13:57], 'r', markerfmt='ro',)
+    pylab.xlabel('Time (sec)')
+    pylab.ylabel('Signal value')
+
+    #resample
+    fs1 = 10.0
+    t1 = numpy.arange(0, 1 + 1.0/fs1, 1.0/fs1)
+    x = t1
+    y = resample(x, 3, 2)
+    t2 = numpy.arange(0,(len(y)))*2.0/(3.0*fs1)
+    pylab.figure()
+    pylab.plot(t1, x, '*')
+    pylab.plot(t2, y, 'o')
+    pylab.plot(numpy.arange(-0.5,1.5, 0.01), numpy.arange(-0.5,1.5, 0.01), ':')
+    pylab.legend(('original','resampled'))
+    pylab.xlabel('Time')
+
+    x = numpy.hstack([numpy.arange(1,11), numpy.arange(9,0,-1)])
+    y = resample(x,3,2)
+    pylab.figure()
+    pylab.subplot(2, 1, 1)
+    pylab.title('Edge Effects Not Noticeable')
+    pylab.plot(numpy.arange(19)+1, x, '*')
+    pylab.plot(numpy.arange(29)*2/3.0 + 1, y, 'o')
+    pylab.legend(('original', 'resampled'))
+    x = numpy.hstack([numpy.arange(10, 0, -1), numpy.arange(2,11)])
+    y = resample(x,3,2)
+    pylab.subplot(2, 1, 2)
+    pylab.plot(numpy.arange(19)+1, x, '*')
+    pylab.plot(numpy.arange(29)*2/3.0 + 1, y, 'o')
+    pylab.title('Edge Effects Very Noticeable')
+    pylab.legend(('original', 'resampled'))
+
+    pylab.show()
+    return 0
+
+if __name__ == '__main__':
+    sys.exit(main())