4

I am trying to determine the most dominant frequency of a signal. However, when artificially creating a 50 Hz signal and applying sufficient zeropadding to enhance fft resolution, I get a top frequency of 49,997 Hz. For my application this is a significant difference. Did I do something wrong here?

import numpy as np
import matplotlib.pyplot as plt

fs = 2**12

x = np.linspace(0,1,fs+1)

signal = np.sin(50*2*np.pi*x)
spect = abs(np.fft.fft(np.append(signal,np.zeros(999*fs))))

plt.figure('Four Coef')
plt.plot(spect)
plt.axis([49995,49999,2048.01,2048.05])

plt.show()

Note that coefficient 49997 corresponds to a frequency of 49,997 Hz due to the zero-padding.

Edits: The array represents exactly 1 seconds of 50 Hz signal. The last 999 seconds are zeros in order to increase the fft "resolution" to 1 mHz. I have only 1 second of signal available, from which i need the top frequency, accurate up to the mHz

Changing the the sample rate fs = 2**8 gives a maximum of 49.999 so i suppose the way of sampling is critical here...

Tom Baksteen
  • 41
  • 6

2 Answers2

2

You are not taking the FFT of 1000 s of a 50 Hz wave: the array you pass to np.fft.fft is 1 second of signal followed by 999 seconds of silence zeros). So your clipped signal FFTs to a funky, multi-peaked thing.

When I do the following with a continuous signal, I see the peak at index 50000 as expected:

import numpy as np
import matplotlib.pyplot as plt

fs = 2**12

x = np.linspace(0,1000,fs*1000)

signal = np.sin(50*2*np.pi*x)
spect = abs(np.fft.fft(signal))

plt.figure('Four Coef')
plt.plot(spect)
print np.argmax(spect), np.max(spect)

plt.show()

Output:

50000 2047497.79244

NB1/ just repeating your array won't work properly either, because the ends won't "match up" (the signal will jump from the end of one 1 s array to the beginning of th next).

NB2/ You might consider using rfft and rfftfreq to get the frequencies here.

xnx
  • 24,509
  • 11
  • 70
  • 109
  • The 999 seconds of silence are essential to investigate at mHz resolution. I have to deal with signals of only 1 second, not 1000 seconds unfortunately. The funky peaks are called spectral leakage. – Tom Baksteen Jul 16 '15 at 13:11
  • Ah - I misunderstood what you were doing, then. You might look to the other answer instead. By the way, `linspace(0,1,fs+1)` will return an array of `fs+1` points, which may (or may not) be what you intended here. – xnx Jul 16 '15 at 13:22
  • 1
    @xnx For a sample rate of `fs`, `linspace(0, 1, fs+1)` is correct. Alternatively, `linspace(0, 1, fs, endpoint=False)` could be used. – Warren Weckesser Jul 16 '15 at 13:31
  • I notice I should have explained better, thanks! Indeed I want `fs+1` since the endpoint is also included. – Tom Baksteen Jul 16 '15 at 13:33
0

The frequency of the peak magnitude FFT result bin will indicate the frequency of a signal only if the length of the FFT is an exact integer multiple of the period of the frequency of that signal. For any other frequency, try using a frequency estimation algorithms (such as parabolic or Sinc interpolation of the FFT result, but there are many other estimation methods) instead of just the raw single peak magnitude bin frequency.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153
  • Using the quadratic interpolation technique I got the same results. I intentionally left it out for simplicity. – Tom Baksteen Jul 19 '15 at 07:01
  • Print out the length of your fft as a number. – hotpaw2 Jul 19 '15 at 07:28
  • The length of spect is 2048001. Index 0 represents 0 Hz and index 2048000 is 2048 Hz which is exactly fs/2 The simpler DFT implementation as a matrix vector mulitplication yields the same result, both with and without quadratic interpolation. – Tom Baksteen Jul 20 '15 at 08:19
  • Try an FFT length of N=2048000, where index N/2 represents Fs/2 and index 0 is Fs (and DC). That length will be an integer multiple of the period of 50 Hz. – hotpaw2 Jul 20 '15 at 12:31
  • No difference. Also, when using a 1 Hz signal in stead of 50, the estimation error is about 50 times larger. Maybe that info helps? – Tom Baksteen Jul 20 '15 at 13:15