0
$\begingroup$

I am trying to find out the dominating frequency of a signal with a frequency of 50 Hz (sampled at 200 Hz - every 5 milliseconds).

The python code I am using to do this is the following (based on this):

import numpy
from numpy import sin
from math import pi

t = numpy.linspace(0, 1, 201)  # 200 Hz sampling rate
y = sin(2*pi*t*50)

fourier = numpy.fft.fft(y)
frequencies = numpy.fft.fftfreq(len(t), 0.005)  # where 0.005 is the inter-sample time difference
positive_frequencies = frequencies[numpy.where(frequencies > 0)]
magnitudes = abs(fourier[numpy.where(frequencies > 0)])  # magnitude spectrum

peak_frequency = numpy.argmax(magnitudes)

The peak frequency I get is 49 Hz, not 50 Hz (what I think I should obtain).

Is this the correct way to calculate the peak frequency? If yes, then can you please explain why this happens? Where is my mistake?

EDIT Code modified according to SergV's answer:

import numpy
from numpy import sin
from math import pi

Fs=200.
F=50.
t = [i*1./Fs for i in range(200)]
y = sin(2*pi*numpy.array(t)*F)

fourier = numpy.fft.fft(y)
frequencies = numpy.fft.fftfreq(len(t), 0.005)  # where 0.005 is the inter-sample time difference
positive_frequencies = frequencies[numpy.where(frequencies > 0)]  
magnitudes = abs(fourier[numpy.where(frequencies > 0)])  # magnitude spectrum

peak_frequency = numpy.argmax(magnitudes)

The calculated peak frequency is still 49 Hz...

EDIT 2: It seems that the problem was that the positive frequencies were not including 0. After changing the code to

positive_frequencies = frequencies[numpy.where(frequencies >= 0)]  
magnitudes = abs(fourier[numpy.where(frequencies >= 0)])

the calculated dominating frequency was 50 Hz (as expected).

This is the magnitude spectrum I get: http://imgur.com/JxhnNc5 It shows that the dominating frequency is 50 Hz. Interestingly, even when, because of the bug, the calculated peak frequency was 49 Hz, the magnitude spectrum was showing the correct dominating frequency (50 Hz)

$\endgroup$
2
$\begingroup$

The size of $y$ is 201. So frequency bin of fft is not 1 Hz. Change code for calculating $t$:

Fs=200.
F=50.
t = [i*1./Fs for i in range(200)]

Now you have only 1 peak on 50 Hz.

P.S.

  • Plot magnitude spectrum
  • Try to modify F (49, 49.5, 49.67). Enjoy!

Welcome to DSP and new questions!

$\endgroup$
  • $\begingroup$ Thank you for your answer. I modified the code according to your suggestion (see the edit), but the calculated peak frequency is still 49 Hz. $\endgroup$ – al1na Apr 22 '15 at 12:00
  • $\begingroup$ @al1na. I have not numpy now on computer. Possible suggestions - you have to convert from index ( peak_frequency) to frequency, so may be condition for positive_frequencies must be (frequencies >= 0), but not (frequencies > 0). And PLOT magnitude spectrum! $\endgroup$ – SergV Apr 22 '15 at 12:22
  • $\begingroup$ Seems that changing the condition for the positive frequencies to include 0 fixed the problem. I also added to the question the magnitude spectrum. Thank you! $\endgroup$ – al1na Apr 22 '15 at 13:10
  • $\begingroup$ @al1na, If you want to get new knowledge of FFT, plot magnitude spectrum for F=49.67. But it will be another question. $\endgroup$ – SergV Apr 22 '15 at 14:00
2
$\begingroup$

To get the frequency of an FFT result bin, you need to multiply the bin number by the sample rate divided by the length of the FFT.

Note that a sinusoid of a frequency between FFT result bins will get its energy distributed among several FFT result bins, not just the one with the peak magnitude or the one closest. So some interpolation may be required.

$\endgroup$
  • $\begingroup$ Regarding the last point in the answer: Some information about how to interpolate can be found here and here, for example. Furthermore, there's another method to increase the accuracy of the frequency estimation, that, however, is based on certain assumptions about the signal. $\endgroup$ – applesoup Aug 25 '15 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.