Covariance Matrix by Simulation

Hello,

I want to check the covariance matrix of a colored noise if it is correct. I started to see how to find the covariance matrix of a white noise first. So, I did the following:

 SNRdB=10;
 SNR=10.^(SNRdB./10);

 N=1000;
 Noise=zeros(100,N);
 for ii=1:100
    n=(1/sqrt(2*SNR))*(randn(N,1)+1i.*randn(N,1));
    Noise(ii,:)=n;
 end

 NoiseCov=cov(Noise);

Theoretically, NoiseCov should be a diagonal matrix with its diagonal elements be 1/SNR. But this is not the case. Why? Did I do it wrongly?

Thanks

Answer 1

[Youssef Khmou]

hi S. David

What you described is correct, and as Wayne King said, one needs to to perform N times and averages them to get the result,

Try this version for your experience :

 SNR=10;
sigma=10^(-SNR/20);
N=400;
r=sigma*randn(N);
V=cov(r);
figure, imagesc(V);
mean(diag(V))

You get 1/SNR=0.1 and E[diag(Cov(r))]=0.0991 .

[S. David]

Are you saying that I need to find cov(Noise) (in my code) several times (say L times), and then average them by L to get the covariance matrix?

[Youssef Khmou]

David, yes, you see the result above the diag was 0.0991, now run 100 times and average the result you get exactly 1/SNR :

 SNR=10;
 sigma=10^(-SNR/20);
 N=400;
 L=100; %  Number of runs 
 S=0;
 for n=1:L

r=sigma*randn(N);
V=cov(r);
S=S+V/L;
end

 figure, imagesc(S);
 mean(diag(S))

Answer 2

[Wayne King]

You have to consider that you are looking at samples from a sampling distribution of a statistic, the sample covariance.

Consider the following.

I'll set the random number generator to its default settings for reproducible results.

rng default;
X = randn(100,2);
plot(X(:,1),X(:,2),'.')

You can see in the scatter plot that there is no obvious relationship between the two columns. The scatter plot is very much what you expect for uncorrelated random variables. Now, execute:

Y = cov(X);

The off diagonal covariance is 0.0881. The variances of the two random variables are 1.3512 and 1.010, close to the expected value of 1, but not exact.

If you have the Statistics Toolbox, you can use mvnrnd() to model correlation.

    mu = [0 0]; Sigma = [1 .6; .6 1];
    X = mvnrnd(mu,Sigma,100);
    plot(X(:,1),X(:,2),'.')

The scatter plot shown above shows a clear linear relationship between the two random variables.

   Y = cov(X);

Now the off diagonal elements of the matrix are close to 0.6.

Finally, to look at a histogram indicative of the sampling distribution for the uncorrelated case, let's repeat the experiment 500 times.

    for nn = 1:500
      X = randn(100,2);
      Y = cov(X);
      cov_values(nn) = Y(1,2);
   end 
   hist(cov_values)
   mean(cov_values)

You should see that the sampling distribution looks Gaussian with a mean that is very close to the theoretical value of 0.

So bottom line, for any given realization you cannot expect the off diagonal covariance to be zero, but in repeated sampling, the statistic (sample covariance) will have zero mean. That is always the case when dealing with statistics, you have to take into account the sampling distribution.

[S. David]

Thanks for your reply. So, how to find the covariance matrix of an AWGN vector of size N-by-1, then? I need the simulated result to compare it with I have theoretically to make sure I have derived the noise covariance matrix correctly.

[Wayne King]

The covariance matrix implies that you have a bivariate sample, not a univariate sample. If you have a random vector, then cov() will just give you an estimate of the variance.

You can compute the autocovariance sequence. This gives you the covariance between lagged values of the random vector.

If you have the Signal Processing Toolbox, you can compute the cross-covariance with xcov(), or the cross correlation with xcorr(), see

https://www.mathworks.com/help/signal/ug/confidence-intervals-for-sample-autocorrelation.html

For example how to detect whether a sequence is white noise using autocorrelation.