raDIYo/swfftcalc.cpp

/***************************************************************************

    source::worx raDIYo
    Copyright © 2020-2022 c.holzheuer
    chris@sourceworx.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 2 of the License, or
    (at your option) any later version.

***************************************************************************/


#include <QtDebug>
//#include <QtAlgorithms>
//#include <qalgorithms.h>


#include <qmath.h>
#include <swfftcalc.h>


//const qint16  PCMS16MaxValue     =  32767;
const quint16 PCMS16MaxAmplitude =  32768; // because minimum is -32768


qreal xpcmToReal(qint16 pcm)
{
    return qreal(pcm) / PCMS16MaxAmplitude;
}


SWFFTCalc::SWFFTCalc( QObject* parent )
    : QObject( parent )
{

    _timer.start();
    // window functions are used to filter some undesired
    // information for fft calculation.
    _hannWindow.resize( SWNumSamples );

    // the complex frame that is sent to fft function
    _complexFrame.resize( SWNumSamples );

    // only half spectrum is used because of the simetry property
    _spectrum.resize( SWNumSamples/2 );

    // logscale is used for audio spectrum display
    //_logScale.resize( SWNumSamples/2+1 );

    // by default, spectrum is log scaled (compressed)
    //compressed = false;

    // window function (HANN)
    for( int i=0; i<SWNumSamples; i++ )
        _hannWindow[i] = 0.5 * ( 1 - cos( (2*PI*i) /  (SWNumSamples-1) ) );

    //double hannWindow = 0.5 * (1 - cos((2 * M_PI * i) / (n - 1)));
    // x = 0.5 * (1 - qCos((2 * M_PI * i) / (SWNumSamples - 1)));
    //  case HannWindow:
    //      x = 0.5 * (1 - qCos((2 * M_PI * i) / (SWNumSamples - 1)));

    // the log scale
    //for( int i=0; i<=SWNumSamples/2; i++ )
    //  _logScale[i] = powf( SWNumSamples/2, (float) 2*i / SWNumSamples ) - 0.5f;


}


SWFFTCalc::~SWFFTCalc()
{

}



/*

    Also: Wir bekommen zu willkürlichen Zeitpunkten im millisekunden
    Bereich einen Buffer mit Rohdaten im Format X.

    TODO:
        - Format X PeakValue rausfinden
        - einen Buffer<double> anlegen
        - mit Peak(X) auf [0.0 ... 1.0 [ normalisieren
        - FFT drüber

*/


/**
 * https://stackoverflow.com/questions/46947668/draw-waveform-from-raw-data-using-qaudioprobe
 * @brief Track::getPeakValue
 * @param format
 * @return The peak value
 */


void SWFFTCalc::onAudioProbed( const QAudioBuffer &audiobuffer )
{


    //
    // Step 1. Daten anreichern bis ca. 100 ms
    //

    QAudioFormat format = audiobuffer.format();

    if( !format.isValid()
      || format.channelCount() != 2
      || format.sampleType() == QAudioFormat::Unknown )
        return;

    //int duration = format.durationForFrames( audiobuffer.frameCount() ) / 1000;
    // _timer.start();

    //if( duration < 70 )
        //return;

    // Initialize data array
    _complexFrame = CArray( 0.0, SWNumSamples );
    const char* ptr = (char*) audiobuffer.constData();
    int to = qMin( SWNumSamples, audiobuffer.frameCount() );

    //for( int i=0; i<SWNumSamples; ++i )
    for( int i=0; i<to; ++i )
    {
        const qint16 pcmSample = *reinterpret_cast<const qint16*>( ptr );
        // Scale down to range [-1.0, 1.0]
        //_input[i] = pcmToReal(pcmSample) * _hannWindow[i];
        _complexFrame[i] = Complex( _hannWindow[i] *  xpcmToReal(pcmSample), 0 );
        ptr += 4; //bytesPerSample;
    }

    // do the magic
    bareFFT( _complexFrame );

    // Analyze output to obtain amplitude and phase for each frequency
    //for( int i=2; i < SWNumSamples/2; ++i )
    for( int i=0; i < SWNumSamples/2; ++i )
    {

        const qreal magnitude = abs( _complexFrame[i] );

        //??
        //qreal amplitude = qLn( magnitude );
        //qreal amplitude = 0.27 * qLn( magnitude );
        qreal amplitude = SWFudgeFactor * qLn( magnitude );
        //qreal amplitude = SpectrumAnalyserMultiplier * qLn( magnitude );
        //qreal amplitude = SpectrumAnalyserMultiplier * magnitude;
        //qreal amplitude = SWFudgeFactor * magnitude;

        //qreal amplitude = magnitude;

        _spectrum[i] = qBound( 0.0, amplitude, 1.0 );

    }

    emit spectrumReady( _spectrum );


}

void SWFFTCalc::bareFFT( CArray& x )
{

    // fft function.
    //
    // ATTENTION ///////////
    //
    // size of x must be a power of two
    //
    ////////////////////

    const size_t N = x.size();
    if (N <= 1)
        return;

    // divide
    CArray even = x[ std::slice( 0, N/2, 2 ) ];
    CArray  odd = x[ std::slice( 1, N/2, 2 ) ];

    // conquer
    bareFFT( even );
    bareFFT( odd );

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar( 1.0, -2 * PI * k / N ) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }

}