Gauss-Markov process generation: how to filter a white noise sequence properly

Hi all,

I'm trying to generate a first order Gauss-Markov process (https://en.wikipedia.org/wiki/Gauss%E2%80%93Markov_process).

From signal processing theory, it is known that this could be performed filtering a white gaussian noise properly (noise shaping filter).

This filter should have a transfer function equal to:

with the two parameters %sigma, %beta that have to be tuned and that will described completely the Gauss-Markov process.

fs = 100;
Nsamples = 2000;

rng('default');
gaussianNoise = randn(Nsample, 1);
whiteGaussianNoise = (gaussianNoise - mean(gaussianNoise)) / std(gaussianNoise);

beta = 0.01;
sigma = 0.1;

b = sqrt(2 * beta * sigma^2);
a = [1, beta];

gaussMarkov = filter(b,a,whiteGaussianNoise);

This is the code I used for the Gauss Markov synthetis. Unfortunately, something is wrong because it has a similar power spectral density and a similiar autocorrelation of the white noise. In the link provided above, the correct auocorrelation and psd are defined.

What's wrong with this code? Why I cannot generate a Gauss Markov process?

Thank you all.

Answer 1

[Adil]

I'm not sure if that's the reason for your problem, but you're using a transfer function in the s-domain G(s) and it looks like the MATLAB filter function uses a transfer function in the z-domain G(z).

So maybe it will work after applying this method: Matched Z-transform method - Wikipedia

EDIT:

Ok, now I'm sure, just watch this video: https://www.youtube.com/watch?v=88tWmyBaKIQ and you will see that you can just use the bilinear transform and get (for your sampling frequency) the following parameters:

b = [((sqrt(2 * beta * sigma^2))*1/fs)/(beta*1/fs + 2), (sqrt(2 * beta * sigma^2)*1/fs)/(beta*1/fs + 2)];

a = [1,(beta*1/fs - 2)/(beta*1/fs + 2)];

[Johnny Scalera]

Thank you @Adil for your answer. I am taking into account your observation.

Using the bilinear transform described in the section 9.2 Converting S Domain to Z Domain of this page

https://www.sciencedirect.com/topics/engineering/domain-equation

I obtain this transfer function in the z domain

So my two vectors of coefficients are

k = sqrt(2*beta*sigma^2) / (beta+1)

b = [k, k]
a = [1, (beta-1) / (beta+1)]

Am I right? Why do you use different coefficients?

[Adil]

Hi Johnny,

I'm not sure why your source uses a different bilinear transform without dependency of dt, but as you suggested in your entry post: You can check the generated signal via its autocorrelation function, psd or Allan deviation plot.

So if I'm checking the gauss markov signal which was created with your parameters, then I can unfortunately see no exponentially time correlated autocorrelation plot, as expected.

But with the parameters which I posted above, I can see it.

[Rxx, lags] = xcorr(gaussMarkov);
figure();
plot(lags, Rxx);

This shape is what you would expect: https://www.researchgate.net/profile/Spiros_Pagiatakis/publication/232711015/figure/download/fig14/AS:668420496453664@1536375273795/Autocorrelation-function-of-the-first-order-Gauss-Markov-process.png

I calculated the coeffitients as described here: Bilinear transform - Wikipedia with s = (2 - 2z^(-1))/(dt + dtz^(-1))

I hope this helps