endolith/frequency_estimator.py

## readme.md

      
    Raw
  

              readme.md
            
          
    A few simple frequency estimation methods in Python
See also https://github.com/endolith/waveform-analyzer/blob/master/frequency_estimator.py, which is mostly the same thing, maybe more up-to-date.  I need to keep them both in sync with each other or delete one.
These are the methods that everyone recommends when someone asks about
frequency estimation or pitch detection.  Such as here:

Music - How do you analyse the fundamental frequency of a PCM or WAV sample
CCRMA Pitch detection methods review
Pitch Detection Algorithms (Middleton)

So these are my attempts at implementation.  Initially I was trying to measure the frequency of long sine waves with high accuracy (to indirectly measure clock frequency), then added methods for other types of signals later.
None of them work well in all situations, these are "offline", not real-time, and I am sure there are much better methods "in the literature", but here is some sample code for the simple methods at least.
Count zero-crossings, divide average period by time to get frequency


Works well for long low-noise sines, square, triangle, etc.
Supposedly this is how cheap guitar tuners work
Using interpolation to find a "truer" zero-crossing gives better accuracy
Pro: Fast
Pro: Accurate (increasing with signal length)
Con: Doesn't work if there are multiple zero crossings per cycle, low-frequency baseline shift, noise, etc.

Do FFT and find the peak


Using parabolic interpolation to find a truer peak gives better accuracy
Accuracy also increases with signal/FFT length
Con: Doesn't find the right value if harmonics are stronger than fundamental, which is common.  Better method would try to be smarter about identifying the fundamental, like template matching using the "two-way mismatch" (TWM) algorithm.
Pro: Accurate, usually even more so than zero crossing counter (1000.000004 Hz for 1000 Hz, for instance).  Due to parabolic interpolation being a very good fit for windowed log FFT peaks?

Do autocorrelation and find the peak


Pro: Best method for finding the true fundamental of any repetitive wave, even with weak or missing fundamental (finds GCD of all harmonics present)
Con: Inaccurate result if waveform isn't perfectly repeating, like inharmonic musical instruments (piano, guitar, ...), however:
Pro: This inaccurate result more closely matches the pitch that humans perceive :)
Con: Not as accurate as other methods for precise measurement of sine waves
Con: This implementation has trouble with finding the true peak

Calculate harmonic product spectrum and find the peak


Pro: Good at finding the true fundamental even if weak or missing


## frequency_estimator.py
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 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

    # harmonic product spectrum:
    c = abs(rfft(windowed))
    maxharms = 8
    subplot(maxharms, 1, 1)
    plot(log(c))
    for x in range(2, maxharms):
        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)[0]
        print('Pass %d: %f Hz' % (x, fs * true_i / len(windowed)))
        c *= a
        subplot(maxharms, 1, x)
        plot(log(c))
    show()


filename = sys.argv[1]

print('Reading file "%s"\n' % filename)
try:
    signal, fs = sf.read(filename)
except NameError:
    signal, fs, enc = flacread(filename)

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.md

      
    Raw
  

              parabolic.md
            
          
    Motivation:

Note that the Gaussian window transform magnitude is precisely a parabola on a dB scale. As a result, quadratic spectral peak interpolation is exact under the Gaussian window. Of course, we must somehow remove the infinitely long tails of the Gaussian window in practice, but this does not cause much deviation from a parabola


Quadratic Interpolation of Spectral Peaks

Apparently this produces error if the peak is not exactly at the edge or center of a bin.  Interpolating by zero-padding before the FFT does not produce this kind of error, but is more computationally expensive.  So a good trade-off is to do some zero-padding interpolation and then follow with parabolic interpolation:

While we could choose our zero-padding factor large enough to yield any desired degree of accuracy in peak frequency measurements, it is more efficient in practice to combine zero-padding with parabolic interpolation (or some other simple, low-order interpolator). In such hybrid schemes, the zero-padding is simply chosen large enough so that the bias due to parabolic interpolation is negligible.


Bias of Parabolic Peak Interpolation 

Here's a Matlab function that does the same thing.
Actually this function could be done with the polyfit() function, instead, but it requires more steps:
f = [2, 3, 1, 6, 4, 2, 3, 1]

In [8]: parabolic(f, argmax(f))
Out[8]: (3.2142857142857144, 6.1607142857142856)

In [9]: a, b, c = polyfit([2,3,4],[1,6,4],2)
In [10]: x = -0.5*b/a

In [11]: x
Out[11]: 3.2142857142857695

In [12]: a*x**2 + b*x + c
Out[12]: 6.1607142857143025

parabolic() could be updated to use this, and then it will work for non-uniform sampled input as well.
Alternative estimators with even lower error: How to Interpolate the Peak Location of a DFT or FFT if the Frequency of Interest is Between Bins

  
## parabolic.py
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 [3]: f = [2, 3, 1, 6, 4, 2, 3, 1]

    In [4]: parabolic(f, argmax(f))
    Out[4]: (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')
	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 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

	# harmonic product spectrum:
	c = abs(rfft(windowed))
	maxharms = 8
	subplot(maxharms, 1, 1)
	plot(log(c))
	for x in range(2, maxharms):
	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)[0]
	print('Pass %d: %f Hz' % (x, fs * true_i / len(windowed)))
	c *= a
	subplot(maxharms, 1, x)
	plot(log(c))
	show()


	filename = sys.argv[1]

	print('Reading file "%s"\n' % filename)
	try:
	signal, fs = sf.read(filename)
	except NameError:
	signal, fs, enc = flacread(filename)

	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))
	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 [3]: f = [2, 3, 1, 6, 4, 2, 3, 1]

	In [4]: parabolic(f, argmax(f))
	Out[4]: (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')