**Define the Power Spectral Density Curve** freq = np.linspace(1, 60, 1000) psd = np.zeros_like(freq) psd[:100] = np.linspace(2, 2.4, 100) psd[100:200] = 2.4 psd[200:] = np.linspace(2.4, 0.001, 800) 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 **Perform the inverse Fourier transform (IFFT)** s = np.fft.ifft(DFT)   s_real = s.real.squeeze()   s_imag = s.imag.squeeze()   sqrt_sum = np.sqrt(s_real**2 + s_imag**2) 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)   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)   plt.yscale('log')   plt.show() 

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

I am trying to generate a time serie 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)

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

Define the Power Spectral Density Curve

freq = np.linspace(1, 60, 1000)

psd = np.zeros_like(freq)

psd[:100] = np.linspace(2, 2.4, 100)

psd[100:200] = 2.4

psd[200:] = np.linspace(2.4, 0.001, 800)

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]

t = np.arange(0, t_fin+1/Fs, 1/Fs)

t = t[1:-1]

Descritize the 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

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

plampl = np.zeros_like(PSD)

PSD= control.db2mag(PSD)

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)

plampl[0:-1] = plampl[0:-1] / 2

plreal = plampl * np.cos(plphase)

plimag = plampl * np.sin(plphase)

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

s_imag = s.imag.squeeze()

sqrt_sum = np.sqrt(s_real2 + s_imag2)

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)

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)

plt.yscale('log')

plt.show()

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

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

**Define the Power Spectral Density Curve** freq = np.linspace(1, 60, 1000) psd = np.zeros_like(freq) psd[:100] = np.linspace(2, 2.4, 100) psd[100:200] = 2.4 psd[200:] = np.linspace(2.4, 0.001, 800) 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 **Perform the inverse Fourier transform (IFFT)** s = np.fft.ifft(DFT) s_real = s.real.squeeze() s_imag = s.imag.squeeze() sqrt_sum = np.sqrt(s_real**2 + s_imag**2) 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) 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) plt.yscale('log') plt.show()