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
from 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))[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))
maxharms = 6
resarray = []
for x in range(2, maxharms):
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()
return mean(resarray[1:3])