How to properly scale frequency axis in Fast Fourier Transform?

Question

I am trying some sample code taking the FFT of a simple sinusoidal function. Below is the code

import numpy as np
from matplotlib import pyplot as plt

N = 1024
limit = 10
x = np.linspace(-limit, limit, N)
dx = x[1] - x[0]
y = np.sin(2 * np.pi * 5 * x) + np.sin(2 * np.pi * x)
Y = np.abs(np.fft.fft(y) ** 2)
z = fft.fftshift(np.fft.fftfreq(N, dx))
plt.plot(z[int(N/2):], Y[int(N/2):])
plt.show()

From the function that is given, , it is clear there should be two spikes at frequencies 1 and 5. However, when I run this code, I get the following plot.

                                         

Clearly the spikes are not where they should be. Additionally, I have noticed that the frequency scaling is sensitive to the number of points N as well as the interval limits that I make limit. As an example, setting N = 2048 gives the following plot.

                                         

As you can see, the locations of the spikes have changed. Now keeping N = 1024 and setting limit = 100 also changes the result.

                                         

How can I make it so the frequency axis is properly scaled at all times?

Answer 1

fftfreq returns the frequency range in the following order: the positive frequencies from lowest to highest, then the negative frequencies in reverse order of absolute value. (You usually only want to plot one half, as you do in your code.) Note how the function actually needs to know very little about the data: just the number of samples and their spacing in the time domain.

fft performs the actual (Fast) Fourier transformation. It makes the same assumption about the input sampling, that it's equidistant, and outputs the Fourier components in the same order as fftfreq. It doesn't care about the actual frequency values: the sampling interval is not passed in as a parameter.

It does however accept complex numbers as input. In practice, this is rare. Input is usually samples of real numbers, as in the above example. In that case, the Fourier transform has a special property: it's symmetric in the frequency domain, i.e. has the same value for f and −f. For that reason, it often doesn't make sense to plot both halves of the spectrum, as they contain the same information.

There is one frequency that stands out: f = 0. It's a measure of the average value of the signal, its offset from zero. In the spectrum returned by fft and the frequency range from fftfreq, it's at the very first array index. If plotting both halves, it may make sense to shift the frequency spectrum around, so that the negative half is to the left of the zero-component, the positive half to its right, meaning all values are in ascending order and ready to be plotted.

fftshift does exactly that. However, if you plot only half of the spectrum anyway, then you may as well not bother doing this at all. Though if you do, you must shift both arrays: the frequencies and the Fourier components. In your code, you only shifted the frequencies. That's how the peaks ended up on the wrong side of the spectrum: You plotted the Fourier components referring to the positive half of the frequencies with respect to the negative half, so the peaks on the right are actually meant to be close to zero, not at the far end.

You don't really need to rely on any of those functions operating on the frequencies. It is straightforward to generate their range based on the documentation of fftfreq alone:

from numpy.fft import fft
from numpy import arange, linspace, sin, pi as π
from matplotlib import pyplot

def FFT(t, y):
    n = len(t)
    Δ = (max(t) - min(t)) / (n-1)
    k = int(n/2)
    f = arange(k) / (n*Δ)
    Y = abs(fft(y))[:k]
    return (f, Y)

t = linspace(-10, +10, num=1024)
y = sin(2*π * 5*t) + sin(2*π * t)
(f, Y) = FFT(t, y)
pyplot.plot(f, Y)
pyplot.show()

Note that NumPy also offers dedicated functions, rfft and rfftfreq, for the common use case of real-valued data.

Answer 2

import numpy as np
from matplotlib import pyplot as plt

N = 1024
limit = 10

x = np.linspace(-limit, limit, N)
dx = x[1] - x[0]

y = np.sin(2 * np.pi * 5 * x) + np.sin(2 * np.pi * x)
yhat = np.fft.fft(y)

Y = np.abs(yhat)
freq = np.linspace(0.0, 1.0/(2*dx), N//2)

plt.plot(freq, Y[0:N//2]*(2/N))
plt.show()