Scipy FFT and Numpy FFT disagree on pulse train spectrum?

Question

I am doing an FFT on a series of pulses. The series is one pulse of amplitude 1 every 7 days over a total of 367 days.

When I run the following code:

import numpy as np
import pandas as pd
from scipy.fft import fft, fftfreq, fftshift, ifft
from scipy.signal import blackman
from matplotlib import pyplot as plt
import random

## Signal 
num_samples = 367
# time in days
t = np.arange(int(num_samples))
# Amplitude and position of pulse. Amplitude here is 0 or 1 but can generate random values
# Position here is every 7th day
signal = [random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#np.sin(2*np.pi*5*t/N)#[random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#

# FFT and IFFT using Numpy

sr = 367
X = np.fft.fft(signal)
n = np.arange(num_samples)
T = num_samples/sr
freq = n/T 

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Numpy')
plt.stem(freq, np.abs(X), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')

plt.subplot(122)
plt.title('IFFT using Numpy')
plt.plot(t, np.fft.ifft(X), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# FFT and IFFT using Scipy

sp = fft(signal)
freq = fftfreq(t.shape[-1])

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Scipy')
plt.stem(freq, np.abs(sp), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |sp(freq)|')

plt.subplot(122)
plt.title('IFFT using Scipy')
plt.plot(t, ifft(sp), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

I get the following:

Clearly there are shifting and scaling issues but more importantly I expected the fft of a pulse train to be a series of uniform peaks in the frequency spectrum. I don't understand the peaks that result which means I'm probably misunderstanding how the functions are interpreting the signal. Any guidance would be appreciated.

user_14062403 · Answer

You should have used mp.fft.fftshift to plot the numpy fft. Also, there is an easier way to make the frequ axis the numpy case.

Your code with modifications indicated by "<<<" and *>>>":

import numpy as np
import pandas as pd
from scipy.fft import fft, fftfreq, fftshift, ifft
from scipy.signal import blackman
from matplotlib import pyplot as plt
import random

## Signal 
num_samples = 367
# time in days
t = np.arange(int(num_samples))
# Amplitude and position of pulse. Amplitude here is 0 or 1 but can generate random values
# Position here is every 7th day
signal = [random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#np.sin(2*np.pi*5*t/N)#[random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#

# FFT and IFFT using Numpy

sr = 367

X = np.fft.fft(signal)  
freq = np.fft.fftfreq(len(t), d=1) # <<< Let numpy build the frequency axis for you >>>

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Numpy')

# <<< Shift the zero-frequency component to the center of the spectrum on frequ and X >>>
# <<< Note: use_line_collection=True in plt.stem() removes a warning (not important) >>>
plt.stem(np.fft.fftshift(freq), np.fft.fftshift(np.abs(X)), 'b', markerfmt=" ", basefmt="-b", use_line_collection=True)
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')

plt.subplot(122)
plt.title('IFFT using Numpy')
plt.plot(t, np.fft.ifft((X)), 'r') 
plt.plot(t, np.fft.ifft((X)), 'r') 

plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# FFT and IFFT using Scipy

sp = fft(signal)  

freq = fftfreq(len(t))
print(freq.shape)

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Scipy')
plt.stem(freq, np.abs(sp), 'b', markerfmt=" ", basefmt="-b",use_line_collection=True)
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |sp(freq)|')

plt.subplot(122)
plt.title('IFFT using Scipy')
plt.plot(t, ifft(sp), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

OUTPUT:

Scipy FFT and Numpy FFT disagree on pulse train spectrum?

1 Answers