working implementation with basic functionalities

Cpp/lib/waveform.cpp

@@ -152,19 +152,21 @@ void ArrayWaveform::ArrayParam::printParam() {
     std::cout << "amplitudes:" << "\n" << this->amplitudes << "\n";
 }
 
-EigenVectorXi16 ArrayWaveform::getStaticWaveform(const ArrayParam& param) {
-    auto minSampleLen = 2 * param.samplingRate
-        / std::gcd(param.samplingRate, param.freqResolution);
-    Eigen::VectorXd t = Eigen::VectorXd::LinSpaced(
-        minSampleLen,
-        0,
-        minSampleLen - 1
-    ) / param.samplingRate;
-    Eigen::MatrixXd wfmMatrix = Eigen::sin(
-        ((param.freqTones * t.transpose() * M_PI * 2).array().colwise()
-        + param.phases.array())
-    ).array().colwise() * param.amplitudes.array();
-    return wfmMatrix.colwise().sum().cast<int16_t>();
+ulong ArrayWaveform::getMinSampleLen(
+    ulong samplingRate,
+    ulong freqResolution
+) {
+    return 2 * samplingRate / std::gcd(samplingRate, freqResolution);
+}
+
+ulong ArrayWaveform::getSampleLen(
+    double tau,
+    ulong samplingRate,
+    ulong freqResolution
+) {
+    ulong sampleLen = ulong(tau * samplingRate);
+    ulong minLen = getMinSampleLen(samplingRate, freqResolution);
+    return sampleLen % minLen == 0 ? sampleLen : (sampleLen / minLen + 1) * minLen;
 }
 
 double ArrayWaveform::setStaticSegment(
@@ -173,11 +175,11 @@ double ArrayWaveform::setStaticSegment(
     double initPhase,
     double amp
 ) {
-    Eigen::Ref<Eigen::ArrayXd> t = timeSeries.array();
-    t -= t(0);
-    double last = 2*M_PI * f * (t(t.size() - 1) + t(1)) + initPhase;
-    timeSeries.array() = amp * Eigen::sin(2*M_PI * f * t + initPhase);
-    return last;
+  Eigen::Ref<Eigen::ArrayXd> t = timeSeries.array();
+  t -= t(0);
+  double nextPhase = 2 * M_PI * f * (t(t.size() - 1) + t(1)) + initPhase;
+  timeSeries.array() = amp * Eigen::sin(2 * M_PI * f * t + initPhase);
+  return nextPhase;
 }
 
 double ArrayWaveform::setMovingSegment(
@@ -195,21 +197,101 @@ double ArrayWaveform::setMovingSegment(
     unsigned int midIdx = int(sampleLen / 2.0);
     Eigen::Ref<Eigen::ArrayXd> segl = timeSeries(Eigen::seq(0, midIdx));
     Eigen::Ref<Eigen::ArrayXd> segr = timeSeries(Eigen::seq(midIdx+1, sampleLen-1));
-    segr -= segl(0) + tau / 2;
+    segr -= segl(0) + tau / 2;  // it may seem more intuitive to use
+                                // segr -= segr(0) here, but doing so
+                                // will cause phase jump
     segl -= segl(0);
     segl = initPhase
         + 2 * M_PI * fInit * segl
         + accel / 6 * Eigen::pow(segl,3);
-    // segr.array() -= tau / 2;
     segr = segl(segl.size()-1)
         + (
             2 * M_PI * fInit + accel / 2 * std::pow((tau / 2),2)
         ) * segr
         + accel / 2 * tau / 2 * Eigen::pow(segr, 2)
         - accel / 6 * Eigen::pow(segr, 3);
-    double last = segr(segr.size() - 1) + dt * 2*M_PI * fFinal;
+    double nextPhase = segr(segr.size() - 1) + dt * 2*M_PI * fFinal;
     timeSeries = amp * Eigen::sin(timeSeries.array());
-    return last;
+    return nextPhase;
+}
+
+EigenVectorXi16 ArrayWaveform::getStaticWaveform(const ArrayParam& param) {
+    auto minSampleLen = 2 * param.samplingRate
+        / std::gcd(param.samplingRate, param.freqResolution);
+    Eigen::VectorXd t = Eigen::VectorXd::LinSpaced(
+        minSampleLen,
+        0,
+        minSampleLen - 1
+    ) / param.samplingRate;
+    Eigen::MatrixXd wfmMatrix = Eigen::sin(
+        ((param.freqTones * t.transpose() * M_PI * 2).array().colwise()
+        + param.phases.array())
+    ).array().colwise() * param.amplitudes.array();
+    return wfmMatrix.colwise().sum().cast<int16_t>();
+}
+
+EigenVectorXi16 ArrayWaveform::getTrickWaveform(
+    const ArrayParam& param,
+    std::set<int> siteIndex,
+    double df,
+    double tauMove,
+    double tauStay
+) {
+    if ((tauMove == 0 and tauStay == 0)
+        or siteIndex.empty()) {
+        return getStaticWaveform(param);
+    }
+    if (tauMove != 0 and tauStay != 0) {
+        tauStay = std::ceil((tauMove + tauStay) * df) / df - tauMove;
+    } else if (tauMove != 0) {
+        tauMove = std::ceil(tauMove * df) / df;
+    } else {
+        tauStay = std::ceil(tauStay * df) / df;
+    }
+    auto tauTotal = tauMove * 2 + tauStay;
+    ulong sampleLen = getSampleLen(tauTotal, param.samplingRate, param.freqResolution);
+    ulong moveLen = tauMove * param.samplingRate;
+    ulong stayLen = tauStay * param.samplingRate;
+    auto idxSeg0 = moveLen;
+    auto idxSeg1 = moveLen + stayLen;
+    auto idxSeg2 = moveLen + stayLen + moveLen;
+    auto t = Eigen::ArrayXd::LinSpaced(sampleLen, 0, sampleLen - 1) / param.samplingRate;
+    Eigen::MatrixXd wfmMat(sampleLen, param.freqTones.size()); // site wfm stored column by column
+    for (auto i = 0; i < param.freqTones.size(); i++) {
+        wfmMat.col(i) = t;
+        Eigen::Ref<Eigen::ArrayXd> siteWfm = wfmMat.col(i);
+        auto fInit = param.freqTones[i];
+        auto fFinal = fInit + df;
+        auto phi = param.phases[i];
+        auto amp = param.amplitudes[i];
+        if (siteIndex.contains(i)) {
+            if (tauMove != 0) {
+                phi = setMovingSegment(
+                    siteWfm(Eigen::seq(0, idxSeg0 - 1)),
+                    fInit, fFinal, phi, amp
+                );
+            }
+            if (tauStay != 0) {
+                phi = setStaticSegment(
+                    siteWfm(Eigen::seq(idxSeg0, idxSeg1 - 1)),
+                    fFinal, phi, amp
+                );
+            }
+            if (tauMove != 0) {
+                phi = setMovingSegment(
+                    siteWfm(Eigen::seq(idxSeg1, idxSeg2 - 1)),
+                    fFinal, fInit, phi, amp
+                );
+            }
+            setStaticSegment(
+                siteWfm(Eigen::seq(idxSeg2, siteWfm.size() - 1)),
+                fInit, phi, amp
+            );
+        } else {
+            setStaticSegment(siteWfm, fInit, phi, amp);
+        }
+    }
+    return wfmMat.rowwise().sum().cast<int16_t> ();
 }
 
 void ArrayWaveform::saveCSV(std::string fileName, EigenVectorXi16 wfm) {

Cpp/lib/waveform.h

 #pragma once
 #include <Eigen/Dense>
 #include <string>
+#include <set>
 typedef Eigen::Vector<int16_t, Eigen::Dynamic> EigenVectorXi16;
 
 namespace ArrayWaveform {
@@ -35,8 +36,14 @@ namespace ArrayWaveform {
         void loadParam(std::string fileName);
         void printParam();
     };
-    
-    EigenVectorXi16 getStaticWaveform(const ArrayParam& param);
+
+    ulong getMinSampleLen(ulong samplingRate, ulong freqResolution);
+    ulong getSampleLen(
+        double tau,
+        ulong samplingRate,
+        ulong freqResolution
+    );
+
     double setStaticSegment(
         Eigen::Ref<Eigen::VectorXd> timeSeries,
         double f,
@@ -50,6 +57,16 @@ namespace ArrayWaveform {
         double initPhase,
         double amp
     );
+
+    EigenVectorXi16 getStaticWaveform(const ArrayParam& param);
+    EigenVectorXi16 getTrickWaveform(
+        const ArrayParam& param,
+        std::set<int> siteIndex,
+        double df,
+        double tauMove=0,
+        double tauStay=0
+    );
+
     void saveCSV(std::string fileName, EigenVectorXi16 wfm);
 };