c-gaborator/cgaborator.cpp

#include <iostream>

#include "include/cgaborator.h"
#include "gaborator/gaborator.h"
#include <cmath>
#include <memory>

class Gaborator {

public:
    Gaborator(int block_size, double sampleRate, int bandsPerOctave, double minimumFrequency, double referenceFrequency, double maximumFrequency, int stepSize) :
    parameters(bandsPerOctave, minimumFrequency / sampleRate, referenceFrequency / sampleRate, 1.0, 1e-5),
    analyzer(parameters),
    coefs(analyzer),
    frequencyBinTimeStepSize(stepSize),
    sample_rate((int) sampleRate),
    blockSize(block_size)
    {

        //converts frequency (ff_max) in hertz to the number of bands above the min frequency
        //the ceil is used to end up at a full band
        int interesting_bands = ceil(bandsPerOctave * log(maximumFrequency/minimumFrequency)/log(2.0f));

        //since bands are ordered from high to low we are only interested in lower bands:
        //fs/2.0 is the nyquist frequency
        int total_bands = ceil(bandsPerOctave * log(sampleRate/2.0/minimumFrequency)/log(2.0f));

        latency = (int64_t) ceil(analyzer.analysis_support());
        min_band = total_bands - interesting_bands;
        t_in = 0;

        int max_band = analyzer.bandpass_bands_end();

        std::vector<float> bandcenterCache(max_band);

        for(int i = 0 ; i < max_band ; i++){
            if(i < min_band){
                bandcenterCache[i] = -1;
            }else{
                bandcenterCache[i] = analyzer.band_ff(i) * sample_rate;
            }

            if (bandcenterCache[i] > 0) {
                ++numberOfBandsCache;
                if (firstBandCache == -1) {
                    firstBandCache = i;
                }
            }
        }

        coefficientSize = (latency + 2*blockSize) / frequencyBinTimeStepSize;

        coefficients.resize(coefficientSize);
        for (auto & coefficient : coefficients){
            coefficient.resize(numberOfBandsCache);
        }

        assert(t_in == 0);

    }


	float* gaborTransform(float* audio_block, int64_t audio_block_length, size_t* return_size, size_t* slice_size) {
        resultCache.clear();

        if (audio_block == nullptr || audio_block_length == 0) { //finish
            finish();
        } else {
            analyze(audio_block, audio_block_length);
        }

        if(!resultCache.empty()){
			*return_size = resultCache.size();
			*slice_size = numberOfBandsCache;
			return resultCache.data();
        } else{
			*return_size = 0;
			*slice_size = 0;
			return nullptr;
		}
    }

    int64_t analysisSupport() const {
        return latency;
    }

    int numberOfBands() const {
        return numberOfBandsCache;
    }

    ~Gaborator()= default;


private:
    void analyze(float* audio_block, int64_t audio_block_length){
        analyzer.analyze(audio_block, t_in, t_in + audio_block_length, coefs);

        int64_t st0 = t_in - latency;
        int64_t st1 = t_in - latency + audio_block_length;

		gaborApplySlice(st0, st1);

        t_in += audio_block_length;

        int64_t t_out = t_in - latency;

        forget_before(analyzer, coefs, t_out - audio_block_length);
    }

    void finish(){
        int64_t st0 = t_in - latency;
        int64_t st1 = t_in;

        //flush all till latency spot
		gaborApplySlice(st0, st1);

        //flush remaining
        for (int i = 1; i < coefficientSize; ++i) {
            int64_t circularIndex = (mostRecentCoefficentIndex + i) % coefficientSize;

            auto& currentCoefficient = coefficients[circularIndex];

            resultCache.insert(resultCache.end(), currentCoefficient.begin(), currentCoefficient.end());
            // fill the oldest with zeros, but only the first round
			if(i <= coefficientSize) {
				std::fill(currentCoefficient.begin(), currentCoefficient.end(), 0);
			}
        }
    }

	inline void gaborApplySlice(int64_t st0, int64_t st1) {
		/*
		Following code is equivalent, but it has been inlined for performance

		gaborator::process([&](int band, int64_t audioSampleIndex, std::complex<float>& coef) {
			gaborProcessEntry(band, audioSampleIndex, coef);
		}, min_band, INT_MAX, st0, st1, coefs);
		*/

		gaborator::apply_to_slice(false, [&](int band, int64_t st, int time_step, unsigned len, const std::complex<float> *p0){
			for (unsigned int i = 0; i < len; i++) {
				gaborProcessEntry(band, st, std::abs(*p0++));
				st += time_step;
			}
		}, min_band, INT_MAX, st0, st1, coefs);
	}

    inline void gaborProcessEntry(int band, int64_t sampleIndex, float coef) {
		int64_t coefficientIndex = sampleIndex / frequencyBinTimeStepSize;
        int bandIndex = band - firstBandCache;

        // The first results have a negative audio sample index
        // ignore these
        if (coefficientIndex > 0 && bandIndex < numberOfBandsCache) {

			int64_t circularIndex = coefficientIndex % coefficientSize;

			auto& currentCoefficient = coefficients[circularIndex];

            // If a new index is reached, save the old (fixed) coefficients in the history
            // Fill the array with zeros to get the max
            if (coefficientIndex > mostRecentCoefficentIndex && coefficientIndex > coefficientSize) {
                // keep the new maximum
                mostRecentCoefficentIndex = coefficientIndex;
                // "copy" the oldest data to the history
                // the slice can be reused thanks to the oldest being filled with zeros just after
                resultCache.insert(resultCache.end(), currentCoefficient.begin(), currentCoefficient.end());
                // fill the oldest with zeros
                std::fill(currentCoefficient.begin(), currentCoefficient.end(), 0);
            }
            // due to reduction in precision (from audio sample accuracy to steps) multiple
            // magnitudes could be placed in the same stepIndex, bandIndex pair.
            // We take the maximum magnitudes value.
            currentCoefficient[bandIndex] = std::max(currentCoefficient[bandIndex], coef);
        }

    }

private:


    std::vector<float> resultCache;
    std::vector<std::vector<float>> coefficients;
    int firstBandCache = -1;
    int numberOfBandsCache = 0;
    int64_t mostRecentCoefficentIndex = 0;

    const int blockSize;
    const int64_t frequencyBinTimeStepSize;
    int64_t t_in;
    int min_band;
    const int sample_rate;
    int64_t latency;
    int64_t coefficientSize;

private:
    const gaborator::parameters parameters;
    gaborator::analyzer<float> analyzer;
    gaborator::coefs<float> coefs;
};

uintptr_t gaborator_initialize(int blockSize, double sampleRate, int bandsPerOctave, double minimumFrequency, double referenceFrequency, double maximumFrequency, int stepSize){
    return reinterpret_cast<uintptr_t>(new Gaborator(blockSize, sampleRate, bandsPerOctave, minimumFrequency, referenceFrequency,
                                                     maximumFrequency, stepSize));
}

int64_t gaborator_analysis_support(uintptr_t ptr) {
    return reinterpret_cast<Gaborator*>(ptr)->analysisSupport();
}

int gaborator_number_of_bands(uintptr_t ptr) {
    return reinterpret_cast<Gaborator*>(ptr)->numberOfBands();
}

float* gaborator_transform(uintptr_t ptr, float* audio_block, int64_t audio_block_length, size_t* return_size, size_t* slice_size){
	return reinterpret_cast<Gaborator*>(ptr)->gaborTransform(audio_block, audio_block_length, return_size, slice_size);
}

void gaborator_release(uintptr_t ptr) {
    auto g = reinterpret_cast<Gaborator*>(ptr);
    delete g;
}