How to generate time-series from a given one-sided PSD?

Question

Generating time series from a given PSD has been discussed in this and some other forums already (which I've referred some of them below) but almost all of them are mostly descriptive and none of them gives the exact details or code that implements and verifies it. As always the devil is in the details. So I went through the discussions in those pages and tried to follow them, but at the end I didn't get what I wanted.

Two methods of calculating time series from PSD is discussed in the following page, which I tried to implement both of them in my MATLAB code below but none of them seems to be working properly. I think the conversion from PSD to amplitude response is not precise in one of the methods and scaling can be the issue in the other one.

This page mentions a more correct formula for converting PSD to amplitude response.

Also, for verifying the methods, I extract PSD out of the generated time series to make sure it matches the original PSD. For that, I'm following the instructions given in the following page from Mathworks, also I am using periodogram() built-in function from Matlab and I compare them (but I get different outcomes).

I see that in the methods I have used, the estimated PSD does not match the original PSD. So I am not sure where I am exactly wrong in my code.

The Matlab code below needs a .mat file to read the one-sided PSD as input to get the 'freq' and 'psd_measured' vectors. I have provided those two vector values in the following link. Simply copy and paste them at he beginning of the code to replace load('psd_freq.mat') line so that you have those values ready (but generally you should be able to use any PSD as the input)

psd_freq.mat on pastebin.

clc
close all
clear all

load('psd_freq.mat')

asd = psd_measured;

% Since the data is given in v/sqrt(Hz) (i.e. it's in amplitude), we need to square it to get the psd. 
psd = asd .^ 2; 

noise_freq  = 100*1e9;

i_noise_req = find(freq <=  noise_freq, 1, 'last');
i_noise_req = i_noise_req(1);

fsym         = 53.125e9; % baud rate
oversampling = 16;       % number of samples per ui
block_size   = 1024;

si    = 1/fsym/oversampling; % sampling interval
fs    = 1/si;                % sampling frequency (two-sided)
fmax  = 1/si/2;              % half of sampling frequency = fs/2 (one-sided)

Nfft        = 2048*block_size;   % Number of two-sided points [from 0 to fs-df]
M           = Nfft/2 + 1;        % Number of one-sided points [from 0 to fmax (or fs/2)]
df          = fmax/(M-1);        % or = fs/Nfft
ff_oneSided = (0:df:fmax);
ff_twoSided = (0:df:fs-df);
tt          = (0:si:Nfft*si-si); % time span is fro 0 to 1/df = M/fmax = M*2*si

psd_itpl = interp1(freq(1:i_noise_req), psd(1:i_noise_req), ff_oneSided, 'linear', 'extrap');
psd_itpl(1) = psd(1);

% Make negative values equal to zero (it's required for he second method in which we take sqrt of psd)   
psd_itpl(psd_itpl<0) = 0;

% Plot PSD and interpolated PSD
figure(12) 
clf;
semilogx(freq(1:i_noise_req)/1e9, psd(1:i_noise_req), 'LineWidth', 4, 'disp', 'one-sided PSD') 
hold on; 
grid on; 
semilogx(ff_oneSided/1e9, psd_itpl, 'o--r', 'LineWidth', 4, 'disp', 'one-sided PSD (interpolated)')
xlabel('frequency (GHz)')

%% generate time samples from psd (method 1)

% get the one-side amplitude response from psd
asd_oneSided = sqrt(psd_itpl*length(psd_itpl)*fs);

% Generate some random phases uniformly distributed from 0 to 2*pi
rand_phase = 2*pi*rand(size(asd_oneSided));

% Form a complex signal (one-sided) in frequency domain from amplitude and phase above
asd_with_phase = asd_oneSided .* exp(sqrt(-1)*rand_phase);

% add the complex cnjugate to the complex signal to make it two-sided
asd_twoSided = [
    real(asd_with_phase(1)),... 
    1* asd_with_phase(2:+1:end-1),... 
    real(asd_with_phase(end)),...
    1* conj(asd_with_phase(end-1:-1:2))
];

sig_freq = asd_twoSided;

% Generate random time series by taking ifft of the frequency response (which should give you real values due to frequency response symmetry)
sig_time = ifft(sig_freq, Nfft); 

%% Checking the result of method 1

% calculate the psd by using fft and averaging

% number of sets for averaging periodogram
N_set = floor(length(sig_time)/block_size);

sum_PG = 0;
for i=1:N_set
    Y = fft(sig_time(1+(i-1)*block_size:i*block_size));
    Y = Y(1:block_size/2+1);
    PG = 1/length(Y)/fs * abs(Y).^2;
    PG(2:end-1) = 2*PG(2:end-1);
    sum_PG = sum_PG + PG;
end
%averaging
psd_hat = sum_PG / N_set;

ff = (0:fs/block_size:fs-fs/block_size);

% get the psd of the generated signal using matlab priodogram function
[psd_hat_os,f1]=periodogram(sig_time,boxcar(length(sig_time)),block_size,fs,'onesided');
[psd_hat_ts,f2]=periodogram(sig_time,boxcar(length(sig_time)),block_size,fs,'twosided');

figure(12) 
semilogx(ff(1:block_size/2+1)/1e9, (psd_hat(1:block_size/2+1)), 'k', 'LineWidth', 4, 'disp', 'one-sided PSD estimation usig fft')
semilogx(f1/1e9, psd_hat_os, 'c', 'LineWidth', 4, 'disp', 'one-sided PSD estimation usig periodogram')
semilogx(f2/1e9, psd_hat_ts, 'm', 'LineWidth', 4, 'disp', 'two-sided PSD estimation usig periodogram')
legend('show', 'Location', 'NorthWest')
title('psd vs estimated psd (method 1)')
set(gca, 'fontsize', 18)

%% generate time samples from psd (Method 2: using white noise)

% Generate double-side band PSD
psd_dsb = [
    real(psd_itpl(1))         ,... 
    1/2* psd_itpl(2:+1:end-1) ,... 
    real(psd_itpl(end))       ,...
    1/2*(psd_itpl(end-1:-1:2))
];

% Generte a white noise in time domain
wn_time = randn(size(psd_dsb));

% Get the frequency response of the white noise
wn_freq = fft(wn_time);

wn_filtered = wn_freq .* sqrt(psd_dsb);

time_series = ifft(wn_filtered);

% get the psd with the generated signal
[psd_est_os,f1]=periodogram(time_series,boxcar(length(time_series)),block_size,fs,'onesided');
[psd_est_ts,f2]=periodogram(time_series,boxcar(length(time_series)),block_size,fs,'twosided');


figure(13) 
clf;
semilogx(freq(1:i_noise_req)/1e9, psd(1:i_noise_req), 'LineWidth', 4, 'disp', 'one-sided PSD') 
hold on; 
grid on; 
semilogx(ff_oneSided/1e9, psd_itpl, 'o--r', 'LineWidth', 4, 'disp', 'one-sided PSD (interpolated)')
semilogx(f1/1e9, psd_est_os, 'c', 'LineWidth', 4, 'disp', 'one-sided PSD estimation usig periodogram')
semilogx(f2/1e9, psd_est_ts, 'm', 'LineWidth', 4, 'disp', 'two-sided PSD estimation usig periodogram')
xlabel('frequency (GHz)')
legend('show', 'Location', 'NorthWest')
title('psd vs estimated psd (method 1)')
set(gca, 'fontsize', 18)

Any help or feedback is appreciated from anyone who knows about this. Thanks

Answer 1

You can use the fact that the Periodogram is calculated using the conjugate square of the Fourier transform to back out a time series from any PSD. Let $\mathcal{F}_x(\omega)$ be the Fourier transform of $x$. Because the Fourier transform is complex, $$\mathcal{F}_x(\omega)=|\mathcal{F}_x(\omega)|e^{i\phi(\omega)}.$$ If we further constrain $x$ to be real, then we can say $$\phi(\omega)=-\phi(-\omega).$$ Note that because $P$ is energy, it does not contain any phase information. By defintion,

$$P(\omega)=\mathcal{F_x(\omega)}\mathcal{F}^*_x(\omega)=|\mathcal{F}_x(\omega)|e^{i\phi(\omega)}|\mathcal{F}^*_x(\omega)|e^{-i\phi(\omega)}.$$

If we solve for $\mathcal{F(x)}$, we must assume some phase distribution. In the code below, I used a uniform random sampling. Whatever phase you choose, this will not affect the resulting time series' PSD. $$\mathcal{F}_x(\omega)=\sqrt{P(\omega)}e^{i\phi(\omega)}$$ $$x=\mathcal{F}^{-1}\left(\sqrt{P(\omega)}e^{i\phi(\omega)}\right)$$ Here is my code, which accounts for the special scalings that Python (and most other languages) put on PSD's and FFT's. The first portion demonstrates calculating a PSD manually and verifies that this matches the built-in Periodogram method. The second portion uses the above derivation and undoes the manual Periodogram calculation to retrieve a new, different time series that has a PSD that matches the first generated signal.

import numpy as np
import scipy.signal
import scipy.fft

N_gen = 8192

# Generate white noise and filter it to give it an interesting PSD
w = np.random.normal(0, 1, (N_gen,))
x = scipy.signal.filtfilt(np.array([1, -2, -3, -4, 5]), np.array([1]), w) + 2

# Find PSD using built-in function and manually:
# Built-in function
f, P = scipy.signal.periodogram(x, 1, 'boxcar', nfft=8192, detrend=False)

# Manually
N = len(x)
X = scipy.fft.fft(x) / np.sqrt(8192) * np.sqrt(2)
P_fft = np.abs(X*np.conj(X))
P_fft_one_sided = P_fft[0:int(N/2)+1]  # P_fft_one_sided is identical to P
P_fft_one_sided[0] = P_fft_one_sided[0] / 2
P_fft_one_sided[-1] = P_fft_one_sided[-1] / 2
print(np.allclose(P, P_fft_one_sided))  # True, P matches P_fft_one_sided


# Now undo the manual operation of P_fft_one_sided to get back to a time series, insert your PSD here
N_P = len(P_fft_one_sided)  # Length of PSD
N = 2*(N_P - 1)

# Because P includes both DC and Nyquist (N/2+1), P_fft must have 2*(N_P-1) elements
P_fft_one_sided[0] = P_fft_one_sided[0] * 2
P_fft_one_sided[-1] = P_fft_one_sided[-1] * 2
P_fft_new = np.zeros((N,), dtype=complex)
P_fft_new[0:int(N/2)+1] = P_fft_one_sided
P_fft_new[int(N/2)+1:] = P_fft_one_sided[-2:0:-1]

X_new = np.sqrt(P_fft_new)

# Create random phases for all FFT terms other than DC and Nyquist
phases = np.random.uniform(0, 2*np.pi, (int(N/2),))

# Ensure X_new has complex conjugate symmetry
X_new[1:int(N/2)+1] = X_new[1:int(N/2)+1] * np.exp(2j*phases)
X_new[int(N/2):] = X_new[int(N/2):] * np.exp(-2j*phases[::-1])
X_new = X_new * np.sqrt(8192) / np.sqrt(2)

# This is the new time series with a given PSD
x_new = np.real(scipy.fft.ifft(X_new))

# Verify that P matches P_new
f_new, P_new = scipy.signal.periodogram(x_new, 1, 'boxcar', nfft=8192, detrend=False)
print(np.allclose(P, P_new))  # True, P matches P_new