From a1e359231e867425812f4a025d442065e71b6d1e Mon Sep 17 00:00:00 2001
From: sunwen <sunwen@realcan.cn>
Date: Thu, 18 May 2023 11:42:48 +0800
Subject: [PATCH] Add xcorr function.

---
 src/Function1D.cpp | 34 ++++++++++++++++++++++++++++++++++
 src/Function1D.h   |  2 ++
 2 files changed, 36 insertions(+)

diff --git a/src/Function1D.cpp b/src/Function1D.cpp
index bee37ef..16e2768 100644
--- a/src/Function1D.cpp
+++ b/src/Function1D.cpp
@@ -735,6 +735,40 @@ Matrix Aurora::intersect(const Matrix& aMatrix1, const Matrix& aMatrix2, Matrix&
     return result;
 }
 
+Matrix Aurora::xcorr(const Matrix& aMatrix1, const Matrix& aMatrix2)
+{
+    //not support for complex
+    if (aMatrix1.isNull() || aMatrix2.isNull() || aMatrix1.getDataSize() != aMatrix2.getDataSize())
+    {
+        return Matrix();
+    }
+    size_t matrixSize = aMatrix1.getDataSize();
+    size_t resultSize = 2 * matrixSize -1;
+    double* resultData = Aurora::malloc(resultSize);
+    for(int i=0;i<matrixSize;++i)
+    {
+        double data = 0;
+        for(int j=0;j<i+1;++j)
+        {
+            data+= aMatrix1[j] * aMatrix2[matrixSize-i-1+j];
+        }
+        resultData[i] = data;
+    }
+
+    for(int i=0;i<matrixSize-1;++i)
+    {
+        double result = 0;
+        for(int j=0;j<i+1;++j)
+        {
+            result+= aMatrix1[matrixSize-i-1+j]*aMatrix2[j];
+
+        }
+        resultData[resultSize - 1 - i] = result;
+    }
+
+    return Matrix::New(resultData,resultSize);
+}
+
 Matrix Aurora::reshape(const Matrix& aMatrix, int aRows, int aColumns, int aSlices)
 {
     if(aMatrix.isNull() || (aMatrix.getDataSize() != aRows * aColumns * aSlices))
diff --git a/src/Function1D.h b/src/Function1D.h
index 27fb08b..75cce0b 100644
--- a/src/Function1D.h
+++ b/src/Function1D.h
@@ -80,6 +80,8 @@ namespace Aurora {
     Matrix intersect(const Matrix& aMatrix1, const Matrix& aMatrix2);
 
     Matrix reshape(const Matrix& aMatrix, int aRows, int aColumns, int aSlices);
+
+    Matrix xcorr(const Matrix& aMatrix1, const Matrix& aMatrix2);
     /**
      * 并集
      * @param aIa, [C,ia,~] = intersect(A,B)用法中ia的返回值