aboutsummaryrefslogtreecommitdiffstats log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
-rw-r--r--AlignSample.cpp50
1 files changed, 45 insertions, 5 deletions
 diff --git a/AlignSample.cpp b/AlignSample.cppindex 7de4459..0227d41 100644--- a/AlignSample.cpp+++ b/AlignSample.cpp@@ -233,7 +233,8 @@ static std::vector 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 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 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; }