Add ifft_symmetric and its unit test.

2023-05-05 13:22:15 +08:00
parent c967788302
commit 5353aa53d5
4 changed files with 70 additions and 30 deletions
--- 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));
--- 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
--- 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;
    }
--- 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,