# Random gaussian noise estimation from a given trajectory

I wish to estimate the noise of given trajectory from a measurement (data which is provided).

I consider that the noise is random gaussian noise.

I add here examples of two trajectories, one of the original data, and the other smoothened.

The trajectories are the measured velocity of some particles.

The model which is used to calculate the trajectories is quite complicated, but in the simplest form it is a langevin equation.

The results I would like to obtain at the end to have an estimation of the amplitude of the noise (considering it is an uncorrelated random gaussian noise).

I add here a csv file here

Note that there are some NaN values in the file, and the plot is done by normalizing each time series with its average.

This is the plot code for the smoothened data in matlab

    figure
hold on
plot(t,smooth(vf,25)/mean(vf,'omitnan'),'-b');
plot(t,smooth(vb,25)/mean(vb,'omitnan'),'-r');
xlabel('time'); ylabel('mean normalized intensity');
box on
legend('vf','vb')


After reading in the data (and "correcting" the NaN values with the previous valid number):

import csv
with open('Q83666vfvb.csv', newline='') as csvfile:
next(rdr)
data = list(rdr)

idx = 0
vf = np.zeros(len(data))
vb = np.zeros(len(data))
for tup in data:
vf[idx] = vf[idx-1] if tup[1] == "NaN" else tup[1]
vb[idx] = vf[idx-1] if tup[2] == "NaN" else tup[2]
idx = idx + 1


I'd first try using a simple DC blocking filter:

$$\displaystyle y[n] = x[n] - x[n-1] + \alpha\, y[n-1]$$

from scipy import signal

alpha = 0.9

vf_dc = np.zeros(len(vf))
vb_dc = np.zeros(len(vb))
vf_dc[0] = 0
vb_dc[0] = 0
for idx in np.arange(len(vf)-1) + 1:
vf_dc[idx] = vf[idx] - vf[idx-1] + alpha*vf_dc[idx-1]
vb_dc[idx] = vb[idx] - vb[idx-1] + alpha*vb_dc[idx-1]


and then just take the standard deviation of the vf_dc and vb_dc measurements.

Comparing the two results (blue) with their originals (orange):

The next thing I'd look at doing is removing the "noise" component from the original signal:

vf_corrected = vf - vf_dc
vb_corrected = vb - vb_dc


and then run a first order smoother over them:

$$\displaystyle y[n] = \alpha\, y[n-1] + (1-\alpha) x[n]$$

alpha = 0.9

vf_corrected_smooth = np.zeros(len(vf_corrected))
vb_corrected_smooth = np.zeros(len(vb_corrected))
vf_corrected_smooth[0] = vf_corrected[0]
vb_corrected_smooth[0] = vb_corrected[0]
for idx in np.arange(len(vf)-1) + 1:
vf_corrected_smooth[idx] = alpha*vf_corrected_smooth[idx-1] + (1-alpha)*vf_corrected[idx]
vb_corrected_smooth[idx] = alpha*vb_corrected_smooth[idx-1] + (1-alpha)*vb_corrected[idx]

plt.figure(1)
plt.plot(vf_corrected_smooth)
plt.plot(vf_corrected)

plt.figure(2)
plt.plot(vb_corrected_smooth)
plt.plot(vb_corrected)


Note that the above is completely disregarding the other information you have given: that this data is generated from a system described by the Langevin equation.

If playing with the two alpha values in what I've done above doesn't give satisfactory results, the next step will be to do what Filipe Pinto says in their answer and write down the differential equations governing your measurements and derive the appropriate Kalman filter (or extended Kalman filter if the system ends up being nonlinear).

Please give us more detailed information about the model your measurements follow and we can probably derive the appropriate discrete-time Kalman filter for your problem. There is some literature already that seems to apply the Kalman filter equations to Langevin systems. Perhaps have a look at that if the simple approach I'm suggesting doesn't work for you.

• Thanks, what would be your second order approximation as a filter? in case there is also a drift and not just diffusion. Jul 3, 2022 at 14:04