Determine orientation by correlation

Workspace describing a method to extract a vector measurement from directional measurements by means of forming the cross correlation

Initialization

Off[General::spell, General::spell1]

LSolve = Last[Solve[##]] &

Needs["Calculus`VectorAnalysis`"]
Needs["Statistics`MultiDescriptiveStatistics`"]
Needs["Statistics`ContinuousDistributions`"]

PeS = PowerExpand[Simplify[#]] &
FpPeS = FixedPoint[PeS, #] &

Geometric sum (magnitude of a vector in Euclidean metrics)

\(GeomSum = Function[{v}, \@Plus @@ \((\(#\^2 &\) /@ v)\)]\)

Basic principle: Reconstructing a vector from its autocorrelation

We define a 3 component vector, 
in this particular case in sperical coordinates

Bvec = 
      CoordinatesToCartesian[{B, θ, ϕ}, Spherical]

The correlation of the vector is obtained by forming its outer product with itself

MatrixForm[Bcor = Outer[Times, Bvec, Bvec]]

The eigensystem of the correlation matrix yields

Simplify[Eigensystem[Bcor]]

two eigenvectors to the eigenvalue 0,
and one non-zero eigenvector in the direction of the original vector, with the eigenvalue being the square of the magnitude of the vector.

TableForm[
      Transpose[{Eigenvalues[Bcor], 
          MatrixForm /@ Eigenvectors[Bcor]}]]

The direction of the eigenvector becomes more apparent as we multiply the last eigenvector by Cos[θ] and compare the result to the original vector

TableForm[
      MatrixForm /@ {Simplify[Last[Eigenvectors[Bcor]] Cos[θ]], 
          Bvec}]

As an aside, the trace of the correlation matrix also yields the square of the magnitude

Tr[Bcor]
Simplify[%]

Application: Signal with random amplitude, but fixed orientation

First, we generate a "sufficiently long" scalar random sequence, 
in this particular case uniformly distributed in the interval [-1,1]

rvec = Table[Random[Real, {-1, 1}], {n, 4096}];

and generate a sequence of vectors oriented along a given direction

rndcorrel = Transpose[rvec # & /@ Bvec];

The covariance matrix of the components of the random sequence is, unsurprisingly, the matrix obtained above, 
times the variance of the random sequence

CovarianceMatrix[rndcorrel] // MatrixForm

as we verify by calculating the variance as \(TraditionalForm\`⟨x\^2⟩ - ⟨x⟩\^2\)

{Integrate[ρ, {ρ, -1, 1}], Integrate[(ρ^2), {ρ, -1, 1}]}/Integrate[1, {ρ, -1, 1}]
%[[2]] - %[[1]]^2

Diagonalization of the random sequence (with numerical values for {B, θ, ϕ} because even Mathematica will hiccup with the symbolic solution) allows us to recover the vector orientation of the random sequence

Eigensystem[
      CovarianceMatrix[
        rndcorrel /. {B -> 
              1, θ ->  .5, ϕ ->  .6}]]
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
            CovarianceMatrix[
              rndcorrel /. {B -> 
                    1, θ ->  .5, ϕ ->  .6}]]]]]

as we can verify by comparing it to the eigensystem of the correlation of the orientation vector

Eigensystem[
      Bcor /. {B -> 
            1, θ ->  .5, ϕ ->  .6}]
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
            Bcor /. {B -> 
                  1, θ ->  .5, ϕ ->  .6}]]]]

Application to uncorrelated random noise

Unbiased noise

First, we generate a sequence of random vectors, 
each component of which is uniformly distributed in the interval [-1,1] in this particular case

rvec3 = 
        Table[Table[Random[Real, {-1, 1}], {m, 3}], {n, 
            4096}];

The covariance matrix is essentially diagonal, with the diagonal elements essentially equal to the variance of the individual distributions

CovarianceMatrix[rvec3] // MatrixForm

Consequently, the eigenvalues are essentially equal

Eigensystem[CovarianceMatrix[rvec3]]
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
            CovarianceMatrix[rvec3]]]]]

Obviously, the trace of the covariance matrix is the sum of the noise variances

Tr[CovarianceMatrix[rvec3]]

Biased random noise

Now we generate biased random noise

rvec3b = 
        Table[Table[Random[Real, {0, 1}], {m, 3}], {n, 4096}];

Mean, second moment
variance
of the random sequence.

{Integrate[ρ, {ρ, 0, 1}], Integrate[(ρ^2), {ρ, 0, 1}]}/Integrate[1, {ρ, 0, 1}]
%[[2]] - %[[1]]^2

Again, the covariance matrix is essentially diagonal

CovarianceMatrix[rvec3b] // MatrixForm

Since we're calculating the covariance (not the correlation), the bias is removed and we still get eigenvalues which are essentially equal to the variance of the component random sequences.

Eigensystem[CovarianceMatrix[rvec3b]]
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
            CovarianceMatrix[rvec3b]]]]]

N[(1/12)]

However, the bias does show up in the correlation matrix.

While nothing changes between correlation and covariance matrix for zero mean random noise

\(MatrixForm[\(\((Plus @@ #\/\(Length[#] - 1\))\) &\)[\(Outer[
              Times, #, #] &\) /@ rvec3]]\)
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], 
                MatrixForm /@ Eigenvectors[#]} &[%]]]]

the bias of the random sequences sticks out like a signal: 
If all the variances are equal, one eigenvector is equal to \(TraditionalForm\`\((∑⟨x\_n⟩\^2)\) + ⟨x\^2⟩\), the other ones are ⟨(x^2)⟩

\(MatrixForm[\(\((Plus @@ #\/\(Length[#] - 1\))\) &\)[\(Outer[
              Times, #, #] &\) /@ rvec3b]]\)
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], 
                MatrixForm /@ Eigenvectors[#]} &[%]]]]

\(\((∑⟨x\_n⟩\^2)\) + ⟨x\^2⟩\)

(3/4) + (1/12)
N[%]

Another case for which the eigenvectors of the correlation matrix are predictable is the case of uncorrelated, zero mean, orthogonal, noise sources with different variances

\(MatrixForm[\(\((Plus @@ #\/\(Length[#] - 1\))\) &\)[\(Outer[
              Times, #, #] &\) /@ 
          Table[Table[Random[Real, {\(-m\), m}], {m, 3}], {n, 
              16384}]]]\)
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], 
                MatrixForm /@ Eigenvectors[#]} &[%]]]]

Table[Variance[UniformDistribution[-m, m]], {m, 3}]

Application to a noisy measurement

We add zero mean, isotropic noise to the oriented random signal

CovarianceMatrix[(rndcorrel /. {B -> 
                  1, θ ->  .5, ϕ ->  .6}) + rvec3] // 
      MatrixForm
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
            CovarianceMatrix[(rndcorrel /. {B -> 
                        1, θ ->  .5, ϕ ->  .6}) + 
                rvec3]]]]]

As the noise increases, the eigenvector for the largest eigenvalue deviates increasingly from the direction of the oriented signal. 
Note the sign ambiguity of the eigenvector.

Table[
      TableForm[
        Transpose[
          Chop[{Eigenvalues[#], MatrixForm /@ Eigenvectors[#]} &[
              CovarianceMatrix[(rndcorrel /. {B -> 
                          1, θ ->  .5, ϕ ->  .6}) + ν rvec3]]]]], {ν, 0, 4,  .5}]

The deviation of the noisy measurement from the noiseless measurement can be obtained by forming the cross product of the noiseless and the noisy eigenvector. 
This formulation exploits the feature of Mathematica to normalize the eigenvectors to unit length.

Table[
      First[Eigenvectors[
          CovarianceMatrix[(rndcorrel /. {B -> 
                      1, θ ->  .5, ϕ ->  .6}) + ν rvec3]]], {ν, 0, 20,  .5}]
MatrixForm /@ (aux = (Cross[#, 
                  First[%]] & /@ %))

The magnitude of the cross product is a good indicator for the direction error

GeomSum /@ aux

ListPlot[Transpose[{Range[0, 20,  .5], GeomSum /@ aux}], 
        Frame -> True];

Sensitivity to systematic errors

Gain mismatch

Perturb the input signal by different gain factors aξ,  aυ,  aζ.
Show the eigensystem of the correlation matrix of the perturbed vector
Express the error vector as difference of the perturbed and unperturbed signal vectors
Calculate the magnitude of the error vector
Give an upper limit for the magnitude of the error as function of the maximum deviation ϵ of the gain factors aξ,  aυ,  aζ from unity.

Bvec {aξ, aυ, aζ}
Simplify[Outer[Times, %, %]];
TableForm[
      Transpose[
        Chop[{Eigenvalues[#], 
                MatrixForm /@ 
                  Eigenvectors[#]} &[%]]]]
Last[Last[Last[%]]] Cos[θ] aζ
Simplify[% - (% /. {aξ -> 1, aυ -> 1, 
              aζ -> 1})]
FpPeS[GeomSum[%]]
FpPeS[% /. {aξ -> 1 + ϵ, 
          aυ -> 1 + ϵ, 
          aζ -> 1 + ϵ}]

Phase mismatch for a sinusoidal signal

In cartesian coordinates, because spherical polar coordinates get too messy to solve

(Sqrt[2]) {x, y, z} Sin[
        2  π t + {δξ, δυ, δζ}]
MatrixForm[
      CartCor = 
        Simplify[Outer[Integrate[#1  #2, {t, 0, 1}] &, %, %]]]

eigsys = Eigensystem[CartCor];

Chop[eigsys /. {x -> 1, y -> 1, z -> 1} /. 
          Last[Solve[{δξ - δυ == ξυ, δζ - δξ == ζξ, δζ - δυ == υζ}, {ξυ, ζξ, υζ}]] /. {δξ -> (π/6.) , δυ -> 
            0, δζ -> ((-π)/6.) }]

Chop[
      Eigensystem[CartCor /. {x -> 5, y -> 2, z -> 3} /. 
            Last[Solve[{δξ - δυ == ξυ, δζ - δξ == ζξ, δζ - δυ == υζ}, {ξυ, ζξ, υζ}]] /. {δξ -> (π/12.) , δυ -> 
              0, δζ -> ((-π)/8.) }]]

Chop[
      Eigensystem[CartCor /. {x -> 5, y -> 2, z -> 3} /. 
            Last[Solve[{δξ - δυ == ξυ, δζ - δξ == ζξ, δζ - δυ == υζ}, {ξυ, ζξ, υζ}]] /. {δξ -> 
              0. , δυ -> 0. , δζ -> 
              0. }]]

CartCor /. 
      Cos[α_] -> 
        Normal[Series[Cos[α], {α, 0, 3}]]