Merge pull request #375 from tlecomte/fixFftSmooth

fix(Spectrum): copy data in harmonic product estimation

Author: tlecomte

friture/spectrum.py

@@ -18,7 +18,7 @@
 # along with Friture.  If not, see <http://www.gnu.org/licenses/>.
 
 from PyQt5.QtCore import QObject
-from numpy import log10, argmax, zeros, arange, floor, float64
+from numpy import log10, argmax, zeros, arange, floor, float64, ones
 from friture.audioproc import audioproc  # audio processing class
 from friture.spectrum_settings import (Spectrum_Settings_Dialog,  # settings dialog
                                        DEFAULT_FFT_SIZE,
@@ -105,11 +105,21 @@ def harmonic_product_spectrum(self, sp):
         # for pitch detection are good.
         product_count = 3
 
-        # Downsample and multiply
         harmonic_length = sp.shape[0] // product_count
-        res = sp[:harmonic_length]
-        for i in range(2, product_count + 1):
-            res *= sp[::i][:harmonic_length]
+
+        # common case: product_count == 3 — multiply the three slices
+        # directly. This avoids allocating a temporary ones array and still
+        # doesn't make copies of the original `sp` beyond NumPy's routine
+        # temporary results for the elementwise multiplication.
+        if product_count == 3:
+            # Downsample and multiply
+            res = (sp[:harmonic_length]
+                   * sp[::2][:harmonic_length]
+                   * sp[::3][:harmonic_length])
+        else:
+            res = ones(harmonic_length, dtype=sp.dtype)
+            for i in range(1, product_count + 1):
+                res *= sp[::i][:harmonic_length]
         return res
 
     def handle_new_data(self, floatdata):
@@ -191,7 +201,14 @@ def setresponsetime(self, response_time):
         self.response_time = response_time
 
         # an exponential smoothing filter is a simple IIR filter
-        # s_i = alpha*x_i + (1-alpha)*s_{i-1}
+        # y_i     = a*x_i     + (1-a)*y_{i-1}
+        #
+        # we can unroll the recurrence a bit:
+        # y_{i+1} = a*x_{i+1} + (1-a)*y_i
+        #         = a*x_{i+1} + (1-a)*a*x_i + (1-a)^2*y_{i-1}
+        # y_{i+2} = a*x_{i+2} + (1-a)*y_{i+1}
+        #         = a*x_{i+2} + (1-a)*a*x_{i+1} + (1-a)^2*a*x_i + (1-a)^3*y_{i-1}
+        # ...
         # we compute alpha so that the N most recent samples represent 100*w percent of the output
         w = 0.65
         delta_n = self.fft_size * (1. - self.overlap)

friture/spectrum_settings.py

@@ -190,7 +190,7 @@ def responsetimechanged(self, index):
             response_time = 1.
         elif index == 4:
             response_time = 5.            
-        self.logger.info("responsetimechanged slot %d %d", index, response_time)
+        self.logger.info("responsetimechanged slot %d %f", index, response_time)
         self.view_model.setresponsetime(response_time)
 
     # method

friture_extensions/exp_smoothing_conv.pyx

@@ -15,48 +15,85 @@ def pyx_exp_smoothed_value(np.ndarray[dtype_t, ndim=1] kernel, dtype_t alpha, np
 	cdef Py_ssize_t Nk = kernel.shape[0]
 	cdef Py_ssize_t i
 	cdef dtype_t value, conv
+	cdef double a = (1.-alpha)**N
 
 	# as we disable the cython bounds checking, do it by ourselves
 	# it is absolutely needed as the kernel is not of infinite size !!
 	if N > Nk:
 		N = Nk
+		a = 0.
 
 	if N == 0:
 		value = previous
-	else:
-		conv = 0.
-		
-		for i in range(0, N):
-			conv = conv + kernel[Nk - N + i]*data[i]
-		
-		value = alpha * conv + previous*(1.-alpha)**N
+
+	conv = 0.
+	
+	for i in range(0, N):
+		conv = conv + kernel[Nk - N + i]*data[i]
+	
+	value = alpha * conv + previous*a
 
 	return value
 
+
 def pyx_exp_smoothed_value_numpy(np.ndarray[dtype_t, ndim=1] kernel, dtype_t alpha, np.ndarray[dtype_t, ndim=2] data, np.ndarray[dtype_t, ndim=1] previous):
-	cdef Py_ssize_t N = data.shape[1]
+	"""
+	Compute exponential filtering of data with given kernel and alpha.
+	The filter is applied on the 2nd axis (axis=1), independently for each row (axis=0).
+
+	This is equivalent to applying the following IIR filter on each row of data:
+	y_i = alpha*x_i + (1-alpha)*y_{i-1}
+	where x_i is the input data, y_i the filtered data, and y_{i-1} the previous filtered value.
+
+	We can unroll the recurrence a bit:
+    y_{i+1} = a*x_{i+1} + (1-a)*y_i
+            = a*x_{i+1} + (1-a)*a*x_i + (1-a)^2*y_{i-1}
+    y_{i+2} = a*x_{i+2} + (1-a)*y_{i+1}
+            = a*x_{i+2} + (1-a)*a*x_{i+1} + (1-a)^2*a*x_i + (1-a)^3*y_{i-1}
+    ...
+
+	This unrolling can be implemented with a convolution of a precomputed kernel.         
+	
+	By using a precomputed kernel, this function is optimized to process a large number of samples at once.
+
+	Parameters:
+	-----------
+	kernel : 1D array
+		the pre-computed convolution kernel for the exponential smoothing
+	alpha : float
+		the exponential smoothing factor, between 0 and 1
+	data : 2D array
+		the input data to filter
+	previous : 1D array
+		the previous filtered value
+	Returns:
+	--------
+	1D array
+		the filtered value"""
 	cdef Py_ssize_t Nf = data.shape[0]
+	cdef Py_ssize_t Nt = data.shape[1]
 	cdef Py_ssize_t Nk = kernel.shape[0]
 	cdef Py_ssize_t i, j
 	cdef np.ndarray[dtype_t, ndim=1] value, conv
-	cdef double a = (1.-alpha)**N
+	cdef double a = (1.-alpha)**Nt
 
 	# as we disable the cython bounds checking, do it by ourselves
 	# it is absolutely needed as the kernel is not of infinite size !!
-	if N > Nk:
-		N = Nk
+	if Nt > Nk:
+		Nt = Nk
+		a = 0.
 
-	if N == 0:
+	if Nt == 0:
 		return previous
-	else:
-		conv = np.zeros(Nf)
-		value = np.zeros(Nf)
-		
-		for i in range(0, N):
-			for j in range(Nf):
-				conv[j] = conv[j] + kernel[Nk - N + i]*data[j, i]
-		
+
+	conv = np.zeros(Nf, dtype=dtype)
+	value = np.zeros(Nf, dtype=dtype)
+	
+	for i in range(0, Nt):
 		for j in range(Nf):
-			value[j] = alpha * conv[j] + previous[j]*a
+			conv[j] = conv[j] + kernel[Nk - Nt + i]*data[j, i]
 
-		return value
+	for j in range(Nf):
+		value[j] = alpha * conv[j] + previous[j]*a
+
+	return value