0
$\begingroup$

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.

This is the time series of all the data.


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

elevation=np.loadtxt('/content/S1S2S3.txt')

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[1]-time[0]
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. enter image description here

enter image description here

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Welcome to SE.SP! Please do not add "answers" that do not answer your question. Please edit your question to add new information. $\endgroup$
    – Peter K.
    Apr 13, 2022 at 16:08

1 Answer 1

1
$\begingroup$

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.

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

Your Answer

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

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