How to calculate frequency of a give wave and time

Question

I have data for Velocity vs time. The time steps are not uniform, but the Velocity data is a wave. How do I calculate the principal frequency of the velocity using FFT of Python? Most of the examples I am seeing online are for uniform time stepping.

My data is like

7.56683038E+02  2.12072850E-01 
7.56703750E+02  2.13280844E-01
7.56724461E+02  2.14506402E-01
7.56745172E+02  2.15748934E-01
7.56765884E+02  2.17007907E-01
7.56786595E+02  2.18282753E-01

10000 lines like that.

Seeing some online responses, I wrote a code like the following, but it is not working:

#!/usr/bin/env python

import numpy as np
import scipy as sy
import scipy.fftpack as syfp
import pylab as pyl

# Calculate the number of data points
length = 0
for line in open("data.dat"):
    length = length + 1

# Read in data from file here
t = np.zeros(shape=(length,1))
u = np.zeros(shape=(length,1))


length = 0
for line in open("data.dat"):
    columns = line.split(' ')
    t[length] = float(columns[0])
    u[length] = float(columns[1])
    length = length + 1

# Do FFT analysis of array
FFT = sy.fft(u)

# Getting the related frequencies
freqs = syfp.fftfreq(len(u))

# Create subplot windows and show plot
pyl.subplot(211)
pyl.plot(t, u)
pyl.xlabel('Time')
pyl.ylabel('Amplitude')
pyl.subplot(212)
pyl.plot(freqs, sy.log10(FFT), 'x')
pyl.show()

---------------------- edit ------------------------

with this code I am getting an output like the following figure. I am not sure what this figure shows. I was expecting just to see one peak in the FFT diagram

---------------------- edit ------------------------

My results with the mock data with the sin functions suggested in the comments below are here:

Demicanton

What do you mean by "it is not working"? Also, what do you mean by "principal frequency"? – Ephemeron

In the FFT plot, what do you see if you plot sy.log10(np.abs(FFT)) (ie, note the use of the abs there). – Ephemeron

that does not change anything in the output. btw, I checked and the time steps look almost uniform with dt= 0.02071 and I have 100000 data points in my file. – Demicanton

It's hard to diagnose the result with data that we don't have. How about trying with some fake data, where you know what the result should look like: t = 0.02071*np.arange(100000) and u = np.sin(2*np.pi*t*.01)? – Ephemeron

thanks for the idea. I will try that. Meanwhile, I have made the data available at drive.google.com/file/d/0B9w1Yb4SFBxrWEFGTERVX1dZSGM/…. The first coloumn is time, and the second is velocity. Don't bother with the other coloumns. – Demicanton

Personally, I'm going to wait for you to try with the mock data. It will help you, and it's easier for me. – Ephemeron

@Ephemeron I have added your mock data results at the top. There also I don't see the frequency. I feel like I am missing something fundamental here :) – Demicanton

Answer 1

From what I can see, your code is basically fine, but missing a few details. I think your issues are mostly about interpretation. Because of this, the mock data is the best to look at now, and here's an example with the mock data I suggested in the comments (and I've added comments about the important lines, and ## for changes):

import numpy as np
import scipy as sy
import scipy.fftpack as syfp
import pylab as pyl

dt = 0.02071 
t = dt*np.arange(100000)             ## time at the same spacing of your data
u = np.sin(2*np.pi*t*.01)            ## Note freq=0.01 Hz

# Do FFT analysis of array
FFT = sy.fft(u)

# Getting the related frequencies
freqs = syfp.fftfreq(len(u), dt)     ## added dt, so x-axis is in meaningful units

# Create subplot windows and show plot
pyl.subplot(211)
pyl.plot(t, u)
pyl.xlabel('Time')
pyl.ylabel('Amplitude')
pyl.subplot(212)
pyl.plot(freqs, sy.log10(abs(FFT)), '.')  ## it's important to have the abs here
pyl.xlim(-.05, .05)                       ## zoom into see what's going on at the peak
pyl.show()

As you can see, there are two peaks, at + and - the input frequency (.01 Hz), as expected.

Edit:
Puzzled why this approach didn't work for the OP's data, I took a look at that too. The problem is that the sample times aren't evenly spaced. Here's a histogram of the times (code below).

So the time between samples is roughly evenly split between a short time and a long time. I took a quick look for a pattern here and nothing was obvious.

To do an FFT, one needs evenly spaced time samples, so I interpolated to get the following:

which is reasonable (a DC offset, a primary peak and a small harmonic). Here's the code:

data = np.loadtxt("data.dat", usecols=(0,1))
t = data[:,0]
u = data[:,1]

tdiff = t[1:]-t[:-1]
plt.hist(tdiff)

new_times = np.linspace(min(t), max(t), len(t))
new_data = np.interp(new_times, t, u)

t, u = new_times, new_data

FFT = sy.fft(u)
freqs = syfp.fftfreq(len(u), dt)

# then plot as above

Ephemeron

thanks for your reply. It does look like the code is doing the right thing. I made the changes you recommended. But still, after I am inputting my data, the output of the FFT I am getting is like a wave (as shown in the first figure of my question). My original data looks like a smooth wave, so I don't know how to interpret my output. – Demicanton

@jhaprade: I ended up being puzzled by this, so I did it and included the answer to your data issue in the edit. – Ephemeron