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.

Answer 1

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: