correlator transform utilities

lib/Physics/CorrelatorFitter.cpp

@@ -227,6 +227,77 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr,
     return init;
 }
 
+/******************************************************************************
+ *                         Correlator utilities                               *
+ ******************************************************************************/
+DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
+{
+    if (ts != 0)
+    {
+        const Index nt  = c[central].rows();
+        DMatSample  buf = c;
+
+        FOR_STAT_ARRAY(buf, s)
+        {
+            for (Index t = 0; t < nt; ++t)
+            {
+                buf[s]((t - ts + nt)%nt) = c[s](t);
+            }
+        }
+
+        return buf;
+    }
+    else
+    {
+        return c;
+    }
+}
+
+DMatSample CorrelatorUtils::fold(const DMatSample &c)
+{
+    const Index nt  = c[central].rows();
+    DMatSample  buf = c;
+
+    FOR_STAT_ARRAY(buf, s)
+    {
+        for (Index t = 0; t < nt; ++t)
+        {
+            buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
+        }
+    }
+
+    return buf;
+}
+
+DMatSample CorrelatorUtils::fourierTransform(const DMatSample &c, FFT &fft, 
+                                             const unsigned int dir)
+{
+    const Index nSample   = c.size();
+    const Index nt        = c[central].rows();
+    bool        isComplex = (c[central].cols() > 1);
+    CMatSample  buf(nSample, nt, 1);
+    DMatSample  out(nSample, nt, 2);
+
+    fft.resize(nt);
+    FOR_STAT_ARRAY(buf, s)
+    {
+        buf[s].real() = c[s].col(0);
+        if (isComplex)
+        {
+            buf[s].imag() = c[s].col(1);
+        }
+        else
+        {
+            buf[s].imag() = DVec::Constant(nt, 0.);
+        }
+        fft(buf[s], dir);
+        out[s].col(0) = buf[s].real();
+        out[s].col(1) = buf[s].imag();
+    }
+
+    return out;    
+}
+
 /******************************************************************************
  *                      CorrelatorFitter implementation                       *
  ******************************************************************************/

lib/Physics/CorrelatorFitter.hpp

@@ -22,6 +22,7 @@
 
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Functional/Model.hpp>
+#include <LatAnalyze/Numerical/FFT.hpp>
 #include <LatAnalyze/Statistics/XYSampleData.hpp>
 
 BEGIN_LATAN_NAMESPACE
@@ -49,6 +50,17 @@ namespace CorrelatorModels
     DVec        parameterGuess(const DMatSample &corr, const ModelPar par);
 };
 
+/******************************************************************************
+ *                         Correlator utilities                               *
+ ******************************************************************************/
+namespace CorrelatorUtils
+{
+    DMatSample shift(const DMatSample &c, const Index ts);
+    DMatSample fold(const DMatSample &c);
+    DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
+                                const unsigned int dir = FFT::Forward);
+};
+
 /******************************************************************************
  *                       Correlator fit utility class                         *
  ******************************************************************************/

physics/2pt-fit.cpp

@@ -105,28 +105,14 @@ int main(int argc, char *argv[])
     DMatSample tmp, corr;
     Index      nSample, nt;
     
-    tmp     = Io::load<DMatSample>(corrFileName);
+    corr    = Io::load<DMatSample>(corrFileName);
-    nSample = tmp.size();
+    nSample = corr.size();
-    nt      = tmp[central].rows();
+    nt      = corr[central].rows();
-    tmp     = tmp.block(0, 0, nt, 1);
+    corr    = corr.block(0, 0, nt, 1);
-    corr    = tmp;
+    corr    = CorrelatorUtils::shift(corr, shift);
-    FOR_STAT_ARRAY(corr, s)
-    {
-        for (Index t = 0; t < nt; ++t)
-        {
-            corr[s]((t - shift + nt)%nt) = tmp[s](t);
-        }
-    }
     if (fold)
     {
-        tmp = corr;
+        corr = CorrelatorUtils::fold(corr);
-        FOR_STAT_ARRAY(corr, s)
-        {
-            for (Index t = 0; t < nt; ++t)
-            {
-                corr[s](t) = 0.5*(tmp[s](t) + tmp[s]((nt - t) % nt));
-            }
-        }
     }
     
     // make model //////////////////////////////////////////////////////////////