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