1
$\begingroup$

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. The code below is what I run:

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()

This results in the following plots:

enter image description here

enter image description here

So I'm confused as to what is happening:

  • I did not expect any peaks but rather a 'comb' b/c it's fft on pulse train
  • if these are correct peaks then I would expect the closest to 0 frequency (7-day period) to taller than the ones to its right
  • it looks like there is a scaling/shifting issue

Any guidance would be appreciated.

$\endgroup$
2
  • 2
    $\begingroup$ You should read the docs on the modules, it's all there. Just some hints: Hz and s as units makes no sense in the context of days. Use the rfft modules, as your input is real. Your x axis is incorrect for numpy output, but correct for scipy (aside from Hz which should read days). $\endgroup$
    – Max
    Commented Dec 16, 2021 at 8:52
  • 1
    $\begingroup$ That's the advice I always give :-D But I did read the docs and scipy's code (yay github and OS) and can figure out the scaling, etc. But the overall response was still a mystery. I did this experiment because I've been seeing some 'blind' use of the libraries in-house and began questioning results folks were getting ... so I tried blind use of the libraries myself. As it always has been GIGO. $\endgroup$
    – MikeB2019x
    Commented Dec 16, 2021 at 13:52

1 Answer 1

2
$\begingroup$

If you change your line:

plt.stem(freq, np.abs(X), 'b', markerfmt=" ", basefmt="-b")

to

plt.stem(freq, fftshift(np.abs(X)), 'b', markerfmt=" ", basefmt="-b")

the graphs will be the same for both.

The reason for the "extra" spikes is that the length of the signal (367) is not divisible by 7. If that is increased to 371, then I get the following graphs is is probably closer to what you're expecting.

First picture, redone

Second picture, redone


Code update to see interpolation effect

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 = 371
# time in days
t = np.arange(int(num_samples))
t3 = np.arange(int(num_samples)*3)
# 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, num_samples*3)
n = np.arange(len(X))
T = num_samples/sr
freq = n/T 

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Numpy')
plt.stem(freq, fftshift(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)[0:num_samples], 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# FFT and IFFT using Scipy

sp = fft(signal, num_samples*3)
freq = fftfreq(t.shape[-1]*3)

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)[0:num_samples], 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

Interpolation 1

Interpolation 2

$\endgroup$
3
  • $\begingroup$ So padding is required to avoid windowing effects. I guess a challenge is if the pulse train is irregular ie. the time between pulses varies occasionally. tnx! $\endgroup$
    – MikeB2019x
    Commented Dec 16, 2021 at 15:59
  • $\begingroup$ @MikeB2019x It's a little more than windowing effects. It's more like "interpolation" or "zero padding" effects, but yes, if the pulse train is irregular, then things will definitely be different. See added example for 3 x interpolation. $\endgroup$
    – Peter K.
    Commented Dec 16, 2021 at 16:35
  • 1
    $\begingroup$ I have asked a follow-up question regarding varying periodicity here: dsp.stackexchange.com/questions/80632/… $\endgroup$
    – MikeB2019x
    Commented Dec 16, 2021 at 17:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.