PSD from autocorrelation in MATLAB

Question

I am trying to simulate a simple stochastic process defined by the equation: \begin{equation} \frac{1}{v}\frac{db}{dt} +\Gamma_0 b= \sqrt{\sigma}R(t), \end{equation} where $R(t)$ is a zero-mean white noise (delta-correlated in time), and $v$, $\Gamma_0$ and $\sigma$ are all real positive parameters. I am simulating this using an Euler-mayurama scheme, namely \begin{equation} b(t+\Delta t) = b(t) -v \Gamma_0 b(t) \Delta t + v\sqrt{\sigma} \sqrt{\Delta t} N(0,1), \end{equation} with $N(0,1)$ being a standard normal random number, for a finite interval $t\in [0,T]$. I am trying to compute the power spectral density using the Wiener-Khintchine theorem \begin{equation} S[f] = \int_{-\infty}^{\infty} e^{-2\pi f\tau} \left\langle b(t)b^*(t+\tau)\right\rangle d\tau, \end{equation} from the autocorrelation function $\left\langle b(t)b^*(t+\tau)\right\rangle$. Analytically I can compute this and it gives me \begin{equation} S[f] = \frac{v\sigma}{2\Gamma_0}\frac{2v\Gamma_0}{(v\Gamma_0)^2 + (2\pi f)^2}, \end{equation} But when I try to compute the same numerically, my function looks completely different. Am I implementing this correctly? (see Matlab code below)

%% Autocorrelation test
close all; clear all; clc;

% Inputs
dt = 1e-3;
T = 2;
tv = 0:dt:T;
Lt = length(tv);

v = 3.5;
Gamma0 = 20.3;
sigma = 0.75;
a = v*Gamma0;

% Frequency spectrum properties
Fs = 1/dt; % (Hz) sampling frequency, based on Ben's paper

sims = 1; % number of simulations
U_Walks = zeros(Lt,sims);
for nn = 1:sims
    b = zeros(Lt,1);
    for ii = 1:Lt-1
        b(ii+1) = b(ii) - v*Gamma0*b(ii)*dt + ...
            v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);
    end
    Rxx(:,nn) = xcorr(b);
    U_Walks(:,nn) = b;
    disp(nn);
end

%% Stats
LR = length(Rxx);
N = 2^nextpow2(LR);
Y = fft(Rxx,N); % Taking only one side of the y-axis symmetric FFT
Y = Y(1:N/2+1);
mY = (Y).^2/N;
F = Fs*(0:(N/2))/N; % calculating the frequency range for the x-axis

S_f = 10*log10(mY);

S_f_analytic = 10*log10((v*sigma/2/Gamma0.*(2*a./(a^2 + (2*pi*F).^2))));

%% PLOTS
close all;

figure('units','normalized','outerposition',[0 0 1 1])
set(gcf,'Color','w');
plot(F,S_f,'r',F,S_f_analytic,'b','LineWidth',3);
legend('Numeric','Analytic');
title(['PSD']);
xlabel('f (Hz)'); ylabel('S[f] (dB)');
set(gca,'FontSize',20);

Answer 1

There are two mistakes in your code/method. The first is the term $\sqrt{\Delta t}$ in your second formula; it should be replaced by $\Delta t$. The second is in the computation of the power spectrum from the estimated auto-correlation. What you do is square the result of the FFT Y to obtain mY, but that's not correct. First of all, Y is complex-valued, and second, the Fourier transform of the auto-correlation is the power spectrum, so there's no need to square the outcome of the FFT, you just need to compute its magnitude.

In sum, replace the line

v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);

by

v*sqrt(sigma)*dt.*normrnd(0,1);

and

mY = (Y).^2/N;

by

mY = abs(Y)/N;

This should give you the desired result, as shown in the figure below: