### suppression git frequency estimator, récupération de ses .py importants

parent f92492b2
frequency_estimator @ f94aa1d3
 Subproject commit f94aa1d3d63a8ec3c0e80f30c6f9f4e375ceea0f
 from __future__ import division from numpy.fft import rfft from numpy import argmax, mean, diff, log, nonzero from scipy.signal import blackmanharris, correlate from time import time import sys try: import soundfile as sf except ImportError: from scikits.audiolab import flacread from frequency_estimator.parabolic import parabolic def freq_from_crossings(sig, fs): """ Estimate frequency by counting zero crossings """ # Find all indices right before a rising-edge zero crossing indices = nonzero((sig[1:] >= 0) & (sig[:-1] < 0)) # Naive (Measures 1000.185 Hz for 1000 Hz, for instance) # crossings = indices # More accurate, using linear interpolation to find intersample # zero-crossings (Measures 1000.000129 Hz for 1000 Hz, for instance) crossings = [i - sig[i] / (sig[i+1] - sig[i]) for i in indices] # Some other interpolation based on neighboring points might be better. # Spline, cubic, whatever return fs / mean(diff(crossings)) def freq_from_fft(sig, fs): """ Estimate frequency from peak of FFT """ # Compute Fourier transform of windowed signal windowed = sig * blackmanharris(len(sig)) f = rfft(windowed) # Find the peak and interpolate to get a more accurate peak i = argmax(abs(f)) # Just use this for less-accurate, naive version true_i = parabolic(log(abs(f)), i) # Convert to equivalent frequency return fs * true_i / len(windowed) def freq_from_autocorr(sig, fs): """ Estimate frequency using autocorrelation """ # Calculate autocorrelation and throw away the negative lags corr = correlate(sig, sig, mode='full') corr = corr[len(corr)//2:] # Find the first low point d = diff(corr) start = nonzero(d > 0) # Find the next peak after the low point (other than 0 lag). This bit is # not reliable for long signals, due to the desired peak occurring between # samples, and other peaks appearing higher. # Should use a weighting function to de-emphasize the peaks at longer lags. peak = argmax(corr[start:]) + start px, py = parabolic(corr, peak) return fs / px def freq_from_HPS(sig, fs): """ Estimate frequency using harmonic product spectrum (HPS) """ windowed = sig * blackmanharris(len(sig)) from pylab import subplot, plot, log, copy, show import matplotlib.pyplot as plt # harmonic product spectrum: c = abs(rfft(windowed)) maxharms = 8 resarray = [] for x in range(2, 6): a = copy(c[::x]) # Should average or maximum instead of decimating # max(c[::x],c[1::x],c[2::x],...) c = c[:len(a)] i = argmax(abs(c)) true_i = parabolic(abs(c), i) res = fs * true_i / len(windowed) resarray.append(res) #print('Pass %d: %f Hz' % (x, res)) c *= a #plt.subplot(maxharms, 1, x) #plt.plot(log(c)) show() return mean(resarray[1:3]) def hello(signal, fs): print('Calculating frequency from FFT:', end=' ') start_time = time() print('%f Hz' % freq_from_fft(signal, fs)) print('Time elapsed: %.3f s\n' % (time() - start_time)) print('Calculating frequency from zero crossings:', end=' ') start_time = time() print('%f Hz' % freq_from_crossings(signal, fs)) print('Time elapsed: %.3f s\n' % (time() - start_time)) print('Calculating frequency from autocorrelation:', end=' ') start_time = time() print('%f Hz' % freq_from_autocorr(signal, fs)) print('Time elapsed: %.3f s\n' % (time() - start_time)) # print('Calculating frequency from harmonic product spectrum:') # start_time = time() # freq_from_HPS(signal, fs) # print('Time elapsed: %.3f s\n' % (time() - start_time))
parabolic.py 0 → 100644
 from __future__ import division from numpy import polyfit, arange def parabolic(f, x): """Quadratic interpolation for estimating the true position of an inter-sample maximum when nearby samples are known. f is a vector and x is an index for that vector. Returns (vx, vy), the coordinates of the vertex of a parabola that goes through point x and its two neighbors. Example: Defining a vector f with a local maximum at index 3 (= 6), find local maximum if points 2, 3, and 4 actually defined a parabola. In : f = [2, 3, 1, 6, 4, 2, 3, 1] In : parabolic(f, argmax(f)) Out: (3.2142857142857144, 6.1607142857142856) """ xv = 1/2. * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x yv = f[x] - 1/4. * (f[x-1] - f[x+1]) * (xv - x) return (xv, yv) def parabolic_polyfit(f, x, n): """Use the built-in polyfit() function to find the peak of a parabola f is a vector and x is an index for that vector. n is the number of samples of the curve used to fit the parabola. """ a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2) xv = -0.5 * b/a yv = a * xv**2 + b * xv + c return (xv, yv) if __name__ == "__main__": from numpy import argmax import matplotlib.pyplot as plt y = [2, 1, 4, 8, 11, 10, 7, 3, 1, 1] xm, ym = argmax(y), y[argmax(y)] xp, yp = parabolic(y, argmax(y)) plot = plt.plot(y) plt.plot(xm, ym, 'o', color='silver') plt.plot(xp, yp, 'o', color='blue') plt.title('silver = max, blue = estimated max')