# How to draw the PSD from a time series

I try to draw the spectral density of a time series in order to compare it with the theoretical one. Please can any one help me to do this.

Thanks Dan Boschen for your response.

About "scipy.signal.fftpack as fft" is not worked.

import scipy.signal.fftpack as fft


ModuleNotFoundError: No module named 'scipy.signal.fftpack'

So I use Scipy.fft instead of scipy.signal.fftpack, it worked but the result was different from the expected.

Please can you check the code.

enter preformatted text here

import scipy
import numpy as np
import matplotlib.pyplot as plt

time=elevation[200:11070,0] # to extract the first column's value between the ligne 200 to 11070
H=elevation[200:11070,1] # to extract the second column's value between the ligne 200 to 11070

T = time-time
fs = 1/T
N = time.size
M=10870
freq=scipy.fft.fftshift(np.fft.fftfreq(M,T))
x=H
X = 1/fs * scipy.fft.fftshift(np.fft.fft(x,M))
P = np.abs(X)**2 * fs /N

# to plot
plt.plot(freq,P)
#plt.plot(P)
plt.xlim(0,0.4)
plt.title('spectral density')
plt.xlabel('frequency,f,[Hz]')
plt.ylabel('Spectral density [m**2 s]')
plt.show()


the link for the data file is here.

Best regard.

There is my result with blue curve, and it will be the superposition of the two figures with the red curve.   – Peter K.
Apr 13, 2022 at 16:08

The power spectral density is given as the Fourier Transform of the autocorrelation function. This can be determined using an FFT as the complex conjugate product of the FFT result. Below is an example in Python using scipy.signal.fftpack as fft:

X = 1/fs * fft.fftshift(fft.fft(x))
P = np.abs(X)**2 * fs / N


Where $$fs$$ is the sampling rate, and $$N$$ is the total number of samples. Zero padding $$x$$ will serve to interpolate more samples in the spectrum and can be done using fft.fft(x, M) with $$M>>N$$.

Another approach is to use Welch's method which results in less noise in the power estimate at the expense of frequency resolution. See scipy.signal.welch and this link for further details.

• Thanks so much Dan Boschen. Apr 12, 2022 at 11:58
• @DanBoschen wouldn't the PSD be P = np.abs(X)**2 / (fs*N) ?
– Jdip
Sep 2, 2022 at 12:56