FFT based correlation works.

[kdenlive] / src / lib / audio / fftCorrelation.cpp

diff --git a/src/lib/audio/fftCorrelation.cpp b/src/lib/audio/fftCorrelation.cpp
index 229d4ec48e4072aea8cbd09c121e911d29bd8f26..fa31409665ecb7e28ecd6f209ce710fa66c85c8a 100644 (file)
--- a/src/lib/audio/fftCorrelation.cpp
+++ b/src/lib/audio/fftCorrelation.cpp
@@ -1,3 +1,13 @@
+/*
+Copyright (C) 2012  Simon A. Eugster (Granjow)  <simon.eu@gmail.com>
+This file is part of kdenlive. See www.kdenlive.org.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+*/
+
 #include "fftCorrelation.h"
 
 extern "C"
@@ -6,35 +16,86 @@ extern "C"
 }
 
 #include <QTime>
-#include <QDebug>
+#include <iostream>
 #include <algorithm>
 
-void FFTCorrelation::correlate(float *left, int leftSize, float *right, int rightSize,
-                          float **out_correlationData, int &out_size)
+void FFTCorrelation::correlate(const int64_t *left, const int leftSize,
+                               const int64_t *right, const int rightSize,
+                               float **out_correlated, int &out_size)
+{
+    QTime t;
+    t.start();
+
+    float leftF[leftSize];
+    float rightF[rightSize];
+
+    // First the int64_t values need to be normalized to floats
+    // Dividing by the max value is maybe not the best solution, but the
+    // maximum value after correlation should not be larger than the longest
+    // vector since each value should be at most 1
+    int64_t maxLeft = 0;
+    int64_t maxRight = 0;
+    for (int i = 0; i < leftSize; i++) {
+        if (labs(left[i]) > maxLeft) {
+            maxLeft = labs(left[i]);
+        }
+    }
+    for (int i = 0; i < rightSize; i++) {
+        if (labs(right[i]) > maxRight) {
+            maxRight = labs(right[i]);
+        }
+    }
+
+
+    // One side needs to be reverted, since multiplication in frequency domain (fourier space)
+    // calculates the convolution: \sum l[x]r[N-x] and not the correlation: \sum l[x]r[x]
+    for (int i = 0; i < leftSize; i++) {
+        leftF[leftSize-1 - i] = double(left[i])/maxLeft;
+    }
+    for (int i = 0; i < rightSize; i++) {
+        rightF[i] = double(right[i])/maxRight;
+    }
+
+    // Now we can convolve to get the correlation
+    convolute(leftF, leftSize, rightF, rightSize, out_correlated, out_size);
+
+    std::cout << "Correlation (FFT based) computed in " << t.elapsed() << " ms." << std::endl;
+}
+
+void FFTCorrelation::convolute(const float *left, const int leftSize,
+                               const float *right, const int rightSize,
+                               float **out_convolved, int &out_size)
 {
     QTime time;
     time.start();
 
+
+    // To avoid issues with repetition (we are dealing with cosine waves
+    // in the fourier domain) we need to pad the vectors to at least twice their size,
+    // otherwise convolution would convolve with the repeated pattern as well
     int largestSize = leftSize;
     if (rightSize > largestSize) {
         largestSize = rightSize;
     }
 
+    // The vectors must have the same size (same frequency resolution!) and should
+    // be a power of 2 (for FFT).
     int size = 64;
     while (size/2 < largestSize) {
         size = size << 1;
     }
 
-    kiss_fftr_cfg fftConfig = kiss_fftr_alloc(size, false,NULL,NULL);
+    kiss_fftr_cfg fftConfig = kiss_fftr_alloc(size, false, NULL,NULL);
     kiss_fftr_cfg ifftConfig = kiss_fftr_alloc(size, true, NULL,NULL);
     kiss_fft_cpx leftFFT[size/2];
     kiss_fft_cpx rightFFT[size/2];
     kiss_fft_cpx correlatedFFT[size/2];
 
 
+    // Fill in the data into our new vectors with padding
     float leftData[size];
     float rightData[size];
-    *out_correlationData = new float[size];
+    *out_convolved = new float[size];
 
     std::fill(leftData, leftData+size, 0);
     std::fill(rightData, rightData+size, 0);
@@ -42,16 +103,23 @@ void FFTCorrelation::correlate(float *left, int leftSize, float *right, int righ
     std::copy(left, left+leftSize, leftData);
     std::copy(right, right+rightSize, rightData);
 
+    // Fourier transformation of the vectors
     kiss_fftr(fftConfig, leftData, leftFFT);
     kiss_fftr(fftConfig, rightData, rightFFT);
 
+    // Convolution in spacial domain is a multiplication in fourier domain. O(n).
     for (int i = 0; i < size/2; i++) {
         correlatedFFT[i].r = leftFFT[i].r*rightFFT[i].r - leftFFT[i].i*rightFFT[i].i;
         correlatedFFT[i].i = leftFFT[i].r*rightFFT[i].i + leftFFT[i].i*rightFFT[i].r;
     }
 
-    kiss_fftri(ifftConfig, correlatedFFT, *out_correlationData);
+    // Inverse fourier tranformation to get the convolved data
+    kiss_fftri(ifftConfig, correlatedFFT, *out_convolved);
     out_size = size;
 
-    qDebug() << "FFT correlation computed. Time taken: " << time.elapsed() << " ms";
+    // Finally some cleanup.
+    kiss_fftr_free(fftConfig);
+    kiss_fftr_free(ifftConfig);
+
+    std::cout << "FFT convolution computed. Time taken: " << time.elapsed() << " ms" << std::endl;
 }