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

   1 /***************************************************************************
   2  *   Copyright (C) 2010 by Simon Andreas Eugster (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 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  ***************************************************************************/
  10
  11 #include <math.h>
  12 #include <iostream>
  13
  14 #include <QString>
  15
  16 #include "fftTools.h"
  17
  18 // Uncomment for debugging, like writing a GNU Octave .m file to /tmp
  19 //#define DEBUG_FFTTOOLS
  20
  21 #ifdef DEBUG_FFTTOOLS
  22 #include <QDebug>
  23 #include <QTime>
  24 #include <fstream>
  25 #endif
  26
  27 FFTTools::FFTTools() :
  28         m_fftCfgs(),
  29         m_windowFunctions()
  30 {
  31 }
  32 FFTTools::~FFTTools()
  33 {
  34     QHash<QString, kiss_fftr_cfg>::iterator i;
  35     for (i = m_fftCfgs.begin(); i != m_fftCfgs.end(); ++i) {
  36         free(*i);
  37     }
  38 }
  39
  40 const QString FFTTools::windowSignature(const WindowType windowType, const int size, const float param)
  41 {
  42     return QString("s%1_t%2_p%3").arg(size).arg(windowType).arg(param, 0, 'f', 3);
  43 }
  44 const QString FFTTools::cfgSignature(const int size)
  45 {
  46     return QString("s%1").arg(size);
  47 }
  48
  49 // http://cplusplus.syntaxerrors.info/index.php?title=Cannot_declare_member_function_%E2%80%98static_int_Foo::bar%28%29%E2%80%99_to_have_static_linkage
  50 const QVector<float> FFTTools::window(const WindowType windowType, const int size, const float param)
  51 {
  52     Q_ASSERT(size > 0);
  53
  54     // Deliberately avoid converting size to a float
  55     // to keep mid an integer.
  56     float mid = (size-1)/2;
  57     float max = size-1;
  58     QVector<float> window;
  59
  60     switch (windowType) {
  61     case Window_Rect:
  62         return QVector<float>(size+1, 1);
  63         break;
  64     case Window_Triangle:
  65         window = QVector<float>(size+1);
  66
  67         for (int x = 0; x < mid; x++) {
  68             window[x] = x/mid + (mid-x)/mid*param;
  69         }
  70         for (int x = mid; x < size; x++) {
  71             window[x] = (x-mid)/(max-mid) * param + (max-x)/(max-mid);
  72         }
  73         window[size] = .5 + param/2;
  74
  75 #ifdef DEBUG_FFTTOOLS
  76         qDebug() << "Triangle window (factor " << window[size] << "):";
  77         for (int i = 0; i < size; ++i) {
  78             qDebug() << window[i];
  79         }
  80         qDebug() << "Triangle window end.";
  81 #endif
  82
  83         return window;
  84         break;
  85     case Window_Hamming:
  86         // Use a quick version of the Hamming window here: Instead of
  87         // interpolating values between (-max/2) and (max/2)
  88         // we use integer values instead, ranging from -mid to (max-mid).
  89         window = QVector<float>(size+1);
  90
  91         for (int x = 0; x < size; x++) {
  92             window[x] = .54 + .46 * cos( 2*M_PI*(x-mid) / size );
  93         }
  94
  95         // Integrating the cosine over the window function results in
  96         // an area of 0; So only the constant factor 0.54 counts.
  97         window[size] = .54;
  98
  99 #ifdef DEBUG_FFTTOOLS
 100         qDebug() << "Hanning window (factor " << window[size] << "):";
 101         for (int i = 0; i < size; ++i) {
 102             qDebug() << window[i];
 103         }
 104         qDebug() << "Hanning window end.";
 105 #endif
 106
 107         return window;
 108         break;
 109     }
 110     Q_ASSERT(false);
 111     return QVector<float>();
 112 }
 113
 114 void FFTTools::fftNormalized(const QVector<int16_t> audioFrame, const uint channel, const uint numChannels, float *freqSpectrum,
 115                              const WindowType windowType, const uint windowSize, const float param)
 116 {
 117 #ifdef DEBUG_FFTTOOLS
 118     QTime start = QTime::currentTime();
 119 #endif
 120
 121     const uint numSamples = audioFrame.size()/numChannels;
 122
 123     Q_ASSERT((windowSize & 1) == 0);
 124     Q_ASSERT(windowSize > 0);
 125
 126     const QString cfgSig = cfgSignature(windowSize);
 127     const QString winSig = windowSignature(windowType, windowSize, param);
 128
 129
 130     // Get the kiss_fft configuration from the config cache
 131     // or build a new configuration if the requested one is not available.
 132     kiss_fftr_cfg myCfg;
 133     if (m_fftCfgs.contains(cfgSig)) {
 134 #ifdef DEBUG_FFTTOOLS
 135         qDebug() << "Re-using FFT configuration with size " << windowSize;
 136 #endif
 137         myCfg = m_fftCfgs.value(cfgSig);
 138     } else {
 139 #ifdef DEBUG_FFTTOOLS
 140         qDebug() << "Creating FFT configuration with size " << windowSize;
 141 #endif
 142         myCfg = kiss_fftr_alloc(windowSize, false,NULL,NULL);
 143         m_fftCfgs.insert(cfgSig, myCfg);
 144     }
 145
 146     // Get the window function from the cache
 147     // (except for a rectangular window; nothing to to there.
 148     QVector<float> window;
 149     float windowScaleFactor = 1;
 150     if (windowType != FFTTools::Window_Rect) {
 151
 152         if (m_windowFunctions.contains(winSig)) {
 153 #ifdef DEBUG_FFTTOOLS
 154             qDebug() << "Re-using window function with signature " << winSig;
 155 #endif
 156             window = m_windowFunctions.value(winSig);
 157         } else {
 158 #ifdef DEBUG_FFTTOOLS
 159             qDebug() << "Building new window function with signature " << winSig;
 160 #endif
 161             window = FFTTools::window(windowType, windowSize, 0);
 162             m_windowFunctions.insert(winSig, window);
 163         }
 164         windowScaleFactor = 1.0/window[windowSize];
 165     }
 166
 167
 168     // Prepare frequency space vector. The resulting FFT vector is only half as long.
 169     kiss_fft_cpx freqData[windowSize/2];
 170     float data[windowSize];
 171
 172     // Copy the first channel's audio into a vector for the FFT display;
 173     // Fill the data vector indices that cannot be covered with sample data with 0
 174     if (numSamples < windowSize) {
 175         std::fill(&data[numSamples], &data[windowSize-1], 0);
 176     }
 177     // Normalize signals to [0,1] to get correct dB values later on
 178     for (uint i = 0; i < numSamples && i < windowSize; ++i) {
 179         // Performance note: Benchmarking has shown that using the if/else inside the loop
 180         // does not do noticeable worse than keeping it outside (perhaps the branch predictor
 181         // is good enough), so it remains in there for better readability.
 182         if (windowType != FFTTools::Window_Rect) {
 183             data[i] = (float) audioFrame.data()[i*numChannels + channel] / 32767.0f * window[i];
 184         } else {
 185             data[i] = (float) audioFrame.data()[i*numChannels + channel] / 32767.0f;
 186         }
 187     }
 188
 189     // Calculate the Fast Fourier Transform for the input data
 190     kiss_fftr(myCfg, data, freqData);
 191
 192
 193     // Logarithmic scale: 20 * log ( 2 * magnitude / N ) with magnitude = sqrt(r² + i²)
 194     // with N = FFT size (after FFT, 1/2 window size)
 195     for (uint i = 0; i < windowSize/2; ++i) {
 196         // Logarithmic scale: 20 * log ( 2 * magnitude / N ) with magnitude = sqrt(r² + i²)
 197         // with N = FFT size (after FFT, 1/2 window size)
 198         freqSpectrum[i] = 20*log(pow(pow(fabs(freqData[i].r * windowScaleFactor),2) + pow(fabs(freqData[i].i * windowScaleFactor),2), .5)/((float)windowSize/2.0f))/log(10);;
 199     }
 200
 201
 202 #ifdef DEBUG_FFTTOOLS
 203     std::ofstream mFile;
 204     mFile.open("/tmp/freq.m");
 205     if (!mFile) {
 206         qDebug() << "Opening file failed.";
 207     } else {
 208         mFile << "val = [ ";
 209
 210         for (int sample = 0; sample < 256; sample++) {
 211             mFile << data[sample] << " ";
 212         }
 213         mFile << " ];\n";
 214
 215         mFile << "freq = [ ";
 216         for (int sample = 0; sample < 256; sample++) {
 217             mFile << freqData[sample].r << "+" << freqData[sample].i << "*i ";
 218         }
 219         mFile << " ];\n";
 220
 221         mFile.close();
 222         qDebug() << "File written.";
 223     }
 224 #endif
 225
 226 #ifdef DEBUG_FFTTOOLS
 227     qDebug() << "Calculated FFT in " << start.elapsed() << " ms.";
 228 #endif
 229 }
 230
 231 const QVector<float> FFTTools::interpolatePeakPreserving(const QVector<float> in, const uint targetSize, uint left, uint right, float fill)
 232 {
 233 #ifdef DEBUG_FFTTOOLS
 234     QTime start = QTime::currentTime();
 235 #endif
 236
 237     if (right == 0) {
 238         right = in.size()-1;
 239     }
 240     Q_ASSERT(targetSize > 0);
 241     Q_ASSERT(left < right);
 242
 243     QVector<float> out(targetSize);
 244
 245
 246     float x;
 247     float x_prev = 0;
 248     uint xi;
 249     uint i;
 250     if (((float) (right-left))/targetSize < 2) {
 251         for (i = 0; i < targetSize; ++i) {
 252
 253             // i:  Target index
 254             // x:  Interpolated source index (float!)
 255             // xi: floor(x)
 256
 257             // Transform [0,targetSize-1] to [left,right]
 258             x = ((float) i) / (targetSize-1) * (right-left) + left;
 259             xi = (int) floor(x);
 260
 261             if (x > in.size()-1) {
 262                 // This may happen if right > in.size()-1; Fill the rest of the vector
 263                 // with the default value now.
 264                 break;
 265             }
 266
 267
 268             // Use linear interpolation in order to get smoother display
 269             if (xi == 0 || xi == (uint) in.size()-1) {
 270                 // ... except if we are at the left or right border of the input sigal.
 271                 // Special case here since we consider previous and future values as well for
 272                 // the actual interpolation (not possible here).
 273                 out[i] = in[xi];
 274             } else {
 275                 if (in[xi] > in[xi+1]
 276                     && x_prev < xi) {
 277                     // This is a hack to preserve peaks.
 278                     // Consider f = {0, 100, 0}
 279                     //          x = {0.5,  1.5}
 280                     // Then x is 50 both times, and the 100 peak is lost.
 281                     // Get it back here for the first x after the peak (which is at xi).
 282                     // (x is the first after the peak if the previous x was smaller than floor(x).)
 283                     out[i] = in[xi];
 284                 } else {
 285                     out[i] =   (xi+1 - x) * in[xi]
 286                           + (x - xi)   * in[xi+1];
 287                 }
 288             }
 289             x_prev = x;
 290         }
 291     } else {
 292         // If there are more than 2 samples per pixel in average, then use the maximum of them
 293         // since by only looking at the left sample we might miss some maxima.
 294         uint src = left;
 295 //         int xi_prev = 0;
 296         int points;
 297
 298 #ifdef DEBUG_FFTTOOLS
 299         qDebug() << "Interpolation: Ratio over 2; using maximum interpolation";
 300 #endif
 301
 302         for (i = 0; i < targetSize; ++i) {
 303
 304             // x:  right bound
 305             // xi: floor(x)
 306             x = ((float) (i+1)) / (targetSize-1) * (right-left) + left;
 307             xi = (int) floor(x);
 308             points = 0;
 309
 310             out[i] = fill;
 311
 312             for (; src < xi && src < (uint) in.size(); src++) {
 313                 if (out[i] < in[src]) {
 314                     out[i] = in[src];
 315                 }
 316                 points++;
 317             }
 318
 319 //             xi_prev = xi;
 320         }
 321     }
 322     // Fill the rest of the vector if the right border exceeds the input vector.
 323     for (; i < targetSize; ++i) {
 324         out[i] = fill;
 325     }
 326
 327 #ifdef DEBUG_FFTTOOLS
 328     qDebug() << "Interpolated " << targetSize << " nodes from " << in.size() << " input points in " << start.elapsed() << " ms";
 329 #endif
 330
 331     return out;
 332 }
 333
 334 #ifdef DEBUG_FFTTOOLS
 335 #undef DEBUG_FFTTOOLS
 336 #endif