2 Copyright (C) 2012 Simon A. Eugster (Granjow) <simon.eu@gmail.com>
3 This file is part of kdenlive. See www.kdenlive.org.
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
11 #include "fftCorrelation.h"
15 #include "../external/kiss_fft/tools/kiss_fftr.h"
22 void FFTCorrelation::correlate(const int64_t *left, const int leftSize,
23 const int64_t *right, const int rightSize,
24 int64_t *out_correlated)
26 float correlatedFloat[leftSize+rightSize+1];
27 correlate(left, leftSize, right, rightSize, correlatedFloat);
29 // The correlation vector will have entries up to N (number of entries
30 // of the vector), so converting to integers will not lose that much
32 for (int i = 0; i < leftSize+rightSize+1; ++i) {
33 out_correlated[i] = correlatedFloat[i];
37 void FFTCorrelation::correlate(const int64_t *left, const int leftSize,
38 const int64_t *right, const int rightSize,
39 float *out_correlated)
44 float leftF[leftSize];
45 float rightF[rightSize];
47 // First the int64_t values need to be normalized to floats
48 // Dividing by the max value is maybe not the best solution, but the
49 // maximum value after correlation should not be larger than the longest
50 // vector since each value should be at most 1
53 for (int i = 0; i < leftSize; ++i) {
54 if (labs(left[i]) > maxLeft) {
55 maxLeft = labs(left[i]);
58 for (int i = 0; i < rightSize; ++i) {
59 if (labs(right[i]) > maxRight) {
60 maxRight = labs(right[i]);
65 // One side needs to be reverted, since multiplication in frequency domain (fourier space)
66 // calculates the convolution: \sum l[x]r[N-x] and not the correlation: \sum l[x]r[x]
67 for (int i = 0; i < leftSize; ++i) {
68 leftF[i] = double(left[i])/maxLeft;
70 for (int i = 0; i < rightSize; ++i) {
71 rightF[rightSize-1 - i] = double(right[i])/maxRight;
74 // Now we can convolve to get the correlation
75 convolve(leftF, leftSize, rightF, rightSize, out_correlated);
77 std::cout << "Correlation (FFT based) computed in " << t.elapsed() << " ms." << std::endl;
80 void FFTCorrelation::convolve(const float *left, const int leftSize,
81 const float *right, const int rightSize,
88 // To avoid issues with repetition (we are dealing with cosine waves
89 // in the fourier domain) we need to pad the vectors to at least twice their size,
90 // otherwise convolution would convolve with the repeated pattern as well
91 int largestSize = leftSize;
92 if (rightSize > largestSize) {
93 largestSize = rightSize;
96 // The vectors must have the same size (same frequency resolution!) and should
97 // be a power of 2 (for FFT).
99 while (size/2 < largestSize) {
103 kiss_fftr_cfg fftConfig = kiss_fftr_alloc(size, false, NULL,NULL);
104 kiss_fftr_cfg ifftConfig = kiss_fftr_alloc(size, true, NULL,NULL);
105 kiss_fft_cpx leftFFT[size/2];
106 kiss_fft_cpx rightFFT[size/2];
107 kiss_fft_cpx correlatedFFT[size/2];
110 // Fill in the data into our new vectors with padding
111 float leftData[size];
112 float rightData[size];
113 float convolved[size];
115 std::fill(leftData, leftData+size, 0);
116 std::fill(rightData, rightData+size, 0);
118 std::copy(left, left+leftSize, leftData);
119 std::copy(right, right+rightSize, rightData);
121 // Fourier transformation of the vectors
122 kiss_fftr(fftConfig, leftData, leftFFT);
123 kiss_fftr(fftConfig, rightData, rightFFT);
125 // Convolution in spacial domain is a multiplication in fourier domain. O(n).
126 for (int i = 0; i < size/2; ++i) {
127 correlatedFFT[i].r = leftFFT[i].r*rightFFT[i].r - leftFFT[i].i*rightFFT[i].i;
128 correlatedFFT[i].i = leftFFT[i].r*rightFFT[i].i + leftFFT[i].i*rightFFT[i].r;
131 // Inverse fourier tranformation to get the convolved data.
132 // Insert one element at the beginning to obtain the same result
133 // that we also get with the nested for loop correlation.
135 int out_size = leftSize+rightSize+1;
137 kiss_fftri(ifftConfig, correlatedFFT, convolved);
138 std::copy(convolved, convolved+out_size-1, out_convolved+1);
140 // Finally some cleanup.
141 kiss_fftr_free(fftConfig);
142 kiss_fftr_free(ifftConfig);
144 std::cout << "FFT convolution computed. Time taken: " << time.elapsed() << " ms" << std::endl;