Engauge Digitizer  2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
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. More...
 
 ~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. More...
 
void correlateWithoutShift (int N, const double function1[], const double function2[], double &corrMax) const
 Return the correlation of the two functions, without any shift. More...
 

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 ( 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 ( )

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

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 }
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 if defined(QT_NO_DEBUG) &amp;&amp; !defined(QT_FORCE_ASSERTS) define ENGAUGE...
Definition: EngaugeAssert.h:20

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