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.
