Why are FFT results different from theory and how to eliminate the difference?

Question

I am working with time series data, and applying FFT to calculate the power spectrum. The data is something like this:

$y = \sum_{k=1}^{5} \cos(10 k x)$

All the peaks in the FFT output should have identical magnitudes, similar to the pattern shown in the graph below:

However, the FFT output I am getting is different, as shown in the graph below:

The magnitudes of the peaks in my data slightly differ. Increasing the sampling rate or using a Hamming window reduces but doesn't fully eliminate these differences. How can I completely remove these differences?

N = 1000
t = np.linspace(0,2*np.pi,N)
y = np.zeros(N)
for freq in [10,20,30,40,50]:
    y = y + np.cos(freq * t)
h = scipy.signal.windows.hamming(N)
yfft = np.abs(np.fft.fft(y*h) / N)**2
xfft = np.hstack([np.arange(0,N//2), np.arange(-(N//2),0)])
d = np.array(sorted(zip(xfft,yfft))).reshape(-1,2)
xfft = d[:,0]
yfft = d[:,1]
plt.plot(xfft,yfft,color="red")
plt.show()
for freq in [10,20,30,40,50]:
    print(f"freq={freq} magnitude={yfft[xfft==freq]}")

Answer 1

The ideal Fourier transform here would be the sum of 10 Dirac delta functions that has spikes with infinite peaks. It’s not surprising you didn’t get infinite peaks in practice (because that’s impossible), and the width (and so height) of your peaks depends on the specific implementation details and precision of your FFT. If you integrate your FFT’s individual peaks, you’ll probably get totals which are close to each other.

If you really want results that match theory closely, use finite functions. I.e. the (unnormalized) Fourier transform of a limited height pulse is $\sin(x)/x$, so you can use something like $\sum_{k=1}^5 (\sin(x)/x)\cos(10kx)$ which should give spikes with a mostly consistent width/height, but you’ll still get some small differences. You can bound these differences using various theorems related to sampling error in signal analysis.

Answer 2

In this particular ideal case, you need to make sure that your cosines are "commensurate" with your domain. I will explain what this means.

By using np.linspace with endpoint=True (which is its default value) you are sampling the boundary value twice. This is because the FFT models functions as being periodic. For example if your domain has length $T$, then the FFT acts as if $y(t)=y(t+T)$ for all $t$. Your domain should be $[0,T)$, i.e. excluding one endpoint. If you include both endpoints, like in $[0,T]$, you will sample the point $t=0$ twice. This is because $y(0)=y(T)$, again because of periodicity. If you do this, the function is not 'commensurate' with your domain, which just means that it does not fit nicely. You can solve this by having endpoint=False

Here is some code that draw the spectrum of a single sine wave. First I plot it using endpoint=True and then using endpoint=False, the latter of which gives the correct result.

import numpy as np
import matplotlib.pyplot as plt

N = 8
t = np.linspace(0, 2*np.pi, endpoint=False, num=N)
y = np.sin(t)

fig, axes = plt.subplots(ncols=2, figsize=(8,3), dpi=200)

axes[0].plot(t, y, 'o-')
axes[0].set_xlabel('$t$')
axes[0].set_ylabel('$y$')

yfft = np.fft.fft(y)
omega = np.fft.fftfreq(N, d=(t[1]-t[0]))

yfft = np.fft.fftshift(yfft)
omega = np.fft.fftshift(omega)

axes[1].plot(omega, np.abs(yfft)**2, 'o-')
axes[1].set_xlabel('$\omega$')
axes[1].set_ylabel(r'$\hat y$')

plt.tight_layout()
plt.show()

Here the function is incommensurate, with endpoint=True

Here the function is commensurate, with endpoint=False

Notice how nicely the second FFT behaves. Indeed, all the entries that you want to be zero are actually zero.

To drive my point home even more, I tile my function $y$ 2 times, as if it were periodic. Here you see why you get the slightly off Fourier transform.

One more important question: how do you solve this in the real world, where you have no control over your domain? Well, here's where things get ugly and you have to consider your options. You can ensure you sample enough periods, have enough sample points and you could also reduce edge-effects by using a window. You are using all three of these options and it shows. If it was a sample of real data, it would be really good. Even though your functions are not commensurate.