# Estimating Correlation Matrices

I am trying to obtain the correlation matrices of two random signals. Both of them, $$X$$ and $$Y$$, are white Gaussian Noise, with unitary variance. However, they are correlated, with correlation coefficient $$\rho$$. I generate $$Y$$ from $$X$$ using the following expression $$Y = \rho X + \sqrt{1 - \rho^2}Q,$$ where, $$Q$$ is also white Gaussian noise with unitary variance ($$\sigma^2 = 1$$). The correlation matrices I compute as $$R_{X} = X X^{H}/N$$, where $$N$$ is the number of realizations of $$X$$. I do the equivalent for $$R_{YY}$$.

This is my Matlab code:

nSamples = 1e3; % Number of samples in each trial of the process
nTrials = 1e4; % Number of trials
rho = 0.65; % Correlation between two stochastic process

%% Independent samples
X = randn(nSamples, nTrials); % Each trial in each column
Y = rho*X + sqrt(1-rho^2)*randn(nSamples, nTrials);

Rxx = X*X'/nTrials;
Ryy = Y*Y'/nTrials;


From this, I get, for instance:

Rxx(1:5, 1:5)

ans =

1.0031    0.0227    0.0113   -0.0025    0.0130
0.0227    0.9916    0.0005   -0.0094   -0.0030
0.0113    0.0005    1.0065    0.0002   -0.0006
-0.0025   -0.0094    0.0002    1.0093    0.0154
0.0130   -0.0030   -0.0006    0.0154    0.9979


which is what I was expecting. However, when I try to add correlation between the samples of the same signal, things go wrong. For instance, first I tried to use the filter function, with the taps $1, -r$ in the feedback filter, where $r$ is the correlation coefficient of an autoregressive process of first order. Here is the code:

r = 0.8;

% Method 1:
fw = [1 -r];
fw = fw/sum(fw);

X2 = filter(1, fw, randn(nSamples, nTrials));
Y2 = filter(1, fw, randn(nSamples, nTrials));
Y2 = rho*X2 + sqrt(1-rho^2)*Y2;

Rxx2 = X2*X2'/nTrials;
Ryy2 = Y2*Y2'/nTrials;


When I look at the results, for instance of $$R_{XX}$$, I get

Rxx2(1:5, 1:5)

ans =

0.0402    0.0323    0.0256    0.0203    0.0165
0.0323    0.0667    0.0536    0.0431    0.0348
0.0256    0.0536    0.0829    0.0664    0.0532
0.0203    0.0431    0.0664    0.0928    0.0747
0.0165    0.0348    0.0532    0.0747    0.1009


which is not what I was expecting. Although the non-diagonal values seems symmetric, the diagonal itself does not contain consistent values, as before. I understand they would not be close to 1, but that they would be close to each other.

I also tried without the filter function:

X = randn(nSamples, nTrials); % Each trial in each column
X2 = zeros(nSamples, nTrials);
X2(1, :) = X(1, :);

Y = randn(nSamples, nTrials); % Each trial in each column
Y2 = zeros(nSamples, nTrials);
Y2(1, :) = Y(1, :);

for k = 2:nSamples

X2(k, :) = r*X(k-1, :) + sqrt(1-r^2)*randn(1, nTrials);
Y2(k, :) = r*Y(k-1, :) + sqrt(1-r^2)*randn(1, nTrials);

end
Y2 = rho*X2 + sqrt(1-rho^2)*Y2;

Rxx2 = X2*X2'/nTrials;
Ryy2 = Y2*Y2'/nTrials;


which gave me

Rxx2(1:5, 1:5)

ans =

0.0000    0.0000    0.0000   -0.0000    0.0000
0.0000    0.9916   -0.0007    0.0122   -0.0074
0.0000   -0.0007    1.0157   -0.0051   -0.0104
-0.0000    0.0122   -0.0051    0.9996   -0.0003
0.0000   -0.0074   -0.0104   -0.0003    0.9928


which has an acceptable diagonal, except for the first element.

My question is: is this the right way to compute the correlation matrix? Is there anything that I am doing wrong? I am also not sure whether the way I am correlating the samples are correct or not.

Thanks in advance for the help.

Update: I noticed an error in the code using the AR process. Where it was X2(1) = X(1); I changed to X2(1, :) = X(1, :);. This was resulting in the first column and row of the correlation matrix with zeros. With this fix, I am getting

Rxx2(1:5, 1:5)

ans =

1.0167    0.8181    0.0004    0.0109    0.0019
0.8181    1.0111   -0.0044    0.0087   -0.0043
0.0004   -0.0044    0.9928    0.0029    0.0122
0.0109    0.0087    0.0029    0.9989    0.0132
0.0019   -0.0043    0.0122    0.0132    1.0224


which has the correct diagonal, and $$R_{XX}(1)$$ apparently correct. The other elements does not seem to be correct, however. I still have not figured out the problem with the filter method, though.

• $Y$ doesn't have unit variance but a variance of $1-\rho^2+\rho^4$.; the factor in front of $Q$ needs to be $\sqrt{1-\rho^2}$, not just $1-\rho^2$, then the variance of $Y$ becomes $1$. Aug 31, 2021 at 15:52
• Hi @MarcusMüller, thank you for taking the time to point that out. You are correct. I am actually using the square root in my codes. I have just forgot it in the minimal example that I wrote to show here. I will update the question, since correcting it did not change the matrices' behavior I have shown. Aug 31, 2021 at 21:49