Add comments to subsample alignerHEADmaster

1 files changed, 45 insertions, 5 deletions

diff --git a/AlignSample.cpp b/AlignSample.cpp
index 7de4459..0227d41 100644
--- a/AlignSample.cpp
+++ b/AlignSample.cpp
@@ -233,7 +233,8 @@ static std::vector<double> gen_omega(size_t length)
 }
 
 /* Find subsample delay in signal versus the reference signal ref and
- * correct it
+ * correct it. This assumes that the offset between signal and ref
+ * is less than one sample.
  */
 static vec_cf align_subsample(vec_cf& signal, vec_cf& ref)
 {
@@ -246,19 +247,38 @@ static vec_cf align_subsample(vec_cf& signal, vec_cf& ref)
     fftw::plan plan_ifft(N, rotate_vec, corr_sig, FFTW_BACKWARD);
     plan.execute();
 
-    const auto omega = gen_omega(N);
+    const std::vector<double> omega = gen_omega(N);
+
+    // Calculate the correlation for a lag of tau, which
+    // must be between -1 and 1 samples
     auto correlate_point = [&](float tau) {
+
+        // A subsample offset between two signals corresponds, in the frequency
+        // domain, to a linearly increasing phase shift, whose slope
+        // corresponds to the delay.
+        //
+        // Here, we build this phase shift in rotate_vec, and multiply it with
+        // our signal.
         for (size_t i = 0; i < N; i++) {
             const complexf angle(0, tau * i);
-            rotate_vec[i] = std::exp(angle);
+            rotate_vec[i] = std::exp(angle); // e^(j*tau*i)
         }
+        // zero-frequency is rotate_vec[0], so rotate_vec[N/2] is the
+        // bin corresponding to the [-1, 1, -1, 1, ...] time signal, which
+        // is both the maximum positive and negative frequency.
+        // I don't remember why we handle it differently.
         rotate_vec[N/2] = cos(M_PI * tau);
 
         for (size_t i = 0; i < N; i++) {
             rotate_vec[i] *= fft_sig[i];
         }
 
-        plan_ifft.execute();
+        plan_ifft.execute(); // corr_sig = IFFT(rotate_vec)
+
+        // Calculate correlation only the real part of the signal:
+        // The following implements the elementwise vector operation
+        //   corr_real = real(corr_sig) * real(ref)
+        // in an awkward way.
 
         std::vector<float> corr_real(N);
         std::transform(corr_sig.begin(), corr_sig.end(), corr_real.begin(),
@@ -268,18 +288,38 @@ static vec_cf align_subsample(vec_cf& signal, vec_cf& ref)
             complexf f = corr_real[i] * ref[i];
             corr_real[i] = f.real();
         }
+
+        // TODO: replace by the following and verify it's equivalent
+        /*
+        for (size_t i = 0; i < N; i++) {
+            corr_real[i] = corr_sig[i].real() * ref[i].real();
+        }
+        */
+
+        // TODO why do we only look at the real part? Because it's faster than
+        // a complex cross-correlation? Clarify!
+        // Compare with
+        /*
+        for (size_t i = 0; i < N; i++) {
+            corr_real[i] = (corr_sig[i] * std::conj(ref[i])).real();
+        }
+        */
+
+        // The correlation is the sum:
         return std::accumulate(corr_real.begin(), corr_real.end(), 0.0f);
     };
 
     const float start_arg = -1;
     const float end_arg = 1;
 
+    // Let boost find the tau value with the minimal correlation
     auto tau =
         boost::math::tools::brent_find_minima(
                 correlate_point, start_arg, end_arg, 1000).first;
 
     fprintf(stderr, "Fractional delay is %f\n", tau);
 
+    // Prepare rotate_vec = fft_sig with rotated phase
     for (size_t i = 0; i < N; i++) {
         const complexf angle(0, tau * i);
         rotate_vec[i] = std::exp(angle);
@@ -289,7 +329,7 @@ static vec_cf align_subsample(vec_cf& signal, vec_cf& ref)
         rotate_vec[i] *= fft_sig[i];
     }
 
-    plan_ifft.execute();
+    plan_ifft.execute(); // corr_sig = IFFT(rotate_vec)
 
     return corr_sig;
 }