frequency_estimator.py 2.65 KB
 Candre Nathan committed Apr 02, 2021 1 2 3 4 5 ``````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 `````` LucasDANIELE committed Apr 04, 2021 6 `````` `````` Candre Nathan committed Apr 02, 2021 7 `````` `````` LucasDANIELE committed Apr 02, 2021 8 ``````from parabolic import parabolic `````` Candre Nathan committed Apr 02, 2021 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 `````` 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))[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)[0] # 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)[0][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)) `````` LucasDANIELE committed Apr 03, 2021 80 `````` maxharms = 6 `````` Candre Nathan committed Apr 02, 2021 81 `````` resarray = [] `````` LucasDANIELE committed Apr 03, 2021 82 `````` for x in range(2, maxharms): `````` Candre Nathan committed Apr 02, 2021 83 84 85 86 87 88 89 90 `````` a = copy(c[::x]) # Should average or maximum instead of decimating c = c[:len(a)] i = argmax(abs(c)) true_i = parabolic(abs(c), i)[0] res = fs * true_i / len(windowed) resarray.append(res) c *= a show() `````` LucasDANIELE committed Apr 05, 2021 91 `` return mean(resarray[1:3])``