swfftcalc.cpp

/***************************************************************************
 
    source::worx raDIYo
    Copyright © 2020-2022 c.holzheuer
    c.holzheuer@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 pcmToReal(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;
 
 
 
    /*
 
    //sum result into 44 bins for visualization
    int steps = n/44;
    std::vector<double> bins;
    for (int i=0; i<44; i++) {
        double bin = 0;
        for (int q=0; q<steps; q++) {
            //calculate amplitude = sqrt(re*re+im*im)
            double amplitude = sqrt(y[i*steps+q][0]*y[i*steps+q][0]+y[i*steps+q][1]*y[i*steps+q][1]);
            //overwrite NaN with 0. Got a lot of them apparently, without these two lines it crashes immediatly, I guess somewhere downstream
            if (amplitude != amplitude) bin += 0;
            //convert to dB
            else bin += 20*log10(amplitude);
        }
        bins.push_back(bin);
        qDebug() << bins[i];
    }
    */
}
 
 
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
 
 
*/
 
 
void SWFFTCalc::collectFrames( 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;
    //qDebug() << "onAudioProbed:" << duration << " ms frames:" << audiobuffer.frameCount() << " time:" << _timer.elapsed();
    //qDebug() << "zeit: " << _timer.elapsed() << " ms";
    //_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] *  pcmToReal(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.2 * qLn( magnitude );
        //qreal amplitude = SpectrumAnalyserMultiplier * qLn( magnitude );
        //qreal amplitude = SpectrumAnalyserMultiplier * magnitude;
        //qreal amplitude = magnitude;
 
        _spectrum[i] = qBound( 0.0, amplitude, 1.0 );
    }
 
    //for( int i=0; i<10; ++i )
    //    qDebug() << "Ready: " << i << ": " << QString::number( _spectrum[i], 'f', 20 );
 
    emit spectrumReady( _spectrum );
 
    //
    // 3. Hier haben wir in '_input' ein nettes. normalisiertes Array<double>
    //
 
}
 
/*
template<class T>
void SWFFTCalc::prepareInput( QAudioBuffer& audiobuffer, T& dataType, qreal peakValue )
{
    Q_UNUSED( dataType );
 
    int frameCount = audiobuffer.frameCount();
    _input = SWDoubleVec( SWNUMSAMPLES, 0.0 );
     //_input.resize( frameCount );
    const T* frames = audiobuffer.data<T>();
 
    //for( int i=0; i<100; ++i )
    //    qDebug() << "frames: " << i << ": " <<frames[i].average() << " left:" << frames[i].left << " right: " << frames[i].right;
    int limit = qMin( frameCount, SWNUMSAMPLES );
    qDebug() << "peak: " << peakValue << "Frames: " << frameCount << " limit:" << limit;
 
    for( int i=0;  i < limit; i++)
    {
        //_input[i] = frames[i].left / peakValue;
        _input[i] = frames[i].average() / peakValue * _hannWindow[i];
        //levelLeft += abs( data[i].left ) / peakValue;
        //levelRight += abs(data[i].right)/peakValue;
    }
 
}
*/
 
 
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;
    }
 
}
 
 
/*
switch( format.sampleType() )
{
 
    case QAudioFormat::Float:
 
        qDebug() << "QAudioFormat::Float:";
        if( format.sampleSize() != 32 ) // other sample formats are not supported
            prepareInput( audiobuffer, _s32f, 0 );
        else
            prepareInput( audiobuffer, _s32f, 1.00003 );
 
        break;
 
    case QAudioFormat::SignedInt:
 
        qDebug() << "QAudioFormat::SignedInt:";
        if( format.sampleSize() == 32 )
 
#ifdef Q_OS_WIN
            prepareInput( audiobuffer, _s16s, INT_MAX );
#endif
#ifdef Q_OS_UNIX
            prepareInput( audiobuffer, _s16s, SHRT_MAX );
#endif
 
        else if( format.sampleSize() == 16 )
            prepareInput( audiobuffer, _s16s, SHRT_MAX );
        else if( format.sampleSize() == 8 )
            prepareInput( audiobuffer, _s8s, CHAR_MAX );
 
        break;
 
    case QAudioFormat::UnSignedInt:
 
        qDebug() << "QAudioFormat::UnSignedInt:";
        if( format.sampleSize() == 32 )
            prepareInput( audiobuffer, _s16s, UINT_MAX );
        else if( format.sampleSize() == 16 )
            prepareInput( audiobuffer, _s16s, USHRT_MAX );
        else if( format.sampleSize() == 8 )
            prepareInput( audiobuffer, _s16s, UCHAR_MAX );
 
    default:
 
        break;
 
} // switch
*/
 
 
 
/*
qreal getPeakValue(const QAudioFormat& format)
{
    // Note: Only the most common sample formats are supported
    if (!format.isValid())
        return qreal(0);
 
    if (format.codec() != "audio/pcm")
        return qreal(0);
 
    switch (format.sampleType()) {
    case QAudioFormat::Unknown:
        break;
    case QAudioFormat::Float:
        if (format.sampleSize() != 32) // other sample formats are not supported
            return qreal(0);
        return qreal(1.00003);
    case QAudioFormat::SignedInt:
        if (format.sampleSize() == 32)
#ifdef Q_OS_WIN
            return qreal(INT_MAX);
#endif
#ifdef Q_OS_UNIX
            return qreal(SHRT_MAX);
#endif
        if (format.sampleSize() == 16)
            return qreal(SHRT_MAX);
        if (format.sampleSize() == 8)
            return qreal(CHAR_MAX);
        break;
    case QAudioFormat::UnSignedInt:
        if (format.sampleSize() == 32)
            return qreal(UINT_MAX);
        if (format.sampleSize() == 16)
            return qreal(USHRT_MAX);
        if (format.sampleSize() == 8)
            return qreal(UCHAR_MAX);
        break;
    }
 
    return qreal(0);
}
*/