Engauge Digitizer  2
Correlation.cpp
1 #include <complex.h>
2 #include "Correlation.h"
3 #include "EngaugeAssert.h"
4 #include "fftw3.h"
5 #include "Logger.h"
6 #include <QDebug>
7 #include <qmath.h>
8 
10  m_N (N),
11  m_signalA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
12  m_signalB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
13  m_outShifted ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
14  m_outA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
15  m_outB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
16  m_out ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1)))
17 {
18  m_planA = fftw_plan_dft_1d(2 * N - 1, m_signalA, m_outA, FFTW_FORWARD, FFTW_ESTIMATE);
19  m_planB = fftw_plan_dft_1d(2 * N - 1, m_signalB, m_outB, FFTW_FORWARD, FFTW_ESTIMATE);
20  m_planX = fftw_plan_dft_1d(2 * N - 1, m_out, m_outShifted, FFTW_BACKWARD, FFTW_ESTIMATE);
21 }
22 
23 Correlation::~Correlation()
24 {
25  fftw_destroy_plan(m_planA);
26  fftw_destroy_plan(m_planB);
27  fftw_destroy_plan(m_planX);
28 
29  fftw_free(m_signalA);
30  fftw_free(m_signalB);
31  fftw_free(m_outShifted);
32  fftw_free(m_out);
33  fftw_free(m_outA);
34  fftw_free(m_outB);
35 
36  fftw_cleanup();
37 }
38 
40  const double function1 [],
41  const double function2 [],
42  int &binStartMax,
43  double &corrMax) const
44 {
45 // LOG4CPP_DEBUG_S ((*mainCat)) << "Correlation::correlateWithShift";
46 
47  int i;
48 
49  ENGAUGE_ASSERT (N == m_N);
50 
51  // Normalize input functions so that:
52  // 1) mean is zero. This is used to compute an additive normalization constant
53  // 2) max value is 1. This is used to compute a multiplicative normalization constant
54  double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
55  for (i = 0; i < N; i++) {
56 
57  sumMean1 += function1 [i];
58  sumMean2 += function2 [i];
59  max1 = qMax (max1, function1 [i]);
60  max2 = qMax (max2, function2 [i]);
61 
62  }
63 
64  double additiveNormalization1 = sumMean1 / N;
65  double additiveNormalization2 = sumMean2 / N;
66  double multiplicativeNormalization1 = 1.0 / max1;
67  double multiplicativeNormalization2 = 1.0 / max2;
68 
69  // Load length N functions into length 2N+1 arrays, padding with zeros before for the first
70  // array, and with zeros after for the second array
71  for (i = 0; i < N - 1; i++) {
72 
73  m_signalA [i] = 0.0;
74  m_signalB [i + N] = 0.0;
75 
76  }
77  for (i = 0; i < N; i++) {
78 
79  m_signalA [i + N - 1] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
80  m_signalB [i] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
81 
82  }
83 
84  fftw_execute(m_planA);
85  fftw_execute(m_planB);
86 
87  // Correlation in frequency space
88  fftw_complex scale = 1.0/(2.0 * N - 1.0);
89  for (i = 0; i < 2 * N - 1; i++) {
90  m_out[i] = m_outA[i] * conj(m_outB[i]) * scale;
91  }
92 
93  fftw_execute(m_planX);
94 
95  // Search for highest correlation. We have to account for the shift in the index. Specifically,
96  // 0 to N was mapped to the second half of the array that is 0 to 2 * N - 1
97  corrMax = 0.0;
98  for (int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
99 
100  int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
101  fftw_complex shifted = m_outShifted [i0AtCenter];
102  double corr = qSqrt (creal (shifted) * creal (shifted) + cimag (shifted) * cimag (shifted));
103 
104  if ((i0AtLeft == 0) || (corr > corrMax)) {
105  binStartMax = i0AtLeft;
106  corrMax = corr;
107  }
108  }
109 }
110 
112  const double function1 [],
113  const double function2 [],
114  double &corrMax) const
115 {
116 // LOG4CPP_DEBUG_S ((*mainCat)) << "Correlation::correlateWithoutShift";
117 
118  corrMax = 0.0;
119 
120  for (int i = 0; i < N; i++) {
121  corrMax += function1 [i] * function2 [i];
122  }
123 }
void correlateWithShift(int N, const double function1[], const double function2[], int &binStartMax, double &corrMax) const
Return the shift in function1 that best aligns that function with function2.
Definition: Correlation.cpp:39
Correlation(int N)
Single constructor. Slow memory allocations are done once and then reused repeatedly.
Definition: Correlation.cpp:9
void correlateWithoutShift(int N, const double function1[], const double function2[], double &corrMax) const
Return the correlation of the two functions, without any shift.