From 5353aa53d566bea885c29c9ec09c6e4fdbc8e5ad Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Fri, 5 May 2023 13:22:15 +0800
Subject: [PATCH] Add ifft_symmetric and its unit test.

---
 src/Function2D.cpp       | 24 ++++++++++++++-
 src/Function2D.h         | 63 ++++++++++++++++++++++------------------
 src/Matrix.cpp           |  2 +-
 test/Function2D_Test.cpp | 11 +++++++
 4 files changed, 70 insertions(+), 30 deletions(-)

diff --git a/src/Function2D.cpp b/src/Function2D.cpp
index 059a013..e3b8076 100644
--- a/src/Function2D.cpp
+++ b/src/Function2D.cpp
@@ -670,7 +670,7 @@ Matrix Aurora::ifft(const Matrix &aMatrix) {
     //提交 修改配置后的Descriptor(实际上会进行FFT的计算初始化)
     status = DftiCommitDescriptor(my_desc_handle);
     if (status != DFTI_NO_ERROR) goto error;
-
+    
     //执行计算
     status = DftiComputeBackward(my_desc_handle, aMatrix.getData(), output);
     if (status != DFTI_NO_ERROR) goto error;
@@ -685,6 +685,28 @@ Matrix Aurora::ifft(const Matrix &aMatrix) {
     return Matrix();
 }
 
+Matrix Aurora::ifft_symmetric(const Matrix &aMatrix,long length)
+{
+    if(!aMatrix.isVector()){
+        std::cerr<<"ifft_symmetric only support vector!"<<std::endl;
+        return Matrix();
+    }
+    int matrixLength = aMatrix.getDataSize();
+    int resultHalfLength = (int)std::ceil(((double)length*0.5));
+    int copyLength = resultHalfLength<matrixLength?resultHalfLength:matrixLength;
+    double* calcData = malloc(length,true);
+    double zero = 0.0;
+    //所有数据统一置0
+    cblas_dcopy(length*2,&zero,0,calcData,1);
+    //copy前半段数据
+    cblas_zcopy(copyLength,aMatrix.getData(),1,calcData,1);
+    //copy后半段数据,跳过index 0的值,并设置虚部共轭
+    vdAddI(copyLength-1,&zero,0,(aMatrix.getData()+2),2,(calcData+(length-1)*2),-2);
+    vdSubI(copyLength-1,&zero,0,(aMatrix.getData()+2+1),2,(calcData+(length-1)*2+1),-2);
+
+    return real(ifft(Matrix::New(calcData,length,1,1,Complex)));
+}
+
 Matrix Aurora::hilbert(const Matrix &aMatrix) {
     auto x = fft(aMatrix);
     auto h = malloc(aMatrix.getDimSize(0));
diff --git a/src/Function2D.h b/src/Function2D.h
index f87f60a..58bd7e3 100644
--- a/src/Function2D.h
+++ b/src/Function2D.h
@@ -4,19 +4,20 @@
 #include "Matrix.h"
 #include "Function1D.h"
 
-
-namespace Aurora {
-    enum FunctionDirection{
+namespace Aurora
+{
+    enum FunctionDirection
+    {
         Column,
         Row,
         All
     };
-    double immse(const Matrix& aImageA, const Matrix& aImageB);
-    Matrix inv(const Matrix& aMatrix);
-    Matrix inv(Matrix&& aMatrix);
-    Matrix interp2(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod);
-    Matrix interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod);
-    Matrix std(const Matrix& aMatrix);
+    double immse(const Matrix &aImageA, const Matrix &aImageB);
+    Matrix inv(const Matrix &aMatrix);
+    Matrix inv(Matrix &&aMatrix);
+    Matrix interp2(const Matrix &aX, const Matrix &aY, const Matrix &aV, const Matrix &aX1, const Matrix &aY1, InterpnMethod aMethod);
+    Matrix interpn(const Matrix &aX, const Matrix &aY, const Matrix &aV, const Matrix &aX1, const Matrix &aY1, InterpnMethod aMethod);
+    Matrix std(const Matrix &aMatrix);
 
     /**
      * 求矩阵最小值，可按行、列、单元, 目前不支持三维，不支持复数
@@ -24,9 +25,9 @@ namespace Aurora {
      * @param direction 方向，Column, Row, All
      * @return
      */
-    Matrix min(const Matrix& aMatrix,FunctionDirection direction = Column);
+    Matrix min(const Matrix &aMatrix, FunctionDirection direction = Column);
 
-    Matrix min(const Matrix& aMatrix,FunctionDirection direction, long& rowIdx, long& colIdx);
+    Matrix min(const Matrix &aMatrix, FunctionDirection direction, long &rowIdx, long &colIdx);
 
     /**
      * 求矩阵最小值，可按行、列、单元, 目前不支持三维，不支持复数
@@ -34,9 +35,9 @@ namespace Aurora {
      * @param direction 方向，Column, Row, All
      * @return 最大值矩阵
      */
-    Matrix max(const Matrix& aMatrix,FunctionDirection direction = Column);
+    Matrix max(const Matrix &aMatrix, FunctionDirection direction = Column);
 
-    Matrix max(const Matrix& aMatrix,FunctionDirection direction , long& rowIdx, long& colIdx);
+    Matrix max(const Matrix &aMatrix, FunctionDirection direction, long &rowIdx, long &colIdx);
 
     /**
      * 比较两个矩阵，求对应位置的最小值，不支持三维
@@ -45,7 +46,7 @@ namespace Aurora {
      * @param aOther 目标矩阵2
      * @return 最小值矩阵
      */
-    Matrix min(const Matrix& aMatrix,const Matrix& aOther);
+    Matrix min(const Matrix &aMatrix, const Matrix &aOther);
 
     /**
      * 求矩阵和，可按行、列、单元, 目前不支持三维
@@ -53,7 +54,7 @@ namespace Aurora {
      * @param direction 方向，Column, Row, All
      * @return 求和结果矩阵
      */
-    Matrix sum(const Matrix& aMatrix,FunctionDirection direction = Column);
+    Matrix sum(const Matrix &aMatrix, FunctionDirection direction = Column);
 
     /**
      * 求矩阵平均值，可按行、列、单元, 目前不支持三维，不支持复数
@@ -62,49 +63,49 @@ namespace Aurora {
      * @param aIncludeNan 是否包含nan
      * @return 平均值矩阵
      */
-    Matrix mean(const Matrix& aMatrix,FunctionDirection direction = Column, bool aIncludeNan = true);
+    Matrix mean(const Matrix &aMatrix, FunctionDirection direction = Column, bool aIncludeNan = true);
 
     /**
      * 矩阵排序 按列, 目前不支持三维，不支持复数
      * @param aMatrix 目标矩阵
      * @return 排序后矩阵
      */
-    Matrix sort(const Matrix& aMatrix);
+    Matrix sort(const Matrix &aMatrix);
 
     /**
      * 矩阵排序 按列, 目前不支持三维，不支持复数
      * @param aMatrix 目标矩阵
      * @return 排序后矩阵
      */
-    Matrix sort(Matrix&& aMatrix);
+    Matrix sort(Matrix &&aMatrix);
 
     /**
      * 矩阵排序 按行, 目前不支持三维，不支持复数
      * @param aMatrix 目标矩阵
      * @return 排序后矩阵
      */
-    Matrix sortrows(const Matrix& aMatrix);
+    Matrix sortrows(const Matrix &aMatrix);
 
     /**
      * 矩阵排序 按行, 目前不支持三维，不支持复数
      * @param aMatrix 目标矩阵
      * @return 排序后矩阵
      */
-    Matrix sortrows(Matrix&& aMatrix);
+    Matrix sortrows(Matrix &&aMatrix);
 
     /**
      * 对矩阵求中间值 按列, 目前不支持三维，不支持复数
      * @param aMatrix 目标矩阵
      * @return 中值矩阵
      */
-    Matrix median(const Matrix& aMatrix);
+    Matrix median(const Matrix &aMatrix);
 
     /**
      * FFT,支持到2维，输入可以是常数可以是复数，输出必是复数
      * @param aMatrix 目标矩阵
      * @return fft后的复数矩阵
      */
-    Matrix fft(const Matrix& aMatrix);
+    Matrix fft(const Matrix &aMatrix);
 
     /**
      * 逆fft，支持到2维，输入必须是复数，输出必是复数
@@ -112,24 +113,30 @@ namespace Aurora {
      * @param aMatrix
      * @return ifft后的复数矩阵
      */
-    Matrix ifft(const Matrix& aMatrix);
+    Matrix ifft(const Matrix &aMatrix);
+
+    /**
+     * Symmetric逆fft，支持到2维，输入必须是复数，输出必是实数
+     * @param aMatrix
+     * @return ifft后的实数矩阵
+     */
+    Matrix ifft_symmetric(const Matrix &aMatrix,long length);
 
     /**
      * hilbert，支持到2维，输入必须是复数，输出必是复数
      * @param aMatrix
      * @return
      */
-    Matrix hilbert(const Matrix& aMatrix);
+    Matrix hilbert(const Matrix &aMatrix);
 
     /**
      * prod，支持到2维
      * @param aMatrix
      * @return
      */
-    Matrix prod(const Matrix& aMatrix);
+    Matrix prod(const Matrix &aMatrix);
 
-    Matrix dot(const Matrix& aMatrix,const Matrix& aOther,FunctionDirection direction = Column);
+    Matrix dot(const Matrix &aMatrix, const Matrix &aOther, FunctionDirection direction = Column);
 };
 
-
-#endif //AURORA_FUNCTION2D_H
+#endif // AURORA_FUNCTION2D_H
diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index 1bcfb04..b806c5f 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -266,7 +266,7 @@ namespace Aurora {
 
     bool Matrix::isVector() const{
         if (getDimSize(2)>1) return false;
-        if (isScalar) return false;
+        if (isScalar()) return false;
         return getDimSize(0) == 1 ||
                 getDimSize(1) == 1;
     }
diff --git a/test/Function2D_Test.cpp b/test/Function2D_Test.cpp
index 9f62724..c4043fc 100644
--- a/test/Function2D_Test.cpp
+++ b/test/Function2D_Test.cpp
@@ -363,6 +363,17 @@ TEST_F(Function2D_Test, hilbert) {
     EXPECT_DOUBLE_EQ(fourDecimalRound(result[11].imag()),0.3249);
 }
 
+TEST_F(Function2D_Test, ifft_symmetric) {
+    double *input =  new double[18]{10,2,1,3,4,4,16,3,1,2,15,-2,1,-3,4,-4,1,-3};
+    auto ma = Aurora::Matrix::fromRawData(input,9,1,1,Aurora::Complex);
+    auto ret = Aurora::ifft_symmetric(ma,18);
+    auto result = ret.getData();
+    EXPECT_DOUBLE_EQ(fourDecimalRound(result[0]),5.3333);
+    EXPECT_DOUBLE_EQ(fourDecimalRound(result[1]),1.1188);
+    EXPECT_DOUBLE_EQ(fourDecimalRound(result[11]),2.8506);
+    EXPECT_DOUBLE_EQ(fourDecimalRound(result[17]),1.1188);
+}
+
 TEST_F(Function2D_Test, prod) {
     double *dataB = new double[20]{1.1, 2.6, 3.8, 6.2,
                                    4.3, 5.7, 6.9, 10.6,