7 #include "Correlation.h" 8 #include "EngaugeAssert.h" 16 m_signalA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
17 m_signalB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
18 m_outShifted ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
19 m_outA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
20 m_outB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
21 m_out ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1)))
23 m_planA = fftw_plan_dft_1d(2 * N - 1, m_signalA, m_outA, FFTW_FORWARD, FFTW_ESTIMATE);
24 m_planB = fftw_plan_dft_1d(2 * N - 1, m_signalB, m_outB, FFTW_FORWARD, FFTW_ESTIMATE);
25 m_planX = fftw_plan_dft_1d(2 * N - 1, m_out, m_outShifted, FFTW_BACKWARD, FFTW_ESTIMATE);
28 Correlation::~Correlation()
30 fftw_destroy_plan(m_planA);
31 fftw_destroy_plan(m_planB);
32 fftw_destroy_plan(m_planX);
36 fftw_free(m_outShifted);
45 const double function1 [],
46 const double function2 [],
49 double correlations [])
const 55 ENGAUGE_ASSERT (N == m_N);
60 double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
61 for (i = 0; i < N; i++) {
63 sumMean1 += function1 [i];
64 sumMean2 += function2 [i];
65 max1 = qMax (max1, function1 [i]);
66 max2 = qMax (max2, function2 [i]);
70 double additiveNormalization1 = sumMean1 / N;
71 double additiveNormalization2 = sumMean2 / N;
72 double multiplicativeNormalization1 = 1.0 / max1;
73 double multiplicativeNormalization2 = 1.0 / max2;
77 for (i = 0; i < N - 1; i++) {
79 m_signalA [i] [0] = 0.0;
80 m_signalA [i] [1] = 0.0;
81 m_signalB [i + N] [0] = 0.0;
82 m_signalB [i + N] [1] = 0.0;
85 for (i = 0; i < N; i++) {
87 m_signalA [i + N - 1] [0] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
88 m_signalA [i + N - 1] [1] = 0.0;
89 m_signalB [i] [0] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
90 m_signalB [i] [1] = 0.0;
94 fftw_execute(m_planA);
95 fftw_execute(m_planB);
98 fftw_complex scale = {1.0/(2.0 * N - 1.0), 0.0};
99 for (i = 0; i < 2 * N - 1; i++) {
101 fftw_complex term1 = {m_outA [i] [0], m_outA [i] [1]};
102 fftw_complex term2 = {m_outB [i] [0], m_outB [i] [1] * -1.0};
103 fftw_complex term3 = {scale [0], scale [1]};
104 fftw_complex terms12 = {term1 [0] * term2 [0] - term1 [1] * term2 [1],
105 term1 [0] * term2 [1] + term1 [1] * term2 [0]};
106 m_out [i] [0] = terms12 [0] * term3 [0] - terms12 [1] * term3 [1];
107 m_out [i] [1] = terms12 [0] * term3 [1] + terms12 [1] * term3 [0];
110 fftw_execute(m_planX);
115 for (
int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
117 int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
118 fftw_complex shifted = {m_outShifted [i0AtCenter] [0], m_outShifted [i0AtCenter] [1]};
119 double corr = qSqrt (shifted [0] * shifted [0] + shifted [1] * shifted [1]);
121 if ((i0AtLeft == 0) || (corr > corrMax)) {
122 binStartMax = i0AtLeft;
127 correlations [i0AtLeft] = corr;
132 const double function1 [],
133 const double function2 [],
134 double &corrMax)
const 140 for (
int i = 0; i < N; i++) {
141 corrMax += function1 [i] * function2 [i];
void correlateWithoutShift(int N, const double function1 [], const double function2 [], double &corrMax) const
Return the correlation of the two functions, without any shift.
Correlation(int N)
Single constructor. Slow memory allocations are done once and then reused repeatedly.
void correlateWithShift(int N, const double function1 [], const double function2 [], int &binStartMax, double &corrMax, double correlations []) const
Return the shift in function1 that best aligns that function with function2.