diff --git a/data/numpy/XCorrNumpyData.json b/data/numpy/XCorrNumpyData.json new file mode 100644 index 000000000..afb6bfaa2 --- /dev/null +++ b/data/numpy/XCorrNumpyData.json @@ -0,0 +1,300 @@ +{ + "x": [ + -1.254598811526375, + 4.507143064099161, + 2.3199394181140507, + 0.986584841970366, + -3.439813595575635, + -3.4400547966379733, + -4.419163878318005, + 3.6617614577493516, + 1.0111501174320878, + 2.0807257779604544, + -4.7941550570419755, + 4.699098521619943, + 3.324426408004218, + -2.8766088932172384, + -3.181750327928994, + -3.165954901465662, + -1.9575775704046228, + 0.24756431632237863, + -0.6805498135788426, + -2.0877085980195806, + 1.118528947223795, + -3.6050613934795814, + -2.0785535146478185, + -1.336381567063083, + -0.43930015782964027, + 2.8517596139301364, + -3.003262178416403, + 0.14234438413611628, + 0.9241456886204249, + -4.5354958728000225, + 1.0754485190143832, + -3.294758763127085, + -4.349484070147205, + 4.488855372533333, + 4.656320330745594, + 3.0839734811646107, + -1.9538623082662934, + -4.023278859936161, + 1.8423302651215687, + -0.5984750626039865, + -3.7796176515522117, + -0.048230898887298146, + -4.656114788847816, + 4.093204020787821, + -2.412200183999831, + 1.6252228435398202, + -1.8828892391058902, + 0.20068021177810813, + 0.4671027934327965, + -3.1514554447447294, + 4.695846277645586, + 2.7513282336111455, + 4.394989415641891, + 3.9482735042764876, + 0.9789997881108512, + 4.218742350231167, + -4.115074979480805, + -3.040171375808548, + -4.54772711089462, + -1.7466966923673564, + -1.1132271031051797, + -2.286509682261041, + 3.2873750915192943, + -1.432466733064107 + ], + "y": [ + 0.9403061933191574, + 0.9476590216968434, + 0.9545475635032179, + 0.960966514091227, + 0.9669109283698445, + 0.9723762250642471, + 0.9773581906633237, + 0.9818529830509348, + 0.9858571348176298, + 0.9893675562498186, + 0.9923815379936803, + 0.9948967533913998, + 0.9969112604876068, + 0.9984235037041954, + 0.9994323151820118, + 0.9999369157881839, + 0.9999369157881839, + 0.9994323151820118, + 0.9984235037041954, + 0.9969112604876068, + 0.9948967533913998, + 0.9923815379936803, + 0.9893675562498186, + 0.9858571348176298, + 0.9818529830509348, + 0.9773581906633237, + 0.9723762250642471, + 0.9669109283698445, + 0.960966514091227, + 0.9545475635032179, + 0.9476590216968434, + 0.9403061933191574 + ], + "result_valid": [ + -20.327874953772945, + -23.31997659695665, + -23.40624531193906, + -21.213498734189034, + -19.19074781668703, + -17.7360646923138, + -18.26024981538604, + -12.405078452347947, + -16.394233974237434, + -20.870760171745225, + -22.84719683526064, + -22.71838868089418, + -23.277434240184192, + -28.59138991399194, + -24.305197395967138, + -23.030866628178117, + -19.84020035538177, + -17.539113593281897, + -20.69976356425527, + -15.628495665912101, + -11.034341067008846, + -7.873715061102328, + -0.6813686896583899, + 2.298616633098259, + 7.623253306630469, + 4.293877950846721, + -1.1313068425646469, + -2.515212580337982, + -4.256344451433963, + -6.141346277768917, + -4.03703972921134, + -1.9751209916296748, + -0.2411044036181078 + ], + "result_same": [ + -3.6835021145119216, + -5.5628652488084125, + -5.381304213775536, + -6.067887771687595, + -8.07996643394984, + -7.08953178194481, + -10.528470437193244, + -12.554983229986522, + -13.893653047127422, + -14.392343973410926, + -11.792708351730493, + -14.670400116328647, + -14.606390104782424, + -13.798818710027865, + -18.110822790074945, + -17.17319300712088, + -20.327874953772945, + -23.31997659695665, + -23.40624531193906, + -21.213498734189034, + -19.19074781668703, + -17.7360646923138, + -18.26024981538604, + -12.405078452347947, + -16.394233974237434, + -20.870760171745225, + -22.84719683526064, + -22.71838868089418, + -23.277434240184192, + -28.59138991399194, + -24.305197395967138, + -23.030866628178117, + -19.84020035538177, + -17.539113593281897, + -20.69976356425527, + -15.628495665912101, + -11.034341067008846, + -7.873715061102328, + -0.6813686896583899, + 2.298616633098259, + 7.623253306630469, + 4.293877950846721, + -1.1313068425646469, + -2.515212580337982, + -4.256344451433963, + -6.141346277768917, + -4.03703972921134, + -1.9751209916296748, + -0.2411044036181078, + 3.7865468988791626, + -0.4633514984678804, + -4.834133477008061, + -7.699995878935996, + -5.840282601335777, + -2.0632951592104236, + -3.7864088162752196, + -3.217278962065742, + 0.31510946161231823, + 0.3382931581490354, + 4.6577407261346995, + 0.7798465301709483, + 2.9998041131282007, + 1.4345968024678297, + 3.152708456359618 + ], + "result_full": [ + -1.1797070326091048, + 3.0491626549948174, + 5.255113951033274, + 6.222858340206228, + 3.0330887746483532, + -0.1852851367370416, + -4.351135706776737, + -0.9508699411899846, + -0.013940828237067326, + 1.936608302326945, + -2.562038930990207, + 1.8296643655997262, + 4.964562048480669, + 2.292125208665816, + -0.6908454273691205, + -3.6835021145119216, + -5.5628652488084125, + -5.381304213775536, + -6.067887771687595, + -8.07996643394984, + -7.08953178194481, + -10.528470437193244, + -12.554983229986522, + -13.893653047127422, + -14.392343973410926, + -11.792708351730493, + -14.670400116328647, + -14.606390104782424, + -13.798818710027865, + -18.110822790074945, + -17.17319300712088, + -20.327874953772945, + -23.31997659695665, + -23.40624531193906, + -21.213498734189034, + -19.19074781668703, + -17.7360646923138, + -18.26024981538604, + -12.405078452347947, + -16.394233974237434, + -20.870760171745225, + -22.84719683526064, + -22.71838868089418, + -23.277434240184192, + -28.59138991399194, + -24.305197395967138, + -23.030866628178117, + -19.84020035538177, + -17.539113593281897, + -20.69976356425527, + -15.628495665912101, + -11.034341067008846, + -7.873715061102328, + -0.6813686896583899, + 2.298616633098259, + 7.623253306630469, + 4.293877950846721, + -1.1313068425646469, + -2.515212580337982, + -4.256344451433963, + -6.141346277768917, + -4.03703972921134, + -1.9751209916296748, + -0.2411044036181078, + 3.7865468988791626, + -0.4633514984678804, + -4.834133477008061, + -7.699995878935996, + -5.840282601335777, + -2.0632951592104236, + -3.7864088162752196, + -3.217278962065742, + 0.31510946161231823, + 0.3382931581490354, + 4.6577407261346995, + 0.7798465301709483, + 2.9998041131282007, + 1.4345968024678297, + 3.152708456359618, + 2.9116988373103316, + 2.4224558840574915, + 5.310037658607117, + 0.8529368171173142, + -1.7545417989138952, + -5.872268498346905, + -9.536134152837285, + -10.395583181236091, + -14.26328305444989, + -10.319784550973093, + -7.405713997068804, + -3.1059820277852817, + -1.4522025473057558, + -0.402066181980288, + 1.7336491354499195, + -1.34695734082384 + ] +} \ No newline at end of file diff --git a/data/numpy/XCorrNumpyDataShortY.json b/data/numpy/XCorrNumpyDataShortY.json new file mode 100644 index 000000000..b6131ae28 --- /dev/null +++ b/data/numpy/XCorrNumpyDataShortY.json @@ -0,0 +1,271 @@ +{ + "x": [ + -1.254598811526375, + 4.507143064099161, + 2.3199394181140507, + 0.986584841970366, + -3.439813595575635, + -3.4400547966379733, + -4.419163878318005, + 3.6617614577493516, + 1.0111501174320878, + 2.0807257779604544, + -4.7941550570419755, + 4.699098521619943, + 3.324426408004218, + -2.8766088932172384, + -3.181750327928994, + -3.165954901465662, + -1.9575775704046228, + 0.24756431632237863, + -0.6805498135788426, + -2.0877085980195806, + 1.118528947223795, + -3.6050613934795814, + -2.0785535146478185, + -1.336381567063083, + -0.43930015782964027, + 2.8517596139301364, + -3.003262178416403, + 0.14234438413611628, + 0.9241456886204249, + -4.5354958728000225, + 1.0754485190143832, + -3.294758763127085, + -4.349484070147205, + 4.488855372533333, + 4.656320330745594, + 3.0839734811646107, + -1.9538623082662934, + -4.023278859936161, + 1.8423302651215687, + -0.5984750626039865, + -3.7796176515522117, + -0.048230898887298146, + -4.656114788847816, + 4.093204020787821, + -2.412200183999831, + 1.6252228435398202, + -1.8828892391058902, + 0.20068021177810813, + 0.4671027934327965, + -3.1514554447447294, + 4.695846277645586, + 2.7513282336111455, + 4.394989415641891, + 3.9482735042764876, + 0.9789997881108512, + 4.218742350231167, + -4.115074979480805, + -3.040171375808548, + -4.54772711089462, + -1.7466966923673564, + -1.1132271031051797, + -2.286509682261041, + 3.2873750915192943, + -1.432466733064107 + ], + "y": [ + 0.9403061933191574, + 1.0, + 0.9403061933191574 + ], + "result_valid": [ + 5.508889434467941, + 7.485725792601513, + -0.066439782834959, + -5.746826589072059, + -10.82989998849591, + -4.210691731774865, + 0.457185011471409, + 6.4108464302262265, + -1.4764571961116881, + 1.5809557214691905, + 3.317103570431141, + 5.038124692851211, + -2.7424496914840053, + -8.863610487595434, + -7.99849675376669, + -4.701758312139172, + -2.233083201602141, + -2.410848878259554, + -1.6758741061086728, + -4.234417932924314, + -4.507778439966936, + -6.725022934379439, + -3.7039349691650227, + 0.9856192048593386, + -0.3853030717292856, + -0.18788764559305482, + -1.8126617278510286, + -3.206761864359523, + -2.6552650552358994, + -6.287388410413921, + -6.373354669004378, + -3.226667632901103, + 4.777375408708889, + 11.777098202823785, + 5.62511149688694, + -2.8370969733431792, + -4.128153130858685, + -2.503533572240041, + -2.420118390303804, + -4.387719272398887, + -7.980402357778775, + -0.8526015105152269, + -2.5531763245113535, + 2.9648720125701, + -2.4134763418660965, + -0.16598128779007007, + -1.1305925515056185, + -2.307529433218497, + 1.703297542581664, + 4.3196041829011556, + 11.299497338456217, + 10.694666422549425, + 9.001468835395478, + 8.65849537704015, + 1.2698714250714382, + -3.006857392965021, + -11.185857832886702, + -9.048848802002134, + -7.06972699988755, + -4.905676036102568, + -0.24214486356286136, + -0.20960146461882134 + ], + "result_same": [ + 2.9834957258215495, + 5.508889434467941, + 7.485725792601513, + -0.066439782834959, + -5.746826589072059, + -10.82989998849591, + -4.210691731774865, + 0.457185011471409, + 6.4108464302262265, + -1.4764571961116881, + 1.5809557214691905, + 3.317103570431141, + 5.038124692851211, + -2.7424496914840053, + -8.863610487595434, + -7.99849675376669, + -4.701758312139172, + -2.233083201602141, + -2.410848878259554, + -1.6758741061086728, + -4.234417932924314, + -4.507778439966936, + -6.725022934379439, + -3.7039349691650227, + 0.9856192048593386, + -0.3853030717292856, + -0.18788764559305482, + -1.8126617278510286, + -3.206761864359523, + -2.6552650552358994, + -6.287388410413921, + -6.373354669004378, + -3.226667632901103, + 4.777375408708889, + 11.777098202823785, + 5.62511149688694, + -2.8370969733431792, + -4.128153130858685, + -2.503533572240041, + -2.420118390303804, + -4.387719272398887, + -7.980402357778775, + -0.8526015105152269, + -2.5531763245113535, + 2.9648720125701, + -2.4134763418660965, + -0.16598128779007007, + -1.1305925515056185, + -2.307529433218497, + 1.703297542581664, + 4.3196041829011556, + 11.299497338456217, + 10.694666422549425, + 9.001468835395478, + 8.65849537704015, + 1.2698714250714382, + -3.006857392965021, + -11.185857832886702, + -9.048848802002134, + -7.06972699988755, + -4.905676036102568, + -0.24214486356286136, + -0.20960146461882134, + 1.6586724252546174 + ], + "result_full": [ + -1.1797070326091048, + 2.9834957258215495, + 5.508889434467941, + 7.485725792601513, + -0.066439782834959, + -5.746826589072059, + -10.82989998849591, + -4.210691731774865, + 0.457185011471409, + 6.4108464302262265, + -1.4764571961116881, + 1.5809557214691905, + 3.317103570431141, + 5.038124692851211, + -2.7424496914840053, + -8.863610487595434, + -7.99849675376669, + -4.701758312139172, + -2.233083201602141, + -2.410848878259554, + -1.6758741061086728, + -4.234417932924314, + -4.507778439966936, + -6.725022934379439, + -3.7039349691650227, + 0.9856192048593386, + -0.3853030717292856, + -0.18788764559305482, + -1.8126617278510286, + -3.206761864359523, + -2.6552650552358994, + -6.287388410413921, + -6.373354669004378, + -3.226667632901103, + 4.777375408708889, + 11.777098202823785, + 5.62511149688694, + -2.8370969733431792, + -4.128153130858685, + -2.503533572240041, + -2.420118390303804, + -4.387719272398887, + -7.980402357778775, + -0.8526015105152269, + -2.5531763245113535, + 2.9648720125701, + -2.4134763418660965, + -0.16598128779007007, + -1.1305925515056185, + -2.307529433218497, + 1.703297542581664, + 4.3196041829011556, + 11.299497338456217, + 10.694666422549425, + 9.001468835395478, + 8.65849537704015, + 1.2698714250714382, + -3.006857392965021, + -11.185857832886702, + -9.048848802002134, + -7.06972699988755, + -4.905676036102568, + -0.24214486356286136, + -0.20960146461882134, + 1.6586724252546174, + -1.34695734082384 + ] +} \ No newline at end of file diff --git a/src/Numerics.Tests/StatisticsTests/CorrelationTests.cs b/src/Numerics.Tests/StatisticsTests/CorrelationTests.cs index ca1e34e2c..6fd209ef2 100644 --- a/src/Numerics.Tests/StatisticsTests/CorrelationTests.cs +++ b/src/Numerics.Tests/StatisticsTests/CorrelationTests.cs @@ -34,6 +34,8 @@ using MathNet.Numerics.Statistics; using MathNet.Numerics.TestData; using System.Globalization; +using System.Runtime.Serialization; +using System.Runtime.Serialization.Json; namespace MathNet.Numerics.UnitTests.StatisticsTests { @@ -86,6 +88,64 @@ public void AutoCorrelationTest(string fName, double tol) } } + [TestCase("numpy.XCorrNumpyData.json", CorrelationMode.Valid, "result_valid")] + [TestCase("numpy.XCorrNumpyData.json", CorrelationMode.Same, "result_same")] + [TestCase("numpy.XCorrNumpyData.json", CorrelationMode.Full, "result_full")] + [TestCase("numpy.XCorrNumpyDataShortY.json", CorrelationMode.Valid, "result_valid")] + [TestCase("numpy.XCorrNumpyDataShortY.json", CorrelationMode.Same, "result_same")] + [TestCase("numpy.XCorrNumpyDataShortY.json", CorrelationMode.Full, "result_full")] + public void CrossCorrelationTest(string filename, CorrelationMode mode, string resultMemberName) + { + // read the test data + var serializer = new DataContractJsonSerializer(typeof(XCorrTestdata)); + XCorrTestdata testdata; + using (var stream = Data.ReadStream(filename)) + { + testdata = serializer.ReadObject(stream) as XCorrTestdata; + } + Assume.That(testdata, Is.Not.Null); + + // select the test data for this test case + var x = testdata.x; + var y = testdata.y; + var expectedResult = typeof(XCorrTestdata).GetField(resultMemberName).GetValue(testdata) as double[]; + Assume.That(expectedResult, Is.Not.Null); + + var actualResult = Correlation.CrossCorrelation(x, y, mode); + Assert.That(actualResult, Is.EqualTo(expectedResult).Within(1e-10)); + } + + [TestCase(CorrelationMode.Valid)] + [TestCase(CorrelationMode.Same)] + [TestCase(CorrelationMode.Full)] + public void CrossCorrelationInvalid(CorrelationMode mode) + { + // The first sequence must be equally long or longer than the second sequence + Assert.Throws(() => + { + Correlation.CrossCorrelation(new double[5], new double[6], mode); + }); + Assert.DoesNotThrow(() => + { + Correlation.CrossCorrelation(new double[5], new double[5], mode); + }); + } + + [DataContract] + public class XCorrTestdata + { + [DataMember] + public double[] x; + [DataMember] + public double[] y; + [DataMember] + public double[] result_valid; + [DataMember] + public double[] result_same; + [DataMember] + public double[] result_full; + } + [TestCase()] public void AutoCorrelationTest2() { diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index 369acd391..799ab8203 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2019 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -37,6 +37,18 @@ namespace MathNet.Numerics.Statistics { + + /// + /// Correlation modes for calculating the cross correlation + /// + public enum CorrelationMode + { + Valid, + Same, + Full + } + + /// /// A class with correlation measures between two datasets. /// @@ -160,6 +172,95 @@ private static double[] AutoCorrelationFft(double[] x, int kLow, int kHigh) return result; } + /// + /// Cross correlation of two real-valued one-dimensional inputs. + /// + /// First input + /// Second input (cannot be longer than x) + /// Padding of the signal + /// An array with the cross-correlation of the two input signals + public static double[] CrossCorrelation(IList x, IList y, CorrelationMode mode = CorrelationMode.Valid) + { + + if (y.Count > x.Count) + { + throw new ArgumentException("y must not be longer than x"); + } + + var iStart = 0; + var iEnd = x.Count - y.Count + 1; + + if (mode == CorrelationMode.Same) + { + iStart = -(y.Count / 2); + iEnd = iStart + x.Count; + } + else if (mode == CorrelationMode.Full) + { + iStart = -(y.Count - 1); + iEnd = x.Count; + } + + var n = iEnd - iStart; + var result = new double[n]; + + // decide whether to calculate the cross correlation in time-domain or frequency-domain. + // Since FFT is calculated in n * log2(n), we choose log2(x) as the threshold for the + // length of y. Some tweaking of this threshold might help with performance. + var timedomain = y.Count < Math.Log(x.Count) / Math.Log(2); + + if (timedomain) + { + // this is a naive implementation in time-domain. This is actually only fast + // when y is much shorter than x. Whenever x and y are approximately the same + // length an FFT based approach would probably be faster. + Threading.CommonParallel.For(0, n, (a, b) => + { + for (int i = a; i < b; i++) + { + for (int k = 0; k < y.Count; k++) + { + var xi = i + k + iStart; + if (xi < 0 || xi >= x.Count) continue; + result[i] += x[xi] * y[k]; + } + } + }); + } + else + { + // use an FFT based implementation when x and y are within the same order + // of magnitude + var padLeft = -iStart; + var padRight = iEnd - x.Count + y.Count - 1; + var nFFT = Euclid.CeilingToPowerOfTwo(padLeft + x.Count + padRight); + var xFFT = new Complex[nFFT]; + var yFFT = new Complex[nFFT]; + for (int i = 0; i < x.Count; i++) + { + xFFT[i + padLeft] = x[i]; + } + for (int i = 0; i < y.Count; i++) + { + yFFT[i] = y[i]; + } + Fourier.Forward(xFFT, FourierOptions.Matlab); + Fourier.Forward(yFFT, FourierOptions.Matlab); + var zFFT = new Complex[nFFT]; + for (int i = 0; i < nFFT; i++) + { + zFFT[i] = xFFT[i] * yFFT[i].Conjugate(); + } + Fourier.Inverse(zFFT, FourierOptions.Matlab); + for (int i = 0; i < n; i++) + { + result[i] = zFFT[i].Real; + } + } + + return result; + } + /// /// Computes the Pearson Product-Moment Correlation coefficient. /// diff --git a/src/TestData/TestData.csproj b/src/TestData/TestData.csproj index ebbec8f2b..3aaeaf0c9 100644 --- a/src/TestData/TestData.csproj +++ b/src/TestData/TestData.csproj @@ -1,4 +1,4 @@ - + Library net40;netstandard1.1 @@ -28,6 +28,12 @@ Data\numpy\CorrNumpyData_rnd.csv + + Data\numpy\XCorrNumpyData.json + + + Data\numpy\XCorrNumpyDataShortY.json + Data\Codeplex-5667.csv