0
$\begingroup$

I am trying to generate a time series from a defined PSD function, however i tried to do this in python , with the following steps:

  1. Define the Power Spectral Density
  2. Define the time parameters
  3. Define the frequency range
  4. Descritize the PSD
  5. Single-Sided and Normalized DFT
  6. Generate random phase
  7. Calculate Discrete Fourier Transform (DFT)
  8. Perform the inverse Fourier transform (IFFT).

However , at the end when i try the make the comparison between the defined the PSD and the one generated from the time serie , i do get a significant difference , as showed in the following plot :

enter image description here

And this is the code that i am using in python:

import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import interp1d import control from scipy.signal import welch import scipy.signal import random import cmath # ------------------------------------------------------------------------------------------ # Define the Power Spectral Density Curve # ------------------------------------------------------------------------------------------ # Define frequency range # Define frequency range freq = np.linspace(1, 60, 1000) # Initialize PSD array psd = np.zeros_like(freq) # Linearly increasing PSD from 0 to 1 Hz psd[:100] = np.linspace(2, 2.4, 100) # Constant PSD from 1 to 2 Hz psd[100:200] = 2.4 # Linearly decreasing PSD from 2 to 10 Hz psd[200:] = np.linspace(2.4, 0.001, 800) # Plot PSD plt.figure(figsize=(10, 6)) plt.plot(freq, psd, color='blue') plt.title('Power Spectral Density') plt.xlabel('Frequency (Hz)') plt.ylabel('PSD (dB/Hz)') plt.grid(True) plt.show() # ------------------------------------------------------------------------------------------ # Define the time parameters # ------------------------------------------------------------------------------------------ t_fin = 4 Fr = 1 / t_fin # ------------------------------------------------------------------------------------------ # Define the frequency range # ------------------------------------------------------------------------------------------ F0 = 4 F_max = 50 Fs = 2 * freq[-1] # Generate time vector t = np.arange(0, t_fin+1/Fs, 1/Fs) t = t[1:-1] # ------------------------------------------------------------------------------------------ # Descritize the PSD # ------------------------------------------------------------------------------------------ # Sampling PSD f = np.arange(0, Fs/2+ Fr, Fr) PSD = interp1d(freq, psd, fill_value="extrapolate")(f) PSD[f < F0] = -np.inf PSD[f > F_max] = -np.inf # Plot PSD curve plt.figure(figsize=(10, 6)) plt.plot(freq, psd, color='blue', label='Original PSD') plt.scatter(f, PSD, color='red', label='Sampled PSD') plt.title('Power Spectral Density (PSD) of Simplified Seismic Data') plt.xlabel('Frequency (Hz)') plt.ylabel('Power/Frequency') plt.grid(True) plt.xlim([0, 60]) plt.legend() plt.tight_layout() plt.show() # ------------------------------------------------------------------------------------------ # Single-Sided and Normalized DFT # ------------------------------------------------------------------------------------------ # Initialize plampl array plampl = np.zeros_like(PSD) PSD= control.db2mag(PSD) # Calculate plampl plampl[0] = np.sqrt(Fr * PSD[0]) plampl[1:] = np.sqrt(2 * Fr * PSD[1:]) # ------------------------------------------------------------------------------------------ # Generate random phase # ------------------------------------------------------------------------------------------ plphase = np.random.randn(1, len(plampl)) * (2 * np.pi) # From Single-Sided Ampl/Phase to Double-Sided Real/Imag plampl[0:-1] = plampl[0:-1] / 2 plreal = plampl * np.cos(plphase) plimag = plampl * np.sin(plphase) # Construct p2 array p2 = np.concatenate((plreal + 1j * plimag, np.conj(plreal[0:-1] + 1j * plimag[0:-1])[::-1])) # ------------------------------------------------------------------------------------------ # Calculate Discrete Fourier Transform (DFT) # ------------------------------------------------------------------------------------------ DFT = len(p2) * p2 # ------------------------------------------------------------------------------------------ # Perform the inverse Fourier transform (IFFT) # ------------------------------------------------------------------------------------------ s = np.fft.ifft(DFT) s_real = s.real.squeeze() # Remove any singleton dimensions s_imag = s.imag.squeeze() sqrt_sum = np.sqrt(s_real**2 + s_imag**2) #time_vector = np.arange(0, len(s) / Fr, len(s_real)) #print(time_vector) # Extract the real part of s # Plotting plt.plot(np.arange(len(s_real)), s_real) plt.xlabel('Length') plt.ylabel('Real part of s') plt.title('Real part of s vs Length') plt.grid(True) plt.show() # Calculate the Power Spectral Density (PSD) using Welch's method f, Pxx_den = welch(s.real, fs=Fs, nperseg=len(s_real)) Pxx_den = np.squeeze(Pxx_den) # Flatten Pxx_den if necessary plt.semilogy(f, Pxx_den) plt.plot(freq, psd, color='red', label='Original PSD') plt.xlabel('frequency [Hz]') plt.legend(['PSD from Time Series', 'Original PSD']) plt.ylabel('PSD [V**2/Hz]') plt.xlim(0, 100) # Set the frequency range from 0 to 100 Hz plt.yscale('log') # Set the x-axis to logarithmic scale plt.show()

