How to generate a first order Gauss-Markov process in Python

Question

I would like to synthesize a first order Gauss-Markov process from a white Gaussian noise.

I know from signal processing theory that it could be realized using a noise shaping filter designed properly (see Gauss-Markov process).

First order Gauss-Markov processes can be described from two key parameters: $\sigma$, which is the standard deviation of the process, and the time constant $\beta$.

The shaping filter should have a transfer function equal to the one in this figure:

Here it is my code:

import scipy.signal as dsp
import numpy as np

Nsamples = 2000
fs = 100

time = np.arange(Nsamples) / fs

rng = np.random.default_rng()
gaussianNoise = rng.standard_normal(size=time.shape)

whiteGaussianNoise = (gaussianNoise - np.mean(gaussianNoise)) / np.std(gaussianNoise)
print('\n\n\nWGN MEAN: ', np.mean(wgn))
print('WGN STD: ', np.std(wgn))


beta = 0.01
sigma = 0.1
b = np.array([np.sqrt(2 * beta * sigma**2)])
a = np.array([1, beta])

gaussMarkovNoise = dsp.lfilter(b, a, whiteGaussianNoise)

Unfortunately, something is wrong because the gaussMarkovNoise should have an autocorrelation with an exponential decay (see http above); while filtered in this way, it still has a spike in the origin as a white noise sequence. What am I missing?

Answer 1

I believe the problem is that your values for $\beta$ and $\sigma$ are way too small.

If I try

$$ \beta = 0.01$$ $$\sigma = 0.1$$

then do

import matplotlib.pyplot as plt    
acf = plt.acorr(gaussMarkovNoise, maxlags = 1000)

then I get

However, if I do

$$ \beta = 1$$ $$\sigma = 1$$

then I get

which is much more like what you expect.