The following was my original class for computing autocorrelations.
This is written in C# 4.7.0.
using System;
using System.Collections.Generic;
using System.Linq;
namespace MyNameSpace
{
public sealed class TimeSeriesAnalysis
{
private readonly IEnumerable<Vector3> _source;
private int _count;
public TimeSeriesAnalysis(IEnumerable<Vector3> vectorList)
{
_source = vectorList;
_count = vectorList.Count();
}
public double AutoCorr(int lag)
{
int n = _count;
if (n == 0 || lag >= n || lag < 0)
throw new ArgumentException($"Lag must be between 0 and {n - 1}, and the list cannot be empty.");
double sum = 0;
int targetIndex = n - lag;
for (int i = 0; i < targetIndex; i++)
{
var vecA = _source.Skip(i).Take(1).First();
var vecB = _source.Skip(i + lag).Take(1).First();
sum += vecA.Dot(vecB);
}
return sum / (n - lag);
}
public IEnumerable<double> GetCResults(int maxLag, double c0)
{
for (int lag = 0; lag <= maxLag; lag++)
{
double cValue = 0;
try
{
cValue = AutoCorr(lag);
}
catch (ArgumentException ex)
{
Console.WriteLine(ex.Message);
yield break; // Exit the loop if an invalid lag is encountered
}
yield return cValue / c0; // Normalization is done here
}
}
public PairReturn<double> GetAutoCorrelationPoints(int maxLag)
{
int vectorCount = _count;
double c0 = AutoCorr(0); // This is the normalization factor
Console.WriteLine($"Normalization factor: {c0}");
int validMaxLag = Math.Min(maxLag, _count - 1); // Ensure maxLag does not exceed n-1
var values = Enumerable.Range(0, validMaxLag + 1).Select(lag => (double)lag);
return new PairReturn<double>(values, GetCResults(maxLag, c0));
}
}
}
This class works well.
using System;
using System.Collections.Generic;
using Microsoft.VisualStudio.TestTools.UnitTesting;
using MyNameSpace;
namespace MyNameSpace____unit_test
{
[TestClass]
public class TimeSeriesAnalysisVec3AutoCorrelationMemoryLess____unit_test
{
[TestInitialize]
public void Initialize()
{
}
[TestMethod]
public void GetAutoCorrelationTestMethod()
{
IList<Vector3> vect3List = new List<Vector3>()
{
new Vector3(1,1,1),
new Vector3(2,2,2),
new Vector3(3,3,3),
new Vector3(4,4,4),
new Vector3(5,5,5)
};
var timeSeriesAnalysis = new TimeSeriesAnalysis(vect3List);
double c0 = timeSeriesAnalysis.AutoCorr(0);//33
double c1 = timeSeriesAnalysis.AutoCorr(1);//30
double c2 = timeSeriesAnalysis.AutoCorr(2);//26
double c3 = timeSeriesAnalysis.AutoCorr(3);//21
double c4 = timeSeriesAnalysis.AutoCorr(4);//15
Assert.AreEqual(33, c0);
Assert.AreEqual(30, c1);
Assert.AreEqual(26, c2);
Assert.AreEqual(21, c3);
Assert.AreEqual(15, c4);
}
[TestMethod]
public void GetAutoCorrelationPointsTestMethod()
{
List<Vector3> vect3List = new List<Vector3>()
{
new Vector3(1,1,1),
new Vector3(2,2,2),
new Vector3(3,3,3),
new Vector3(4,4,4),
new Vector3(5,5,5)
};
var timeSeriesAnalysis = new TimeSeriesAnalysis(vect3List);
var autocorr = timeSeriesAnalysis.GetAutoCorrelationPoints(10);
IList<double> tausList = new List<double>(autocorr.XData);
IList<double> autocorrList = new List<double>(autocorr.YData);
Assert.AreEqual(5, tausList.Count);
Assert.AreEqual(5, autocorrList.Count);
Assert.AreEqual(1.0, autocorrList[0]);
Assert.AreEqual(0.90909090909090906, autocorrList[1]);
Assert.AreEqual(0.78787878787878785, autocorrList[2]);
Assert.AreEqual(0.63636363636363635, autocorrList[3]);
Assert.AreEqual(0.45454545454545453, autocorrList[4]);
}
}
}
Then, I had to rewrite this class using FFT (fast fourier transform).
using System;
using System.Collections.Generic;
using System.Linq;
using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
namespace MyNameSpace
{
public sealed class TimeSeriesAnalysisFFT
{
private readonly List<Vector3> _source;
private int _count;
public TimeSeriesAnalysisFFT(IEnumerable<Vector3> vectorList)
{
_source = vectorList.ToList();
_count = _source.Count;
}
public double AutoCorr(int lag)
{
int n = _count;
if (n == 0 || lag >= n || lag < 0)
throw new ArgumentException($"Lag must be between 0 and {n - 1}, and the list cannot be empty.");
// Using FFT to compute autocorrelation
var realValues = _source.Select(v => v.X).ToArray();
// Remove mean to center the data
double mean = realValues.Average();
for (int i = 0; i < n; i++)
{
realValues[i] -= mean;
}
var complexValues = new Complex32[n];
for (int i = 0; i < n; i++)
{
complexValues[i] = new Complex32((float)realValues[i], 0);
}
// Forward FFT
Fourier.Forward(complexValues, FourierOptions.Matlab);
// Compute power spectrum (magnitude squared)
for (int i = 0; i < n; i++)
{
complexValues[i] = complexValues[i] * Complex32.Conjugate(complexValues[i]);
}
// Inverse FFT to get the autocorrelation
Fourier.Inverse(complexValues, FourierOptions.Matlab);
// Extract the real part and normalize by number of samples
var autocorr = complexValues.Select(c => c.Real / n).ToArray();
// Return the normalized value for the given lag
double variance = realValues.Sum(x => x * x) / n;
if (variance == 0)
return 0;
return autocorr[lag] / variance;
}
public IEnumerable<double> GetCResults(int maxLag, double c0)
{
for (int lag = 0; lag <= maxLag; lag++)
{
double cValue = 0;
try
{
cValue = AutoCorr(lag);
}
catch (ArgumentException ex)
{
Console.WriteLine(ex.Message);
yield break; // Exit the loop if an invalid lag is encountered
}
yield return cValue / c0; // Normalization is done here
}
}
public PairReturn<double> GetAutoCorrelationPoints(int maxLag)
{
int vectorCount = _count;
double c0 = AutoCorr(0); // This is the normalization factor
Console.WriteLine($"Normalization factor: {c0}");
int validMaxLag = Math.Min(maxLag, _count - 1); // Ensure maxLag does not exceed n-1
var values = Enumerable.Range(0, validMaxLag + 1).Select(lag => (double)lag);
return new PairReturn<double>(values, GetCResults(maxLag, c0));
}
}
}
However, the following unit-test fails:
using System;
using System.Collections.Generic;
using Microsoft.VisualStudio.TestTools.UnitTesting;
using MyNameSpace;
namespace MyNameSpace____unit_test
{
[TestClass]
public class TimeSeriesAnalysisFFT____unit_test
{
[TestInitialize]
public void Initialize()
{
}
[TestMethod]
public void GetAutoCorrelationTestMethod()
{
IList<Vector3> vect3List = new List<Vector3>()
{
new Vector3(1,1,1),
new Vector3(2,2,2),
new Vector3(3,3,3),
new Vector3(4,4,4),
new Vector3(5,5,5)
};
var timeSeriesAnalysis = new TimeSeriesAnalysisFFT(vect3List);
double c0 = timeSeriesAnalysis.AutoCorr(0);//33
double c1 = timeSeriesAnalysis.AutoCorr(1);//30
double c2 = timeSeriesAnalysis.AutoCorr(2);//26
double c3 = timeSeriesAnalysis.AutoCorr(3);//21
double c4 = timeSeriesAnalysis.AutoCorr(4);//15
Assert.AreEqual(33, c0);
Assert.AreEqual(30, c1);
Assert.AreEqual(26, c2);
Assert.AreEqual(21, c3);
Assert.AreEqual(15, c4);
}
[TestMethod]
public void GetAutoCorrelationPointsTestMethod()
{
List<Vector3> vect3List = new List<Vector3>()
{
new Vector3(1,1,1),
new Vector3(2,2,2),
new Vector3(3,3,3),
new Vector3(4,4,4),
new Vector3(5,5,5)
};
var timeSeriesAnalysis = new TimeSeriesAnalysisFFT(vect3List);
var autocorr = timeSeriesAnalysis.GetAutoCorrelationPoints(10);
IList<double> tausList = new List<double>(autocorr.XData);
IList<double> autocorrList = new List<double>(autocorr.YData);
Assert.AreEqual(5, tausList.Count);
Assert.AreEqual(5, autocorrList.Count);
Assert.AreEqual(1.0, autocorrList[0]);
Assert.AreEqual(0.90909090909090906, autocorrList[1]);
Assert.AreEqual(0.78787878787878785, autocorrList[2]);
Assert.AreEqual(0.63636363636363635, autocorrList[3]);
Assert.AreEqual(0.45454545454545453, autocorrList[4]);
}
}
}
What am I doing wrong?