7 #include "Correlation.h" 8 #include "EngaugeAssert.h" 16 m_signalA (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (2 * N - 1)))),
17 m_signalB (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (2 * N - 1)))),
18 m_outShifted (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (2 * N - 1)))),
19 m_outA (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (2 * N - 1)))),
20 m_outB (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (2 * N - 1)))),
21 m_out (static_cast<fftw_complex *> (fftw_malloc(sizeof(fftw_complex) * unsigned (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);
56 ENGAUGE_ASSERT (N > 0);
61 double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
62 for (i = 0; i < N; i++) {
64 sumMean1 += function1 [i];
65 sumMean2 += function2 [i];
66 max1 = qMax (max1, function1 [i]);
67 max2 = qMax (max2, function2 [i]);
79 double additiveNormalization1 = sumMean1 / N;
80 double additiveNormalization2 = sumMean2 / N;
81 double multiplicativeNormalization1 = 1.0 / max1;
82 double multiplicativeNormalization2 = 1.0 / max2;
86 for (i = 0; i < N - 1; i++) {
88 m_signalA [i] [0] = 0.0;
89 m_signalA [i] [1] = 0.0;
90 m_signalB [i + N] [0] = 0.0;
91 m_signalB [i + N] [1] = 0.0;
94 for (i = 0; i < N; i++) {
96 m_signalA [i + N - 1] [0] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
97 m_signalA [i + N - 1] [1] = 0.0;
98 m_signalB [i] [0] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
99 m_signalB [i] [1] = 0.0;
103 fftw_execute(m_planA);
104 fftw_execute(m_planB);
107 fftw_complex scale = {1.0/(2.0 * N - 1.0), 0.0};
108 for (i = 0; i < 2 * N - 1; i++) {
110 fftw_complex term1 = {m_outA [i] [0], m_outA [i] [1]};
111 fftw_complex term2 = {m_outB [i] [0], m_outB [i] [1] * -1.0};
112 fftw_complex term3 = {scale [0], scale [1]};
113 fftw_complex terms12 = {term1 [0] * term2 [0] - term1 [1] * term2 [1],
114 term1 [0] * term2 [1] + term1 [1] * term2 [0]};
115 m_out [i] [0] = terms12 [0] * term3 [0] - terms12 [1] * term3 [1];
116 m_out [i] [1] = terms12 [0] * term3 [1] + terms12 [1] * term3 [0];
119 fftw_execute(m_planX);
124 for (
int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
126 int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
127 fftw_complex shifted = {m_outShifted [i0AtCenter] [0], m_outShifted [i0AtCenter] [1]};
128 double corr = qSqrt (shifted [0] * shifted [0] + shifted [1] * shifted [1]);
130 if ((i0AtLeft == 0) || (corr > corrMax)) {
131 binStartMax = i0AtLeft;
136 correlations [i0AtLeft] = corr;
141 const double function1 [],
142 const double function2 [],
143 double &corrMax)
const 149 for (
int i = 0; i < N; i++) {
150 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.