# PSD from autocorrelation in MATLAB

I am trying to simulate a simple stochastic process defined by the equation: $$$$\frac{1}{v}\frac{db}{dt} +\Gamma_0 b= \sqrt{\sigma}R(t),$$$$ 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 $$$$b(t+\Delta t) = b(t) -v \Gamma_0 b(t) \Delta t + v\sqrt{\sigma} \sqrt{\Delta t} N(0,1),$$$$ 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 $$$$S[f] = \int_{-\infty}^{\infty} e^{-2\pi f\tau} \left\langle b(t)b^*(t+\tau)\right\rangle d\tau,$$$$ from the autocorrelation function $$\left\langle b(t)b^*(t+\tau)\right\rangle$$. Analytically I can compute this and it gives me $$$$S[f] = \frac{v\sigma}{2\Gamma_0}\frac{2v\Gamma_0}{(v\Gamma_0)^2 + (2\pi f)^2},$$$$ 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);

• Shall you please tell me when we need to create PSD from auto correlation? What's the application – FabioSpaghetti Aug 15 at 16:07
• I am modelling stochastic processes (e.g. thermal noise built up inside of a waveguide) and the correct way to do it from the differential equations is to first compute the auto-correlation function, and then calculate the PSD from that to get the frequency response of the noise. Unfortunately I haven't come across much literature on how to numerically do this – OscarNieves Aug 15 at 21:40

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:

• The term $\sqrt(\Delta t)$ term comes from $R(t)$ representing the derivative of a Wiener process (which has variance $\sqrt(\Delta t)$ and mean zero). This is simulated in the Euler-Mayurama method by multiplying $\sqrt(\Delta t)$ by a standard normal random number during each step of the iterative scheme. Also if you change $\Delta t$, you can see that the curves still do not match – OscarNieves Aug 15 at 21:38
• @OscarNieves: $\sqrt{\Delta t}$ is not the variance but the standard deviation, so $\Delta t$ is the variance. With these corrections (correct computation of mY and use of $\Delta t$ instead of its square root) the curves match, as shown in my plot. – Matt L. Aug 16 at 6:55