Frequency resolution issue using FFT in numpy

Question

I use Tektronix oscilloscope to perform some signal acquisition. I get 10.000 measurement points (few signal periods) and I have to do a frequency analysis on that set of data. My signal is 8MHz sine wave. When I use either SciPy or NumPy I get the same result - frequencies are spreaded too wide. The distance between two values is 500kHz and the highest frequency is 2.5GHz (absurd). When I want to measure frequency bandwidth around 8MHz I can only get exact values of 7.5, 8.0 and 8.5 MHz. I tried to change sample spacing determined by (x[1]-x[0]) and I got nothing better.

def CalculateFFT(t_val,p_val):
    x = t_val #Two parameters: [x,y] values
    y = lambda x: p_val
    com_signal = y(x) # Combined signal
    FFT_val = abs(scipy.fft(com_signal))
    freq_val = scipy.fftpack.fftfreq(len(com_signal), x[1]-x[0])
    spec_val = 20*scipy.log10(FFT_val)
    return freq_val, spec_val

Your measurement period should be much longer than a few signal periods to get more accurate frequency bins, what is the sampling frequency of the measurement? — Jan Kuiken, Jul 15 '13 at 17:17
No offense, did you make sure that you exactly understand how the optimal input for a DFFT must be designed and how the output of a DFFT must be interpreted (both is not trivial)? Maybe you might want to read the FFT section of my thesis here: http://gehrcke.de/files/stud/gehrcke_MScThesis_magnetic_particle_imaging.pdf (spontaneously I also do not know where the problem is) — Dr. Jan-Philip Gehrcke, Jul 15 '13 at 17:19
Thank you Jan-Philip Gehrcke, it's helpful (and nice thesis topic as well). I've done some extra simulations and I managed to see that the more signal periods I set in constant set of data (10k points) by changing the time window, the more accurate frequency values I get. — akson128, Jul 15 '13 at 20:52

score 5 · Answer 1 · answered Jul 15 '13 at 19:21

It is worth reading in more depth how DFFTs work but you should always have the following formulae in mind. For a time series with n points and maximum time Tmax, the time resolution is given by dt = Tmax / n

A DFFT will produce n points with

Fmax = 1 / dt

dF = 1 / Tmax

You seem to suggest the maximum frequency is sufficient (so the time resolution is okay) but the frequency resolution isn't good enough: you need to collect more data, at the same time resolution.

score 0 · Answer 2 · edited May 23 '17 at 12:19

If (1) the sampling time is too short, (2) you require higher estimation frequency accuracy, and, (3) you know that your signal is a sine wave, then you can fit the signal to a sine wave. Like in How do I fit a sine curve to my data with pylab and numpy?, with the exception that the frequency needs to be added.

Here is an example figure with a frequency of around 8 MHz:

Figure with fitted sine wave

Below is the example code:

""" Modified from https://stackoverflow.com/a/16716964/6036470 """
from numpy import sin, linspace, pi,average;
from pylab import plot, show, title, xlabel, ylabel, subplot, scatter
from scipy import fft, arange, ifft
import scipy
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import leastsq

ff = 8e6;   # frequency of the signal
Fs = ff*128;  # sampling rate
Ts = 1.0/Fs; # sampling interval

t = arange(0,((1/ff)/128)*(128)*5,Ts) # time vector
A = 2.5;

ff_0 = 8.1456e6
y = A*np.sin(2*np.pi*ff_0*t+15.38654*pi/180) + np.random.randn(len(t))/5

guess_b = 0
guess_a = y.std()*2**0.5;
guess_c = 10*pi/180
guess_d = ff*0.98*2*pi

fig = plt.figure(facecolor="white")
plt.plot(t,y,'.', label='Signal Fred. %0.4f Hz'%(ff_0/1e6))
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.grid(alpha=0.5);

optimize_func = lambda x: (x[0]*np.sin(x[2]*t+x[1]) - y);
est_a,  est_c, est_d = leastsq(optimize_func, [guess_a, guess_c, guess_d])[0]
data_fit = est_a*np.sin(est_d*t+est_c) ;
plt.plot(t,data_fit,label='Fitted Est. Freq. %0.4f Hz'%(est_d/(2*pi)/1e6))
plt.legend()
plt.tight_layout();
plt.show();

fig.save("sinfit.png")

Frequency resolution issue using FFT in numpy

2 Answers2