3
$\begingroup$

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);
$\endgroup$
  • $\begingroup$ Shall you please tell me when we need to create PSD from auto correlation? What's the application $\endgroup$ – FabioSpaghetti Aug 15 at 16:07
  • $\begingroup$ 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 $\endgroup$ – OscarNieves Aug 15 at 21:40
2
$\begingroup$

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:

enter image description here

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – OscarNieves Aug 15 at 21:38
  • $\begingroup$ @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. $\endgroup$ – Matt L. Aug 16 at 6:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.