Engauge Digitizer 2
Loading...
Searching...
No Matches
Correlation.cpp
Go to the documentation of this file.
1/******************************************************************************************************
2 * (C) 2014 markummitchell@github.com. This file is part of Engauge Digitizer, which is released *
3 * under GNU General Public License version 2 (GPLv2) or (at your option) any later version. See file *
4 * LICENSE or go to gnu.org/licenses for details. Distribution requires prior written permission. *
5 ******************************************************************************************************/
6
7#include "Correlation.h"
8#include "EngaugeAssert.h"
9#include "fftw3.h"
10#include "Logger.h"
11#include <QDebug>
12#include <qmath.h>
13
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}
27
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}
43
45 const double function1 [],
46 const double function2 [],
47 int &binStartMax,
48 double &corrMax,
49 double correlations []) const
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}
139
141 const double function1 [],
142 const double function2 [],
143 double &corrMax) const
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}
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT if defined(QT_NO_DEBUG) && !defined(QT_FORCE_ASSERTS) define ENGAUGE...
Definition: EngaugeAssert.h:20
void correlateWithoutShift(int N, const double function1[], const double function2[], double &corrMax) const
Return the correlation of the two functions, without any shift.
Correlation(int N)
Single constructor. Slow memory allocations are done once and then reused repeatedly.
Definition: Correlation.cpp:14
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.
Definition: Correlation.cpp:44