git.sesse.net Git - kdenlive/blob - src/lib/audio/fftCorrelation.cpp

   1  /*
   2 Copyright (C) 2012  Simon A. Eugster (Granjow)  <simon.eu@gmail.com>
   3 This file is part of kdenlive. See www.kdenlive.org.
   4
   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.
   9 */
  10
  11 #include "fftCorrelation.h"
  12
  13 extern "C"
  14 {
  15 #include "../external/kiss_fft/tools/kiss_fftr.h"
  16 }
  17
  18 #include <QTime>
  19 #include <iostream>
  20 #include <algorithm>
  21
  22 void FFTCorrelation::correlate(const int64_t *left, const int leftSize,
  23                                const int64_t *right, const int rightSize,
  24                                int64_t *out_correlated)
  25 {
  26     float correlatedFloat[leftSize+rightSize+1];
  27     correlate(left, leftSize, right, rightSize, correlatedFloat);
  28
  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
  31     // of precision.
  32     for (int i = 0; i < leftSize+rightSize+1; ++i) {
  33         out_correlated[i] = correlatedFloat[i];
  34     }
  35 }
  36
  37 void FFTCorrelation::correlate(const int64_t *left, const int leftSize,
  38                                const int64_t *right, const int rightSize,
  39                                float *out_correlated)
  40 {
  41     QTime t;
  42     t.start();
  43
  44     float leftF[leftSize];
  45     float rightF[rightSize];
  46
  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
  51     int64_t maxLeft = 1;
  52     int64_t maxRight = 1;
  53     for (int i = 0; i < leftSize; ++i) {
  54         if (labs(left[i]) > maxLeft) {
  55             maxLeft = labs(left[i]);
  56         }
  57     }
  58     for (int i = 0; i < rightSize; ++i) {
  59         if (labs(right[i]) > maxRight) {
  60             maxRight = labs(right[i]);
  61         }
  62     }
  63
  64
  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;
  69     }
  70     for (int i = 0; i < rightSize; ++i) {
  71         rightF[rightSize-1 - i] = double(right[i])/maxRight;
  72     }
  73
  74     // Now we can convolve to get the correlation
  75     convolve(leftF, leftSize, rightF, rightSize, out_correlated);
  76
  77     std::cout << "Correlation (FFT based) computed in " << t.elapsed() << " ms." << std::endl;
  78 }
  79
  80 void FFTCorrelation::convolve(const float *left, const int leftSize,
  81                                const float *right, const int rightSize,
  82                                float *out_convolved)
  83 {
  84     QTime time;
  85     time.start();
  86
  87
  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;
  94     }
  95
  96     // The vectors must have the same size (same frequency resolution!) and should
  97     // be a power of 2 (for FFT).
  98     int size = 64;
  99     while (size/2 < largestSize) {
 100         size = size << 1;
 101     }
 102
 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];
 108
 109
 110     // Fill in the data into our new vectors with padding
 111     float leftData[size];
 112     float rightData[size];
 113     float convolved[size];
 114
 115     std::fill(leftData, leftData+size, 0);
 116     std::fill(rightData, rightData+size, 0);
 117
 118     std::copy(left, left+leftSize, leftData);
 119     std::copy(right, right+rightSize, rightData);
 120
 121     // Fourier transformation of the vectors
 122     kiss_fftr(fftConfig, leftData, leftFFT);
 123     kiss_fftr(fftConfig, rightData, rightFFT);
 124
 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;
 129     }
 130
 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.
 134     *out_convolved = 0;
 135     int out_size = leftSize+rightSize+1;
 136
 137     kiss_fftri(ifftConfig, correlatedFFT, convolved);
 138     std::copy(convolved, convolved+out_size-1, out_convolved+1);
 139
 140     // Finally some cleanup.
 141     kiss_fftr_free(fftConfig);
 142     kiss_fftr_free(ifftConfig);
 143
 144     std::cout << "FFT convolution computed. Time taken: " << time.elapsed() << " ms" << std::endl;
 145 }