6
$\begingroup$

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:

expected FFT output

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

actual FFT output

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]}")
$\endgroup$
1
  • 1
    $\begingroup$ Can you share your code? The number of samples can affect what the FFT looks like. $\endgroup$
    – JimmyK4542
    Feb 20 at 3:54

2 Answers 2

12
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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 enter image description here

Here the function is commensurate, with endpoint=False enter image description here

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.

enter image description here

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .