From 24b7d0b3b5eb33a9bafdc82a04fc3ad885c5192d Mon Sep 17 00:00:00 2001
From: Krad <kradchen123@163.com>
Date: Thu, 27 Apr 2023 14:40:11 +0800
Subject: [PATCH] Create match filter test in main.

---
 src/main.cxx | 120 ++++++++++++++++++++++++---------------------------
 1 file changed, 56 insertions(+), 64 deletions(-)

diff --git a/src/main.cxx b/src/main.cxx
index a6c32c3..5fb3b51 100644
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -10,78 +10,70 @@
 #include "Matrix.h"
 #include "Function.h"
 #include "Function1D.h"
+#include "Function2D.h"
+#include "Function3D.h"
 
 #include "spdlog/spdlog.h"
 #include "spdlog/sinks/stdout_color_sinks.h"
 #include "spdlog/sinks/basic_file_sink.h"
+#include "MatlabReader.h"
 
 int main() {
-    {
-        double *dataA = Aurora::malloc(8,true);
-        double *dataB = Aurora::malloc(8);
-        double *dataC = Aurora::malloc(8);
-        for (int i = 0; i < 16; ++i) {
-            dataA[i] = (double) (i + 2);
+    MatlabReader m;
+    Input i;
+    MatchedFilter o;
+
+    m.read(&i,&o,"/mnt/d/TestData/testCreateMatchedFilter.mat");
+    bool measuredCEused = true,findDefects = i.mParams->mFindDefects,removeOutliers = i.mRemoveOutliersFromCEMeasured;
+
+    Aurora::Matrix mFTime = Aurora::Matrix::fromRawData(i.mCe,4000,2304);
+    if (removeOutliers){
+        auto normSTD = Aurora::std(Aurora::abs(mFTime));
+        Aurora::nantoval(normSTD,0.0);
+        auto sortSTD = Aurora::sort(normSTD);
+        int t = (int)std::round(0.4*mFTime.getDimSize(1))-1;
+        t = t<=0?1.0:t;
+        auto absFTime = abs(mFTime);
+        auto maxAbsFTime = max(absFTime);
+        auto maxFlag = maxAbsFTime < (0.1 * max(maxAbsFTime));
+        auto lessFlag = normSTD < sortSTD(0,t).toMatrix();
+        long maxCol,maxRow;
+        max(normSTD,Aurora::Column,maxRow,maxCol);
+        for (int j = 0; j < mFTime.getDimSize(1); ++j) {
+            if((bool)(lessFlag.getData()[j])||(bool)(maxFlag.getData()[j])){
+                mFTime(Aurora::$,j) =  mFTime(Aurora::$,maxCol);
+            }
         }
-        for (int i = 0; i < 8; ++i) {
-            dataB[i] = (double) (i + 2);
-            dataC[i / 2] = (double) (i + 9);
-        }
-        Aurora::Matrix A = Aurora::Matrix::New(dataA, 2, 2, 2,Aurora::ValueType::Complex);
-        printf("A:\r\n");
-        A.printf();
-        Aurora::Matrix B = Aurora::Matrix::New(dataB, 2, 2, 2);
-        printf("B:\r\n");
-        B.printf();
-        printf("A slice 1 = B slice 0\r\n");
-        A(Aurora::$, Aurora::$, 1) = B(Aurora::$, Aurora::$, 0);
-        printf("B:\r\n");
-        B.printf();
-        printf("New A:\r\n");
-        A.printf();
-        printf("A Row slice 1 = B Row slice 0\r\n");
-        A(1, Aurora::$, Aurora::$) = B(0, Aurora::$, Aurora::$);
-        printf("B:\r\n");
-        B.printf();
-        printf("New A:\r\n");
-        A.printf();
-
-        printf("A Col slice 1 = B Col slice 0\r\n");
-        A(Aurora::$, 1, Aurora::$) = B(Aurora::$, 0, Aurora::$);
-        printf("B:\r\n");
-        B.printf();
-        printf("New A:\r\n");
-        A.printf();
-
-        printf("New A col slice 1 toMatrix:\r\n");
-        auto Ds = A(Aurora::$, 1, Aurora::$);
-        auto D = Ds.toMatrix();
-        printf("D:\r\n");
-        D.printf();
-
-
-        Aurora::Matrix C = Aurora::Matrix::New(dataC, 2, 2);
-        printf("C:\r\n");
-        C.printf();
-        A(1, 1, 1) = 1024.0;
-        printf("New A:\r\n");
-        A.printf();
-    }{
-        double * dataA = Aurora::malloc(9);
-        double * dataB = Aurora::malloc(9);
-        for (int i = 0; i < 9; ++i) {
-            dataA[i]=(double)(i-3);
-            dataB[i]=(double)(i+2);
-        }
-        Aurora::Matrix A = Aurora::Matrix::New(dataA,3,3);
-        printf("A:\r\n");
-        A.printf();
-        Aurora::Matrix B = Aurora::Matrix::New(dataB,3,3);
-        printf("B:\r\n");
-        B.printf();
-        printf("sign(A):\r\n");
-        sign(A*B).printf();
+    }
+    auto matchedFilter = fft(mFTime);
+    auto absMatchedFilter = abs(matchedFilter);
+//    auto highNoiseScore = mean(absMatchedFilter) * Aurora::std(absMatchedFilter);
+    auto highNoiseScore = mean( abs(matchedFilter)) * Aurora::std(absMatchedFilter);
+    printf("highNoiseScore[998]: %f4, highNoiseScore[2053]:%f4\r\n",highNoiseScore.getData()[998],highNoiseScore.getData()[2053]);
+    auto medianNoise = median(highNoiseScore);
+    printf("median: %f4\r\n",medianNoise.getScalar());
+    long minCol,row;
+    min(abs(highNoiseScore - median(highNoiseScore)), Aurora::Column, row, minCol);
+    //
+    minCol = 998;
+    auto sumDiff = Aurora::zeros(1, matchedFilter.getDimSize(1));
+    auto maxMatchFilter = matchedFilter(Aurora::$, minCol).toMatrix();
+    for (int k = 0; k < matchedFilter.getDimSize(1); ++k) {
+        sumDiff.getData()[k] = sum(abs(matchedFilter(Aurora::$, k).toMatrix()- maxMatchFilter)).getData()[0];
     }
 
+    auto indexe = sumDiff > (mean(sumDiff)+2.596 * Aurora::std(sumDiff));
+    for (int l = 0; l < indexe.getDataSize(); ++l) {
+        if ((bool)indexe.getData()[l]){
+            matchedFilter(Aurora::$,l) = matchedFilter(Aurora::$,minCol);
+        }
+    }
+    if (measuredCEused){
+        auto mFTime2 = ifft(-matchedFilter);
+        max(Aurora::abs(mFTime2)).printfShape();
+        mFTime2 = mFTime2/repmat(max(Aurora::abs(mFTime2)),mFTime2.getDimSize(0),1);
+        mFTime2 = mFTime2/repmat(sum(Aurora::abs(mFTime2)),mFTime2.getDimSize(0),1);
+        matchedFilter = fft(mFTime2);
+    }
     return 0;
 }