Engauge Digitizer 2
Loading...
Searching...
No Matches
Correlation Class Reference

Fast cross correlation between two functions. More...

#include <Correlation.h>

Collaboration diagram for Correlation:
Collaboration graph

Public Member Functions

 Correlation (int N)
 Single constructor. Slow memory allocations are done once and then reused repeatedly.
 ~Correlation ()
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.
void correlateWithoutShift (int N, const double function1[], const double function2[], double &corrMax) const
 Return the correlation of the two functions, without any shift.

Detailed Description

Fast cross correlation between two functions.

We do not use complex.h along with fftw3.h since then the complex numbers will be native, which would then require platform-dependent code

Definition at line 14 of file Correlation.h.

Constructor & Destructor Documentation

◆ Correlation()

Correlation::Correlation ( int N)

Single constructor. Slow memory allocations are done once and then reused repeatedly.

Definition at line 14 of file Correlation.cpp.

14 :
15 m_N (N),
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))))
22{
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);
26}

◆ ~Correlation()

Correlation::~Correlation ( )

Definition at line 28 of file Correlation.cpp.

29{
30 fftw_destroy_plan(m_planA);
31 fftw_destroy_plan(m_planB);
32 fftw_destroy_plan(m_planX);
33
34 fftw_free(m_signalA);
35 fftw_free(m_signalB);
36 fftw_free(m_outShifted);
37 fftw_free(m_out);
38 fftw_free(m_outA);
39 fftw_free(m_outB);
40
41 fftw_cleanup();
42}

Member Function Documentation

◆ correlateWithoutShift()

void Correlation::correlateWithoutShift ( int N,
const double function1[],
const double function2[],
double & corrMax ) const

Return the correlation of the two functions, without any shift.

The functions are normalized internally.

Definition at line 140 of file Correlation.cpp.

144{
145// LOG4CPP_DEBUG_S ((*mainCat)) << "Correlation::correlateWithoutShift";
146
147 corrMax = 0.0;
148
149 for (int i = 0; i < N; i++) {
150 corrMax += function1 [i] * function2 [i];
151 }
152}

◆ correlateWithShift()

void Correlation::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.

The functions are normalized internally. The correlations vector, as a function of shift, is returned for logging

Definition at line 44 of file Correlation.cpp.

50{
51 // LOG4CPP_DEBUG_S ((*mainCat)) << "Correlation::correlateWithShift";
52
53 int i;
54
55 ENGAUGE_ASSERT (N == m_N);
56 ENGAUGE_ASSERT (N > 0); // Prevent divide by zero errors for additiveNormalization* and scale
57
58 // Normalize input functions so that:
59 // 1) mean is zero. This is used to compute an additive normalization constant
60 // 2) max value is 1. This is used to compute a multiplicative normalization constant
61 double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
62 for (i = 0; i < N; i++) {
63
64 sumMean1 += function1 [i];
65 sumMean2 += function2 [i];
66 max1 = qMax (max1, function1 [i]);
67 max2 = qMax (max2, function2 [i]);
68
69 }
70
71 // Handle all-zero data
72 if (max1 == 0.0) {
73 max1 = 1.0;
74 }
75 if (max2 == 0.0) {
76 max2 = 1.0;
77 }
78
79 double additiveNormalization1 = sumMean1 / N;
80 double additiveNormalization2 = sumMean2 / N;
81 double multiplicativeNormalization1 = 1.0 / max1;
82 double multiplicativeNormalization2 = 1.0 / max2;
83
84 // Load length N functions into length 2N+1 arrays, padding with zeros before for the first
85 // array, and with zeros after for the second array
86 for (i = 0; i < N - 1; i++) {
87
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;
92
93 }
94 for (i = 0; i < N; i++) {
95
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;
100
101 }
102
103 fftw_execute(m_planA);
104 fftw_execute(m_planB);
105
106 // Correlation in frequency space
107 fftw_complex scale = {1.0/(2.0 * N - 1.0), 0.0};
108 for (i = 0; i < 2 * N - 1; i++) {
109 // Multiple m_outA [i] * conj (m_outB) * scale
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];
117 }
118
119 fftw_execute(m_planX);
120
121 // Search for highest correlation. We have to account for the shift in the index. Specifically,
122 // 0 to N was mapped to the second half of the array that is 0 to 2 * N - 1
123 corrMax = 0.0;
124 for (int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
125
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]);
129
130 if ((i0AtLeft == 0) || (corr > corrMax)) {
131 binStartMax = i0AtLeft;
132 corrMax = corr;
133 }
134
135 // Save for, if enabled, external logging
136 correlations [i0AtLeft] = corr;
137 }
138}
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT.

The documentation for this class was generated from the following files: