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);
45 const double function1 [],
46 const double function2 [],
49 double correlations [])
const
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;
void correlateWithoutShift(int N, const double function1[], const double function2[], double &corrMax) const
Return the correlation of the two functions, without any shift.
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.