00001
00002
00003
00004
00005
00006
00007 #include "Correlation.h"
00008 #include "EngaugeAssert.h"
00009 #include "fftw3.h"
00010 #include "Logger.h"
00011 #include <QDebug>
00012 #include <qmath.h>
00013
00014 Correlation::Correlation(int N) :
00015 m_N (N),
00016 m_signalA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
00017 m_signalB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
00018 m_outShifted ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
00019 m_outA ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
00020 m_outB ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1))),
00021 m_out ((fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1)))
00022 {
00023 m_planA = fftw_plan_dft_1d(2 * N - 1, m_signalA, m_outA, FFTW_FORWARD, FFTW_ESTIMATE);
00024 m_planB = fftw_plan_dft_1d(2 * N - 1, m_signalB, m_outB, FFTW_FORWARD, FFTW_ESTIMATE);
00025 m_planX = fftw_plan_dft_1d(2 * N - 1, m_out, m_outShifted, FFTW_BACKWARD, FFTW_ESTIMATE);
00026 }
00027
00028 Correlation::~Correlation()
00029 {
00030 fftw_destroy_plan(m_planA);
00031 fftw_destroy_plan(m_planB);
00032 fftw_destroy_plan(m_planX);
00033
00034 fftw_free(m_signalA);
00035 fftw_free(m_signalB);
00036 fftw_free(m_outShifted);
00037 fftw_free(m_out);
00038 fftw_free(m_outA);
00039 fftw_free(m_outB);
00040
00041 fftw_cleanup();
00042 }
00043
00044 void Correlation::correlateWithShift (int N,
00045 const double function1 [],
00046 const double function2 [],
00047 int &binStartMax,
00048 double &corrMax,
00049 double correlations []) const
00050 {
00051
00052
00053 int i;
00054
00055 ENGAUGE_ASSERT (N == m_N);
00056
00057
00058
00059
00060 double sumMean1 = 0, sumMean2 = 0, max1 = 0, max2 = 0;
00061 for (i = 0; i < N; i++) {
00062
00063 sumMean1 += function1 [i];
00064 sumMean2 += function2 [i];
00065 max1 = qMax (max1, function1 [i]);
00066 max2 = qMax (max2, function2 [i]);
00067
00068 }
00069
00070 double additiveNormalization1 = sumMean1 / N;
00071 double additiveNormalization2 = sumMean2 / N;
00072 double multiplicativeNormalization1 = 1.0 / max1;
00073 double multiplicativeNormalization2 = 1.0 / max2;
00074
00075
00076
00077 for (i = 0; i < N - 1; i++) {
00078
00079 m_signalA [i] [0] = 0.0;
00080 m_signalA [i] [1] = 0.0;
00081 m_signalB [i + N] [0] = 0.0;
00082 m_signalB [i + N] [1] = 0.0;
00083
00084 }
00085 for (i = 0; i < N; i++) {
00086
00087 m_signalA [i + N - 1] [0] = (function1 [i] - additiveNormalization1) * multiplicativeNormalization1;
00088 m_signalA [i + N - 1] [1] = 0.0;
00089 m_signalB [i] [0] = (function2 [i] - additiveNormalization2) * multiplicativeNormalization2;
00090 m_signalB [i] [1] = 0.0;
00091
00092 }
00093
00094 fftw_execute(m_planA);
00095 fftw_execute(m_planB);
00096
00097
00098 fftw_complex scale = {1.0/(2.0 * N - 1.0), 0.0};
00099 for (i = 0; i < 2 * N - 1; i++) {
00100
00101 fftw_complex term1 = {m_outA [i] [0], m_outA [i] [1]};
00102 fftw_complex term2 = {m_outB [i] [0], m_outB [i] [1] * -1.0};
00103 fftw_complex term3 = {scale [0], scale [1]};
00104 fftw_complex terms12 = {term1 [0] * term2 [0] - term1 [1] * term2 [1],
00105 term1 [0] * term2 [1] + term1 [1] * term2 [0]};
00106 m_out [i] [0] = terms12 [0] * term3 [0] - terms12 [1] * term3 [1];
00107 m_out [i] [1] = terms12 [0] * term3 [1] + terms12 [1] * term3 [0];
00108 }
00109
00110 fftw_execute(m_planX);
00111
00112
00113
00114 corrMax = 0.0;
00115 for (int i0AtLeft = 0; i0AtLeft < N; i0AtLeft++) {
00116
00117 int i0AtCenter = (i0AtLeft + N) % (2 * N - 1);
00118 fftw_complex shifted = {m_outShifted [i0AtCenter] [0], m_outShifted [i0AtCenter] [1]};
00119 double corr = qSqrt (shifted [0] * shifted [0] + shifted [1] * shifted [1]);
00120
00121 if ((i0AtLeft == 0) || (corr > corrMax)) {
00122 binStartMax = i0AtLeft;
00123 corrMax = corr;
00124 }
00125
00126
00127 correlations [i0AtLeft] = corr;
00128 }
00129 }
00130
00131 void Correlation::correlateWithoutShift (int N,
00132 const double function1 [],
00133 const double function2 [],
00134 double &corrMax) const
00135 {
00136
00137
00138 corrMax = 0.0;
00139
00140 for (int i = 0; i < N; i++) {
00141 corrMax += function1 [i] * function2 [i];
00142 }
00143 }