# Noise covariance matrix

I am trying to find the noise covariance matrix of a random vector $$\mathbf{x}$$ with itself. The elements of this vector represent image pixels and are treated as random variables i.e, $$\mathbf{x}=\{x_1,x_2 \dots x_n \}$$. For a $$M\times M$$ image we have $$n=M*M$$. The resulting covariance matrix, therefore, should be of the size $$n\times n$$.

The covariance matrix can be given as:

$$COV(\mathbf{x}) = \mathbb{E}((\mathbf{x}-\mathbb{E}(\mathbf{x}))(\mathbf{x}-\mathbb{E}(\mathbf{x}))^\dagger)$$

In element form I can write it as:

$$COV(x_i,x_j) = \mathbb{E}((x_i-\mathbb{E}(x_i))(x_j-\mathbb{E}(x_j))^*), \quad 1\leq i,j\leq n$$

With this definition I think I need to have $$k$$ images or vectors, resulting in $$k$$ samples for each pixel $$x$$, to compute the covariance matrix. I wrote an example code that generates $$50$$ images with Gaussian noise and finds the covariance matrix as defined above. However, I am not sure if this is correct or not.

x = Flatten[Table[Exp[-u^2/5 - v^2/5], {u, -10, 10, 1}, {v, -10, 10, 1}], 1];
Noise[x_] := # + RandomVariate[NormalDistribution[0.5, 0.2]] & /@ x;
NoisyDat = Table[Noise[x], {k, 1, 30}];   (* Generate a list of noisy image vectors *)
MeanDat = Mean[NoisyDat];      (* E(x) *)
OuterProd[list_] := KroneckerProduct[list - MeanDat, Conjugate[list - MeanDat]];  (* Outer product: (x- E(x))((x- E(x))*)^T *)
OuterMeanRemoved = Map[OuterProd, NoisyDat];
COV = Mean[OuterMeanRemoved];
ArrayPlot[COV]


I have $$30$$ recorded images from which I want to find the noise covariance matrix and I am first trying to make a correct code by simulating Gaussian noise. Is this the correct way to do it?.

EDIT: Using Covaiance[NoisyDat] gives different result than the code above.