compare against mag and avg

go.sh

@@ -71,7 +71,7 @@ function BuildMRI {
     # make -j4 
     # make install
     # g++ $path_src/main.cc -o $mri_bin 
-    OMPI_CXXFLAGS=$flags_O2 mpic++ -o $mri_bin $path_src/main.cc $path_src/jade.cc $path_src/generator-rf.cc
+    OMPI_CXXFLAGS=$flags_O2 mpic++ -o $mri_bin $path_src/main.cc $path_src/jade.cc $path_src/generator-rf.cc $path_src/lock-in.cc
     cp $mri_bin $path_bin
     cd $path_tmp
 }  # end of function BuildMRI
@@ -84,4 +84,4 @@ MPI_options=
 #MPI_options="--bind-to-core"
 MPI_size=1
 cd $path_bin
-time mpirun -np $MPI_size $MPI_options ./$mri_bin #mri.config
+time mpirun -np $MPI_size $MPI_options ./$mri_bin > out.txt 

src/generator-rf.cc

@@ -6,6 +6,39 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void GeneratorRF::EvaluateDemodulated(std::vector<double> &mag_sin,
+                                        std::vector<double> &mag_cos) {
+    mag_sin.resize(total_samples_);
+    mag_cos.resize(total_samples_);
+    for (long i = 0; i < total_samples_; ++i) {
+      mag_sin[i] = 0.0;
+      mag_cos[i] = 0.0;
+      double time = timesteps_[i];
+      for (unsigned int voxel = 0; voxel < voxel_num_; ++voxel) {
+        mag_sin[i] +=
+          amplitudes_[voxel]
+          * std::sin(phases_[voxel]
+                     )
+          * std::exp(-time/
+                     (T2s_decays_norm_[voxel]*T2s_scale_
+                      )
+                     )
+          ;
+        mag_cos[i] +=
+          amplitudes_[voxel]
+          * std::cos(phases_[voxel]
+                     )
+          * std::exp(-time/
+                     (T2s_decays_norm_[voxel]*T2s_scale_
+                      )
+                     )
+          ;
+       }
+    }
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   void GeneratorRF::EvaluateSignal(std::vector<double> &signal_rf) {
     signal_rf.resize(total_samples_);
     for (long i = 0; i < total_samples_; ++i) {
@@ -41,12 +74,16 @@ namespace mri {
     T2s_decays_norm_.clear();
     amplitudes_.resize(voxel_num_);
     T2s_decays_norm_.resize(voxel_num_);
+    std::cout << "#";
     for (double &amp:amplitudes_){
       amp = rand(om);
+      std::cout << amp << "\t";
     }
     for (double &decay:T2s_decays_norm_){
       decay = rand(om);
+      std::cout << decay << "\t";
     }
+    std::cout<<std::endl;
   };
 
   // ********************************************************************** //

src/generator-rf.h

@@ -5,14 +5,17 @@
 namespace mri {
   class GeneratorRF {
   public:
+    template<class T> inline T pow2(const T value) {return value*value;}
     const double PI_=3.14159265358979323846;
     void EvaluateSignal(std::vector<double> &signal_rf);
+    void EvaluateDemodulated(std::vector<double> &mag_sin,
+                             std::vector<double> &mag_cos);
     void RandomizeVoxels();
     void SetAmplitudes(std::vector<double> amplitudes);
     void SetT2sDecaysNorm(std::vector<double> T2s_decays_norm);
     void SetTimer(long total_periods, double period_ms, int samples_per_period);
     void SetVoxels(double voxel_num, double phase_init_, double phase_range_, double T2s_scale);
-    const std::vector<double>& GetTimesteps(){return timesteps_};
+    const std::vector<double>& GetTimesteps(){return timesteps_;};
   private:
     std::vector<double> amplitudes_;
     std::vector<double> phases_;

src/lock-in.cc

@@ -0,0 +1,105 @@
+#include "./lock-in.h"
+#include <cmath>
+
+namespace mri {
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void LockIn::Analyse(const std::vector<double> &signal,double period_ms,
+                       const std::vector<double> &timesteps,
+                       int samples_per_period,
+                       std::vector<double> &avg_sin,
+                       std::vector<double> &avg_cos) {
+    period_ms_=period_ms;
+    total_samples_ = signal.size();
+    // Stage 1: Multiply
+    GenerateSinCos(timesteps);
+    std::vector<double> mul_sin, mul_cos;
+    mul_sin.resize(total_samples_);
+    mul_cos.resize(total_samples_);
+    for (unsigned long i = 0; i < total_samples_; ++i) {
+      mul_sin[i] = signal[i]*lockin_sin_[i];
+      mul_cos[i] = signal[i]*lockin_cos_[i];
+    }
+    // Stage 2: Low-pas filter:
+    //y[n] = a0 x[n] + b1 y[n-1]
+    // a0 = 1-x
+    // b1 = x
+    // x = exp (-1/samples)
+    /// a0 = (1-x)^4
+    // b1 = 4x
+    // b2 = -6 x^2
+    // b3 = 4 x^3
+    // b4 = - x^4
+    double x = std::exp(-1.0/(samples_per_period*2));
+    std::vector<double> mag_sin, mag_cos;
+    mag_sin.resize(total_samples_);
+    mag_cos.resize(total_samples_);
+    avg_sin.resize(total_samples_);
+    avg_cos.resize(total_samples_);
+    double mean_sin = 0.0, mean_cos = 0.0;
+    for (int i = 0; i < samples_per_period; ++i){
+      mag_sin[i] = mul_sin[i];
+      mag_cos[i] = mul_cos[i];
+      avg_sin[i] = mul_sin[i];
+      avg_cos[i] = mul_cos[i];
+      mean_sin += mul_sin[i];
+      mean_cos += mul_cos[i];
+    }
+    mean_sin /= samples_per_period;
+    mean_cos /= samples_per_period;
+    avg_sin[samples_per_period] = mean_sin;
+    avg_cos[samples_per_period] = mean_cos;
+    //First order low-pass
+    double a0 = 1-x,  b1 = x;
+    for (unsigned long n = samples_per_period+1; n<total_samples_; ++n) {
+      mag_sin[n] = a0*mul_sin[n] + b1*mag_sin[n-1];
+      mag_cos[n] = a0*mul_cos[n] + b1*mag_cos[n-1];
+      avg_sin[n] = avg_sin[n-1]
+        + (mul_sin[n] - mul_sin[n-samples_per_period])/samples_per_period;
+      avg_cos[n] = avg_cos[n-1]
+        + (mul_cos[n] - mul_cos[n-samples_per_period])/samples_per_period;
+    }
+    // // Forth order low-pass
+    // double a0 = pow2(pow2(1-x)), b1 = 4*x, b2 = -6*pow2(x), b3 = 4*x*x*x, b4 = - pow2(pow2(x));
+    // for (unsigned long n = 1; n<total_samples_; ++n) {
+    //   mag_sin[n] = a0*mul_sin[n] + b1*mag_sin[n-1]
+    //     + b2*mag_sin[n-2] + b3*mag_sin[n-3] + b4*mag_sin[n-4];
+    //   mag_cos[n] = a0*mul_cos[n] + b1*mag_cos[n-1]
+    //     + b2*mag_cos[n-2] + b3*mag_cos[n-3] + b4*mag_cos[n-4];
+    // }
+
+    // Stage 3: Evaluate magnitude and phase
+    std::vector<double> out_magnitude, out_phase;
+    out_magnitude.resize(total_samples_);
+    out_phase.resize(total_samples_);
+    for (unsigned long i = 0; i<total_samples_; ++i) {
+      out_magnitude[i] = std::sqrt(
+                                   pow2(avg_sin[i]) +
+                                   pow2(avg_cos[i])
+                                   );
+      out_phase[i] = std::atan(avg_sin[i]/avg_cos[i]);
+    }
+    
+    
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void LockIn::GenerateSinCos(const std::vector<double> &timesteps) {
+    omega_ = 2.0*PI_/period_ms_;
+    lockin_sin_.resize(total_samples_);
+    lockin_cos_.resize(total_samples_);
+    for (unsigned long i = 0; i < total_samples_; ++i) {
+      double time = timesteps[i];
+      lockin_sin_[i] = std::sin(omega_*time);
+      lockin_cos_[i] = std::cos(omega_*time);
+    }
+
+
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+
+}  // end of namespace mri

src/lock-in.h

@@ -1,9 +1,23 @@
 #ifndef SRC_LOCK_IN_H_
 #define SRC_LOCK_IN_H_
+#include <iostream>
+#include <vector>
+
 namespace mri {
   class LockIn {
   public:
-    void Analyse();
+    template<class T> inline T pow2(const T value) {return value*value;}
+    const double PI_=3.14159265358979323846;
+    void Analyse(const std::vector<double> &signal,double period_ms,
+                 const std::vector<double> &timesteps,
+                 int samples_per_period,
+                 std::vector<double> &avg_sin,
+                 std::vector<double> &avg_cos);
+  private:
+    std::vector<double> lockin_sin_, lockin_cos_;
+    double period_ms_, omega_;
+    unsigned long total_samples_;
+    void GenerateSinCos(const std::vector<double> &timesteps);
   };  // end of class LockIn
 }  // end of namespace mri
 #endif  // SRC_LOCK_IN_H_

src/main.cc

@@ -15,6 +15,7 @@
 const double PI_=3.14159265358979323846;
 
 int voxel_num = 3;
+//double phase_range = 0;
 double phase_range = PI_/2;
 double phase_init = 0.0;
 
@@ -23,7 +24,7 @@ double noise_ratio = 1e8;
 
 // B0=1.5T freq=64Mhz, period = 15.6 ns
 double period_ms = 15.6/1000/1000; //ms
-double total_periods = 10;
+double total_periods = 2000;
 int rf_samples_per_period = 20;
 // double period_ms = 10; //ms
 // double total_periods = 2;
@@ -32,6 +33,7 @@ int rf_samples_per_period = 20;
 //double T2s_scale = 0.01; //ms # need to be 10ms
 //double T2s_min = T2s_scale/1000.0; //
 double T2s_scale = period_ms*total_periods; //ms # need to be 10ms
+template<class T> inline T pow2(const T value) {return value*value;}
 
   
 int main(int argc, char *argv[]) {
@@ -39,16 +41,36 @@ int main(int argc, char *argv[]) {
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   try {
-    std::cout<<"Runing .... "<<std::endl;
+    //std::cout<<"Runing .... "<<std::endl;
     jade::SubPopulation sub_population;  // Optimizer.
     mri::GeneratorRF generator_rf;
     generator_rf.SetVoxels(voxel_num, phase_init, phase_range, T2s_scale);
     generator_rf.SetTimer(total_periods, period_ms, rf_samples_per_period);
     generator_rf.RandomizeVoxels();
     std::vector<double> signal_rf, magnitude_rf, phase_rf;
+    std::vector<double> mag_sin, mag_cos;
     generator_rf.EvaluateSignal(signal_rf);
+    generator_rf.EvaluateDemodulated(mag_sin,mag_cos);
     auto timesteps = generator_rf.GetTimesteps();
-    std::cout<< timesteps.size() <<std::endl;
+    mri::LockIn lock_in;
+    std::vector<double> avg_sin, avg_cos;
+    lock_in.Analyse(signal_rf, period_ms, timesteps,
+                    rf_samples_per_period, avg_sin, avg_cos);
+    for (unsigned long i = 0; i<timesteps.size(); ++i) {
+      if (i > static_cast<unsigned long>(rf_samples_per_period)*2) 
+        std::cout << timesteps[i] << "\t"
+                  << signal_rf[i] << "\t"
+        // << mul_sin[i] << "\t"
+        // << mul_cos[i] << "\t"
+          //         << std::sqrt(pow2(mag_sin[i])+pow2(mag_cos[i]))
+          // /std::sqrt(pow2(avg_sin[i])+pow2(avg_cos[i])) << "\t"
+                  << std::sqrt(pow2(mag_sin[i])+pow2(mag_cos[i])) << "\t"
+                  << 2.0*std::sqrt(pow2(avg_sin[i])+pow2(avg_cos[i])) << "\t"
+                  << avg_sin[i] << "\t"
+                  << avg_cos[i]
+                  << std::endl;
+    }
+    //std::cout<< timesteps.size() <<std::endl;
     // for (auto value:signal_rf){
     //   std::cout << value <<std::endl;
     // }