Katt
Katt

Reputation: 1017

frequency analysis with unevenly spaced data in python

I have a signal generated by a simulation program. Because the solver in this program has a variable time step, I have a signal with unevenly spaced data. I have two lists, a list with the signal values, and another list with the times at which each value occurred. The data could be something like this

npts = 500
t=logspace(0,1,npts)
f1 = 0.5 
f2 = 0.6
sig=(1+sin(2*pi*f1*t))+(1+sin(2*pi*f2*t))

I would like to be able to perform a frequency analysis on this signal using python. It seems I cannot use the fft function in numpy, because this requires evenly spaced data. Are there any standard functions which could help me find the frequencies contained in this signal?

Upvotes: 6

Views: 6745

Answers (3)

kilojoules
kilojoules

Reputation: 10083

An easy way out is to interpolate to evenly-spaced time intervals

Upvotes: 0

DrM
DrM

Reputation: 2525

Very simple, just look up the formula for a Fourier transform, and implement it as a discrete sum over your data values:

given a set of values f(x) over some set of x, then for each frequency k,

F(k) = sum_x ( exp( +/-i * k *x ) )

choose your k's ranging from 0 to 2*pi / min separation in x.

and, you can use 2 * pi / max(x) as the increment size

For a test case, use something for which you known the correct answer, c.f., a single cos( k' * x ) for some k', or a Gaussian.

Upvotes: 2

PearsonArtPhoto
PearsonArtPhoto

Reputation: 39698

The most common algorithm to solve such things is called a Least-Squares Spectral analysis of frequencies. It looks like this will be in a future release of the scipy.signals package. Maybe there is a current version, but I can't seem to find it... In addition, there is some code available from Astropython, which I will not copy in it's entirety, but it essentially creates a lomb class which you can use the following code to get some values out.. What you need to do is the following:

import numpy
import lomb
x = numpy.arange(10)
y = numpy.sin(x)
fx,fy, nout, jmax, prob = lomb.fasper(x,y, 6., 6.)

Upvotes: 7

Related Questions