(to be determined) home

Estimate frequency using Numpy

Estimating frequency of noisy sine curve

Numpy/Scipy doesn’t have function for estimating frequency (pitch) of given 1D time series, but simple functions can do this. I use a very simple method called “Pisarenko Harmonic Decomposition”. On the picture above, I estimate the frequency of sine curve + noise. The estimate (vertical black line on the top panel) is matched to maximum of the power spectral density. Note that argument x of function freq must satisfy x.mean() == 0.

import numpy
PI = numpy.pi

def covariance(x, k):
    N = len(x) - k
    return (x[:-k] * x[k:]).sum() / N

def phd1(x):
    """Estimate frequency using Pisarenko Harmonic Decomposition"""
    r1 = covariance(x, 1)
    r2 = covariance(x, 2)
    a = (r2 + numpy.sqrt(r2 ** 2 + 8 * r1 ** 2)) / 4 / r1
    if a > 1:
        a = 1
    elif a < -1:
        a = -1
    return numpy.arccos(a)

def freq(x, sample_step=1, dt=1.0):
    """Estimate frequency using `phd1`"""
    omega = phd1(x[::sample_step])
    return omega / 2.0 / PI / sample_step / dt

I posted full source code to get the picture above here: gist: 565034 - GitHub.

You can find the formulation of Pisarenko Harmonic Decomposition method here: CiteSeerX — AN UNBIASED PISARENKO HARMONIC DECOMPOSITION ESTIMATOR FOR SINGLE-TONE FREQUENCY.

blog comments powered by Disqus