Generate a Time Series from Power Spectral Density Python

Question

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 :

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 .

Answer 1

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.

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