Any Help will be much appreciated .

Improve this question
$\endgroup$
6
  • $\begingroup$ Hi and welcome to our stack exchange. Please fix your code so that it actually runs. You are missing the import statements and the names DFT and welch are undefined, some lines need a comment character and I don't see a random phase anywhere in there. $\endgroup$ Commented May 13, 2024 at 12:46
  • 1
    $\begingroup$ Does this answer your question? Generate a time series from power spectral density $\endgroup$ Commented May 13, 2024 at 13:28
  • $\begingroup$ @Hilmar Thank you for you reply , i have added the missing commands shoud be working now . $\endgroup$ Commented May 13, 2024 at 14:04
  • $\begingroup$ @Marcus Müller , i think it is the similar procedure but unfortuntly it is not helpfull. i am getting at the end a different PSDs. $\endgroup$ Commented May 13, 2024 at 14:06
  • $\begingroup$ yes, because you're doing nonsensical things, like having infinities in your PSD. That's why I referred you to a more sensical example :) $\endgroup$ Commented May 13, 2024 at 14:07

1 Answer 1

Reset to default
3
$\begingroup$

Sorry, this is a rather tortured implementation of something that's pretty straight forward. It's hard to point out what exactly is wrong in your code since most of it seems nonsensical to me and there aren't enough comments to assess why you choose this design.

The process is simple enough.

  1. Sample the PSD or magnitude spectrum on a FFT frequency grid
  2. If it's an actual PSD, take the square root
  3. Add a random phase
  4. Make sure the spectrum is conjugate symmetric
  5. Take the inverse FFT

Below is on example for noise that's pink above 100Hz, flat below 100Hz and sampled at 48 kHz.

enter image description here

And here is the code

import numpy as np import matplotlib.pyplot as plt from scipy import fft from scipy import signal # %% make soom pink noise, flat below 100 Hz nx = 16384 # FFT length fc = 100 # cut off in Hz fs = 48000 # sample rate fpink = np.ones(nx // 2 + 1, dtype=np.complex128) # init # find the index of of the cutoff frequency i100 = round(fc / (fs / nx)) # make the magnitide specturm proprotional to 1/x above the cutoff fpink[i100:] *= i100 / np.arange(i100, nx / 2 + 1) # add random phase except at DC and Nyquist fpink[1:-1] *= np.exp(1j * np.random.rand(nx // 2 - 1) * 2 * np.pi) # inverse FFT using a conjugate symmetric spectrum xpink = fft.irfft(fpink) # %% plot the stuff plt.clf() fr_axis = np.arange(1, nx / 2 + 1) * fs / nx plt.semilogx(fr_axis, 20 * np.log10(fpink[1:]), linewidth=3, label='Target') # do pwelch analysis frame_size = 2048 fr, fw = signal.welch(xpink, fs=fs, nperseg=frame_size) fw *= fs * nx / 2 plt.semilogx(fr, 10 * np.log10(fw), label='Welch estimate') plt.grid(True) plt.xlabel('Frequency in Hz') plt.ylabel('Level in dB') plt.legend()
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.