Covariance Matrix by Simulation

19 visualizaciones (últimos 30 días)
S. David
S. David el 20 de Jul. de 2013
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

Respuestas (2)

Youssef  Khmou
Youssef Khmou el 20 de Jul. de 2013
Editada: Youssef Khmou el 20 de Jul. de 2013
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 .
  2 comentarios
S. David
S. David el 21 de Jul. de 2013
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
Youssef Khmou el 21 de Jul. de 2013
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))

Iniciar sesión para comentar.


Wayne King
Wayne King el 20 de Jul. de 2013
Editada: Wayne King el 20 de Jul. de 2013
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.
  2 comentarios
S. David
S. David el 20 de Jul. de 2013
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
Wayne King el 20 de Jul. de 2013
Editada: Wayne King el 20 de Jul. de 2013
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
For example how to detect whether a sequence is white noise using autocorrelation.

Iniciar sesión para comentar.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by