1 /*
2 SPDX-FileCopyrightText: 2012 Simon A. Eugster (Granjow) <simon.eu@gmail.com>
3 SPDX-License-Identifier: GPL-3.0-only OR LicenseRef-KDE-Accepted-GPL
4 */
5
6 #include "fftCorrelation.h"
7 #include <QElapsedTimer>
8 extern "C" {
9 #include "../external/kiss_fft/tools/kiss_fftr.h"
10 }
11
12 #include "kdenlive_debug.h"
13 #include <algorithm>
14 #include <vector>
15
correlate(const qint64 * left,const size_t leftSize,const qint64 * right,const size_t rightSize,qint64 * out_correlated)16 void FFTCorrelation::correlate(const qint64 *left, const size_t leftSize, const qint64 *right, const size_t rightSize, qint64 *out_correlated)
17 {
18 auto *correlatedFloat = new float[leftSize + rightSize + 1];
19 correlate(left, leftSize, right, rightSize, correlatedFloat);
20
21 // The correlation vector will have entries up to N (number of entries
22 // of the vector), so converting to integers will not lose that much
23 // of precision.
24 for (size_t i = 0; i < leftSize + rightSize + 1; ++i) {
25 out_correlated[i] = qint64(correlatedFloat[i]);
26 }
27 delete[] correlatedFloat;
28 }
29
correlate(const qint64 * left,const size_t leftSize,const qint64 * right,const size_t rightSize,float * out_correlated)30 void FFTCorrelation::correlate(const qint64 *left, const size_t leftSize, const qint64 *right, const size_t rightSize, float *out_correlated)
31 {
32 QElapsedTimer t;
33 t.start();
34
35 auto *leftF = new float[leftSize];
36 auto *rightF = new float[rightSize];
37
38 // First the qint64 values need to be normalized to floats
39 // Dividing by the max value is maybe not the best solution, but the
40 // maximum value after correlation should not be larger than the longest
41 // vector since each value should be at most 1
42 qint64 maxLeft = 1;
43 qint64 maxRight = 1;
44 for (size_t i = 0; i < leftSize; ++i) {
45 if (qAbs(left[i]) > maxLeft) {
46 maxLeft = qAbs(left[i]);
47 }
48 }
49 for (size_t i = 0; i < rightSize; ++i) {
50 if (qAbs(right[i]) > maxRight) {
51 maxRight = qAbs(right[i]);
52 }
53 }
54
55 // One side needs to be reversed, since multiplication in frequency domain (fourier space)
56 // calculates the convolution: \sum l[x]r[N-x] and not the correlation: \sum l[x]r[x]
57 for (size_t i = 0; i < leftSize; ++i) {
58 leftF[i] = float(left[i]) / maxLeft;
59 }
60 for (size_t i = 0; i < rightSize; ++i) {
61 rightF[rightSize - 1 - i] = float(right[i]) / maxRight;
62 }
63
64 // Now we can convolve to get the correlation
65 convolve(leftF, leftSize, rightF, rightSize, out_correlated);
66
67 qCDebug(KDENLIVE_LOG) << "Correlation (FFT based) computed in " << t.elapsed() << " ms.";
68 delete[] leftF;
69 delete[] rightF;
70 }
71
convolve(const float * left,const size_t leftSize,const float * right,const size_t rightSize,float * out_convolved)72 void FFTCorrelation::convolve(const float *left, const size_t leftSize, const float *right, const size_t rightSize, float *out_convolved)
73 {
74 QElapsedTimer time;
75 time.start();
76
77 // To avoid issues with repetition (we are dealing with cosine waves
78 // in the fourier domain) we need to pad the vectors to at least twice their size,
79 // otherwise convolution would convolve with the repeated pattern as well
80 size_t largestSize = std::max(leftSize, rightSize);
81
82 // The vectors must have the same size (same frequency resolution!) and should
83 // be a power of 2 (for FFT).
84 size_t size = 64;
85 while (size / 2 < largestSize) {
86 size = size << 1;
87 }
88
89 const size_t fft_size = size / 2 + 1;
90 kiss_fftr_cfg fftConfig = kiss_fftr_alloc(int(size), 0, nullptr, nullptr);
91 kiss_fftr_cfg ifftConfig = kiss_fftr_alloc(int(size), 1, nullptr, nullptr);
92 std::vector<kiss_fft_cpx> leftFFT(fft_size);
93 std::vector<kiss_fft_cpx> rightFFT(fft_size);
94 std::vector<kiss_fft_cpx> correlatedFFT(fft_size);
95
96 // Fill in the data into our new vectors with padding
97 std::vector<float> leftData(size, 0);
98 std::vector<float> rightData(size, 0);
99 std::vector<float> convolved(size);
100
101 std::copy(left, left + leftSize, leftData.begin());
102 std::copy(right, right + rightSize, rightData.begin());
103
104 // Fourier transformation of the vectors
105 kiss_fftr(fftConfig, &leftData[0], &leftFFT[0]);
106 kiss_fftr(fftConfig, &rightData[0], &rightFFT[0]);
107
108 // Convolution in spacial domain is a multiplication in fourier domain. O(n).
109 for (size_t i = 0; i < correlatedFFT.size(); ++i) {
110 correlatedFFT[i].r = leftFFT[i].r * rightFFT[i].r - leftFFT[i].i * rightFFT[i].i;
111 correlatedFFT[i].i = leftFFT[i].r * rightFFT[i].i + leftFFT[i].i * rightFFT[i].r;
112 }
113
114 // Inverse fourier transformation to get the convolved data.
115 // Insert one element at the beginning to obtain the same result
116 // that we also get with the nested for loop correlation.
117 *out_convolved = 0;
118 size_t out_size = leftSize + rightSize + 1;
119
120 kiss_fftri(ifftConfig, &correlatedFFT[0], &convolved[0]);
121 std::copy(convolved.begin(), convolved.begin() + int(out_size) - 1, out_convolved + 1);
122
123 // Finally some cleanup.
124 kiss_fftr_free(fftConfig);
125 kiss_fftr_free(ifftConfig);
126
127 qCDebug(KDENLIVE_LOG) << "FFT convolution computed. Time taken: " << time.elapsed() << " ms";
128 }
129