2 #include "Correlation.h"
3 #include "EngaugeAssert.h"
11 m_signalA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
12 m_signalB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
13 m_outShifted ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
14 m_outA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
15 m_outB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
16 m_out ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1)))
18 m_planA = fftw_plan_dft_1d(2 * N - 1, m_signalA, m_outA, FFTW_FORWARD, FFTW_ESTIMATE);
19 m_planB = fftw_plan_dft_1d(2 * N - 1, m_signalB, m_outB, FFTW_FORWARD, FFTW_ESTIMATE);
20 m_planX = fftw_plan_dft_1d(2 * N - 1, m_out, m_outShifted, FFTW_BACKWARD, FFTW_ESTIMATE);
23 Correlation::~Correlation()
25 fftw_destroy_plan(m_planA);
26 fftw_destroy_plan(m_planB);
27 fftw_destroy_plan(m_planX);
31 fftw_free(m_outShifted);
40 const double function1 [],
41 const double function2 [],
43 double &corrMax)
const
49 ENGAUGE_ASSERT (N == m_N);
54 double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
55 for (i = 0; i < N; i++) {
57 sumMean1 += function1 [i];
58 sumMean2 += function2 [i];
59 max1 = qMax (max1, function1 [i]);
60 max2 = qMax (max2, function2 [i]);
64 double additiveNormalization1 = sumMean1 / N;
65 double additiveNormalization2 = sumMean2 / N;
66 double multiplicativeNormalization1 = 1.0 / max1;
67 double multiplicativeNormalization2 = 1.0 / max2;
71 for (i = 0; i < N - 1; i++) {
74 m_signalB [i + N] = 0.0;
77 for (i = 0; i < N; i++) {
79 m_signalA [i + N - 1] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
80 m_signalB [i] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
84 fftw_execute(m_planA);
85 fftw_execute(m_planB);
88 fftw_complex scale = 1.0/(2.0 * N - 1.0);
89 for (i = 0; i < 2 * N - 1; i++) {
90 m_out[i] = m_outA[i] * conj(m_outB[i]) * scale;
93 fftw_execute(m_planX);
98 for (
int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
100 int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
101 fftw_complex shifted = m_outShifted [i0AtCenter];
102 double corr = qSqrt (creal (shifted) * creal (shifted) + cimag (shifted) * cimag (shifted));
104 if ((i0AtLeft == 0) || (corr > corrMax)) {
105 binStartMax = i0AtLeft;
112 const double function1 [],
113 const double function2 [],
114 double &corrMax)
const
120 for (
int i = 0; i < N; i++) {
121 corrMax += function1 [i] * function2 [i];
void correlateWithShift(int N, const double function1[], const double function2[], int &binStartMax, double &corrMax) const
Return the shift in function1 that best aligns that function with function2.
Correlation(int N)
Single constructor. Slow memory allocations are done once and then reused repeatedly.
void correlateWithoutShift(int N, const double function1[], const double function2[], double &corrMax) const
Return the correlation of the two functions, without any shift.