1 #include "Correlation.h"
2 #include "EngaugeAssert.h"
10 m_signalA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
11 m_signalB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
12 m_outShifted ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
13 m_outA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
14 m_outB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
15 m_out ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1)))
17 m_planA = fftw_plan_dft_1d(2 * N - 1, m_signalA, m_outA, FFTW_FORWARD, FFTW_ESTIMATE);
18 m_planB = fftw_plan_dft_1d(2 * N - 1, m_signalB, m_outB, FFTW_FORWARD, FFTW_ESTIMATE);
19 m_planX = fftw_plan_dft_1d(2 * N - 1, m_out, m_outShifted, FFTW_BACKWARD, FFTW_ESTIMATE);
22 Correlation::~Correlation()
24 fftw_destroy_plan(m_planA);
25 fftw_destroy_plan(m_planB);
26 fftw_destroy_plan(m_planX);
30 fftw_free(m_outShifted);
39 const double function1 [],
40 const double function2 [],
43 double correlations [])
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++) {
73 m_signalA [i] [0] = 0.0;
74 m_signalA [i] [1] = 0.0;
75 m_signalB [i + N] [0] = 0.0;
76 m_signalB [i + N] [1] = 0.0;
79 for (i = 0; i < N; i++) {
81 m_signalA [i + N - 1] [0] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
82 m_signalA [i + N - 1] [1] = 0.0;
83 m_signalB [i] [0] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
84 m_signalB [i] [1] = 0.0;
88 fftw_execute(m_planA);
89 fftw_execute(m_planB);
92 fftw_complex scale = {1.0/(2.0 * N - 1.0), 0.0};
93 for (i = 0; i < 2 * N - 1; i++) {
95 fftw_complex term1 = {m_outA [i] [0], m_outA [i] [1]};
96 fftw_complex term2 = {m_outB [i] [0], m_outB [i] [1] * -1.0};
97 fftw_complex term3 = {scale [0], scale [1]};
98 fftw_complex terms12 = {term1 [0] * term2 [0] - term1 [1] * term2 [1],
99 term1 [0] * term2 [1] + term1 [1] * term2 [0]};
100 m_out [i] [0] = terms12 [0] * term3 [0] - terms12 [1] * term3 [1];
101 m_out [i] [1] = terms12 [0] * term3 [1] + terms12 [1] * term3 [0];
104 fftw_execute(m_planX);
109 for (
int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
111 int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
112 fftw_complex shifted = {m_outShifted [i0AtCenter] [0], m_outShifted [i0AtCenter] [1]};
113 double corr = qSqrt (shifted [0] * shifted [0] + shifted [1] * shifted [1]);
115 if ((i0AtLeft == 0) || (corr > corrMax)) {
116 binStartMax = i0AtLeft;
121 correlations [i0AtLeft] = corr;
126 const double function1 [],
127 const double function2 [],
128 double &corrMax)
const
134 for (
int i = 0; i < N; i++) {
135 corrMax += function1 [i] * function2 [i];
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.
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.