Combine WienerProcess and WhiteNoiseProcess

Question

I'm trying to simulate an IMU (inertial measurement unit) using a simple model:

$$ \tilde a(t) = a(t) + b(t) + c(t) $$

where $\tilde a$ is the measured value, $a$ is the true value, $b$ is the random walk bias simulated as Wiener process and $c$ is Gaussian white noise:

$$ \dot b(t) = \sigma_b w(t)\\ c(t) = \sigma_c w(t) $$

where $w(t)$ is Gaussian white noise process with mean 0 and standard deviation 1.

This is modeled as

P=TransformedProcess[b[t]+c[t],
{b\[Distributed]WienerProcess[0,0.1],
c\[Distributed]WhiteNoiseProcess[]},t];

However, MMA complains that the discrete and continuous model cannot be mixed when I plot the data

ListPlot[RandomFunction[P, {0, 10, 0.01}]]

The processes in {[FormalP]1[Distributed]WienerProcess[0,0.1],[FormalP]2[Distributed]WhiteNoiseProcess[NormalDistribution[0,1]]} should be all continuous-time or all discrete-time.

Of course, we can generate the data separately using WienerProcess and WhiteNoiseProcess, but I'm wondering what is the correct way to do this symbolically so that I can utilize all the functions associated with the random process (RandomFunction, CovarianceFunction, etc.) to analyze the model?

Answer 1

Which part of the IMU are you trying to simulate...just the signals of ax, ay, az, gx, gy, gz? Or Maybe them fused together to get an angle of some kind? You didn't post all of your code, so I can't really see how you're combining your discrete and continuous models....

I have done a similar process when designing filters for an inverted pendulum I made once.

\[Alpha] = 0.2
\[CapitalDelta]t = 0.01

dgl = \[Phi]''[t] + g/l Sin[ \[Phi][t]] == -0.07 ArcTan[100 Derivative[1][\[Phi]][t]] - 0.4 Derivative[1][\[Phi]][t] + u[t];

nsys = NonlinearStateSpaceModel[{dgl}, {{\[Phi][t], 0}, {\[Phi]'[t], 0}}, {u[t]}, {\[Phi][t]}, t] /. {g -> 9.81, l -> 2};
ssm = StateSpaceModel[{dgl}, {{\[Phi][t], 0}, {\[Phi]'[t], 0}}, {u[t]}, {\[Phi][t]}, t] /. {g -> 9.81, l -> 2};

A noise function which will create whitenoise when set to 0.

noise[n_Integer /; n >= 1, \[Beta]_?NumericQ] := Re[InverseFourier[(RandomVariate[NormalDistribution[], {n, 2}].{1, I})*Prepend[Range[n - 1]^(\[Beta]/2), 0]]]

Generate discrete signals $\phi(t),\omega(t)$ and add white noise via noise[dataset length, 0] and $\alpha$ was used as a scaling value to reduce noise amplitude.

{sig1, sig2} = StateResponse[{ToDiscreteTimeModel[nsys, \[CapitalDelta]t]}, {-4 UnitBox[3 (t + -1.5)]}, {t, 0, 12}];
data1 = noise[Length[sig1], -0]/15;
data2 = noise[Length[sig2], -0]/15;

Add noise to discrete signals:

phi = Table[{\[CapitalDelta]t i, (sig1[[i + 1]]/(1 10^-\[Alpha]) + data1[[i + 1]])/(1 10^\[Alpha])}, {i, 0, Length[sig1] - 1}];
omega = Table[{\[CapitalDelta]t i, (sig2[[i + 1]]/(1 10^-\[Alpha]) + data2[[i + 1]])/(1 10^\[Alpha])}, {i, 0, Length[sig1] - 1}];

Plot!

ListLinePlot[{phi, omega}, PlotRange -> {{0, 12}, {-1.7, 1.7}}, ImageSize -> Large]

You can make this 'continuous' via Interpolation if needed...

angle = Interpolation[phi]
dotangle = Interpolation[omega]

This would be a noisy measured signal...though but it wouldn't be much different for your system simplying adding your drift b[t]...But this is a typical way to go about adding noise to a simulated signal.

There is also under NonlinearStateSpaceModel -> Applications -> Extended Kalman

An example of how a signal is generated, a noisy signal added, and then filtered via EKF.

Answer 2

You can create a TransformedProcess which time-differences a"Continuous" Wiener Process. Remarkably this matches the Discrete WhiteNoiseProcess exactly!!!

\[CapitalDelta]\[ScriptCapitalW]\[ScriptCapitalP][\[Sigma]_, \
\[Delta]_, k_] := 
 TransformedProcess[
  k (Subscript[B, t][t + \[Delta]] - Subscript[B, t][t]), 
  Subscript[B, t] \[Distributed] WienerProcess[0, \[Sigma]], t]
\[Sigma]\[Sigma] = 1; \[Delta]\[Delta] = 1/100; kk = 1;
SeedRandom[2]; data1 = 
 RandomFunction[\[CapitalDelta]\[ScriptCapitalW]\[ScriptCapitalP][\
\[Sigma]\[Sigma], \[Delta]\[Delta], kk], {0, 10, \[Delta]\[Delta]}, 1];
SeedRandom[2]; data = 
 RandomFunction[
  WhiteNoiseProcess[kk Sqrt[\[Delta]\[Delta]] \[Sigma]\[Sigma]], {0, 
   10/\[Delta]\[Delta]}];
Max[Abs[data["ValueList"][[1]] - data1["ValueList"][[1]]]]

2.22045*10^-16