mix code vision 3. Using the previous algorithm to detect peaks as Nikoleta and Shiyao. Adding overlapping windows

Dependencies:   mpu9250_i2c biquadFilter peakdetection Eigen

Revision:
2:d4c480d17944
Parent:
1:92f42e198925
Child:
4:15d6c7123b09
diff -r 92f42e198925 -r d4c480d17944 main.cpp
--- a/main.cpp	Thu Nov 14 01:28:23 2019 +0000
+++ b/main.cpp	Mon Nov 25 14:28:46 2019 +0000
@@ -14,84 +14,25 @@
 #include "platform/mbed_thread.h"
 #include "stats_report.h"
 #include "MPU9250.h"
-#include <Eigen/Dense.h>
+//#include <Eigen/Dense.h>
 #include <iostream>
+#include <vector>  
+#include <complex>
+#include "BiQuad.h"
+#include "pca.h"
+#include "peak.h"
 
 using namespace std;
 using namespace Eigen;
 
 DigitalOut led1(LED1);
 const int addr7bit = 0x68; // 7bit I2C address,AD0 is 0
-
-#define SLEEP_TIME                  5000 // (msec)
+//the parameter of biquad filter, 40Hz sampling frequence,10Hz cut-off freq, Q:0.719
+BiQuad mybq(0.3403575989782886,0.6807151979565772,0.3403575989782886, -1.511491371967327e-16,0.36143039591315457);
+BiQuadChain bqc;
 
 
-/*
-  * Normalize the Matrix X
-  */
- MatrixXd featurnormail(MatrixXd &X)
-{
-    //I don't know why to use the transpose
-    //compute the mean of every dimension
-    MatrixXd X1 = X.transpose();
-    MatrixXd meanval = X1.colwise().mean();
-    
-    //normalization
-    RowVectorXd meanvecRow = meanval;
-    X1.rowwise() -= meanvecRow;
-    
-    return X1.transpose();
-}
-
- /*
-  * Compute the Covariane Matrix of X, put to C
-  * C = 1/m * X * X.transpose 
-  */
-void ComComputeCov(MatrixXd &X, MatrixXd &C)
-{
-   
-    C = X*X.adjoint();//same as XT*X a
-    //translate to array
-    C = C.array() / X.cols();
-}
- 
- 
-/*
- * Compute the eigenvalue and eigenvector of C
- * val = (first eigenvalue) --smallest --not important
- *              .
- *              .
- *              .
- *       (last eigenvalue)  --largest -- important
- *
- * vec = (first eigenvector, ... , last eigenvector)
- *           not important          important
- */
-void ComputEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
-{
-    //SelfAdjointEigenSolver will sort the values automatically
-    SelfAdjointEigenSolver<MatrixXd> eig(C);
-    vec = eig.eigenvectors();
-    val = eig.eigenvalues();
-}
- 
-/* Compute the dimension need to include enough information of raw data.
- * form large index to small index, since the val is sorted from small to large.
- * in some cases, just decide the number of dimension, instead of compute it.
- */
-int ComputDim(MatrixXd &val)
-{
-    int dim;
-    double sum = 0;
-    for (int i = val.rows() - 1; i >= 0;--i)
-    {
-        sum += val(i, 0);
-        dim = i;
-        if (sum / val.sum()>=0.8)//80% of the information
-            break;
-    }
-    return val.rows() - dim;
-}
+#define SLEEP_TIME                  20 // (msec)
 
 
 // main() runs in its own thread in the OS
@@ -105,24 +46,40 @@
     float AccRead[3];
     float GyroRead[3];
     float TempRead[1];
-    
+    float res_smooth;
+    //vector<float> res_list;
+    float number=0;
     
-    MatrixXd acc_raw(3,0);
+    static MatrixXd acc_raw(3,0);
     Vector3d acc_new;
     MatrixXd C;
     MatrixXd vec, val;
     int dim = 1;    //dimension of PCA
+    //use the class defined in pca.h and peak.h
+    PCA pca;
+    PEAK peak;
+    
+    bqc.add(&mybq);
+    
+    static vector<float> res_list;
     
     while (true) {
         
+        //when i add to 600, there accure faugment error
+        //printf("the %f loop", number);
+        number = number +1;
+        //vector<float> res_list;
         //Blink LED and wait 1 seconds
         led1 = !led1;
         thread_sleep_for(SLEEP_TIME);
         //read and convert date
         mpu->ReadConvertAll(AccRead,GyroRead,TempRead);
+        AccRead[0]= AccRead[0]/1000;
+        AccRead[1]= AccRead[1]/1000;
+        AccRead[2]= AccRead[2]/1000;
         printf("acc value is (%f,%f,%f).\n\r",AccRead[0],AccRead[1],AccRead[2]);
-        printf("gyro value is (%f,%f,%f).\n\r",GyroRead[0],GyroRead[1],GyroRead[2]);
-        printf("temp value is %f.\n\r",TempRead[0]);
+        //printf("gyro value is (%f,%f,%f).\n\r",GyroRead[0],GyroRead[1],GyroRead[2]);
+        //printf("temp value is %f.\n\r",TempRead[0]);
         
         //append new data to matrix acc_raw
         //adding the columns
@@ -130,17 +87,30 @@
         acc_raw.conservativeResize(acc_raw.rows(), acc_raw.cols()+1);
         acc_raw.col(acc_raw.cols()-1) = acc_new;
         
-        cout << "acc_raw:" << acc_raw << endl;
+        //////cout << "acc_raw:" << acc_raw << endl;
         
         //run PCA
-        MatrixXd X1=featurnormail(acc_raw);
-        ComComputeCov(X1, C);
-        ComputEig(C, vec, val);
+        MatrixXd X1=pca.featurnormail(acc_raw);
+        pca.ComComputeCov(X1, C);
+        pca.ComputEig(C, vec, val);
         //select dim num of eigenvector from right to left. right is important
         //compute the result array
         MatrixXd res = vec.rightCols(dim).transpose()*X1;
         
         //show the result after PCA
-        cout << "result" << res << endl;
+        //////cout << "result" << res << endl;
+        res_list.clear();
+        for(int i = 0; i < res.cols(); i++)
+        {
+            res_smooth = bqc.step(res(i));
+            res_list.push_back(res_smooth);
+            //printf("result after filter in for loop %d: %f\n\r",i,res_smooth);
+            //std::cout << "\t" << bqc.step(  ) << std::endl;
+        }
+        int len = res_list.size();
+        printf("len of res:%d\n\r", len);
+        peak.findPeaks(res_list,len,0.1);
+        
+        
     }
 }