]> git.sesse.net Git - kdenlive/blob - src/lib/audio/fftCorrelation.cpp
FFT based correlation works.
[kdenlive] / 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                                float **out_correlated, int &out_size)
25 {
26     QTime t;
27     t.start();
28
29     float leftF[leftSize];
30     float rightF[rightSize];
31
32     // First the int64_t values need to be normalized to floats
33     // Dividing by the max value is maybe not the best solution, but the
34     // maximum value after correlation should not be larger than the longest
35     // vector since each value should be at most 1
36     int64_t maxLeft = 0;
37     int64_t maxRight = 0;
38     for (int i = 0; i < leftSize; i++) {
39         if (labs(left[i]) > maxLeft) {
40             maxLeft = labs(left[i]);
41         }
42     }
43     for (int i = 0; i < rightSize; i++) {
44         if (labs(right[i]) > maxRight) {
45             maxRight = labs(right[i]);
46         }
47     }
48
49
50     // One side needs to be reverted, since multiplication in frequency domain (fourier space)
51     // calculates the convolution: \sum l[x]r[N-x] and not the correlation: \sum l[x]r[x]
52     for (int i = 0; i < leftSize; i++) {
53         leftF[leftSize-1 - i] = double(left[i])/maxLeft;
54     }
55     for (int i = 0; i < rightSize; i++) {
56         rightF[i] = double(right[i])/maxRight;
57     }
58
59     // Now we can convolve to get the correlation
60     convolute(leftF, leftSize, rightF, rightSize, out_correlated, out_size);
61
62     std::cout << "Correlation (FFT based) computed in " << t.elapsed() << " ms." << std::endl;
63 }
64
65 void FFTCorrelation::convolute(const float *left, const int leftSize,
66                                const float *right, const int rightSize,
67                                float **out_convolved, int &out_size)
68 {
69     QTime time;
70     time.start();
71
72
73     // To avoid issues with repetition (we are dealing with cosine waves
74     // in the fourier domain) we need to pad the vectors to at least twice their size,
75     // otherwise convolution would convolve with the repeated pattern as well
76     int largestSize = leftSize;
77     if (rightSize > largestSize) {
78         largestSize = rightSize;
79     }
80
81     // The vectors must have the same size (same frequency resolution!) and should
82     // be a power of 2 (for FFT).
83     int size = 64;
84     while (size/2 < largestSize) {
85         size = size << 1;
86     }
87
88     kiss_fftr_cfg fftConfig = kiss_fftr_alloc(size, false, NULL,NULL);
89     kiss_fftr_cfg ifftConfig = kiss_fftr_alloc(size, true, NULL,NULL);
90     kiss_fft_cpx leftFFT[size/2];
91     kiss_fft_cpx rightFFT[size/2];
92     kiss_fft_cpx correlatedFFT[size/2];
93
94
95     // Fill in the data into our new vectors with padding
96     float leftData[size];
97     float rightData[size];
98     *out_convolved = new float[size];
99
100     std::fill(leftData, leftData+size, 0);
101     std::fill(rightData, rightData+size, 0);
102
103     std::copy(left, left+leftSize, leftData);
104     std::copy(right, right+rightSize, rightData);
105
106     // Fourier transformation of the vectors
107     kiss_fftr(fftConfig, leftData, leftFFT);
108     kiss_fftr(fftConfig, rightData, rightFFT);
109
110     // Convolution in spacial domain is a multiplication in fourier domain. O(n).
111     for (int i = 0; i < size/2; i++) {
112         correlatedFFT[i].r = leftFFT[i].r*rightFFT[i].r - leftFFT[i].i*rightFFT[i].i;
113         correlatedFFT[i].i = leftFFT[i].r*rightFFT[i].i + leftFFT[i].i*rightFFT[i].r;
114     }
115
116     // Inverse fourier tranformation to get the convolved data
117     kiss_fftri(ifftConfig, correlatedFFT, *out_convolved);
118     out_size = size;
119
120     // Finally some cleanup.
121     kiss_fftr_free(fftConfig);
122     kiss_fftr_free(ifftConfig);
123
124     std::cout << "FFT convolution computed. Time taken: " << time.elapsed() << " ms" << std::endl;
125 }