analysistransform.cpp

#include "analysistransform.h"



// Create a constructor for analysis transform
analysistransform::analysistransform()
{
    this->signallength = 20480;
    this->windowsize = 2048;
    this->overlap = 1024;
    this->nofmics = 32;

    int lngth = (windowsize / 2 + 1) * (signallength / (overlap));

    outputsignal = (std::complex <float> *) malloc(sizeof(std::complex <float>) * lngth);

}

// Create a constructor for analysis transform
analysistransform::analysistransform(int windowsize=2048, int signallength=20480, int hopsize=1024,int nofmics=32)
{
    this->signallength = signallength;
    this->windowsize = windowsize;
    this->overlap = hopsize;
    this->nofmics = nofmics;

    int lngth = (windowsize / 2 + 1) * (signallength / (hopsize));

    outputsignal = (std::complex <float> *) malloc(sizeof(std::complex <float>) * lngth);


}

// Create a hamming window of windowLength samples in buffer
void analysistransform::hamming(int windowLength, float* buffer) {

    for (int i = 0; i < windowLength; i++) {

        buffer[i] = 0.54 - (0.46 * cos(2 * PI * (i / ((windowLength - 1) * 1.0))));
    }
}

/*
Signal in : x[0] to x[signalLength - 1]
realpart : fft_real[0] to fft_real[N/2 -1]
imagpart : imagpart[0] to imagpart[N/2 -1]
*/
void* analysistransform::STFT(float* signal, float* realpart, float* imagpart) {


    double* data;
    fftw_complex* fft_result;
    fftw_plan       plan_forward, plan_backward;
    int             i;



    data = (double*)fftw_malloc(sizeof(double) * this->windowsize);
    fft_result = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * this->windowsize);

    int outputcounter = 0;


    plan_forward = fftw_plan_dft_r2c_1d(this->windowsize, data, fft_result, FFTW_ESTIMATE);

    // Create a hamming window of appropriate length
    float* window = (float*)malloc(sizeof(float) * this->windowsize);

    hamming(this->windowsize, window);

    int chunkPosition = 0;

    int readIndex;

    // Should we stop reading in chunks?
    int bStop = 0;

    // Process each chunk of the signal



    while (chunkPosition < this->signallength && !bStop) {

        // Copy the chunk into our buffer
        for (i = 0; i < this->windowsize; i++) {


            readIndex = chunkPosition + i;

            if (readIndex < this->signallength) {

                // Note the windowing!
                data[i] = signal[readIndex] * window[i];

            }
            else {

                // we have read beyond the signal, so zero-pad it!

                data[i] = 0.0;

                bStop = 1;

            }
        }

        // Perform the FFT on our chunk
        fftw_execute(plan_forward);


        // Copy the first (windowSize/2 + 1) data points into your spectrogram.
        // We do this because the FFT output is mirrored about the nyquist
        // frequency, so the second half of the data is redundant. This is how
        // Matlab's spectrogram routine works.


        for (i = 0; i < this->windowsize / 2; i++) {

            *(realpart + outputcounter) = fft_result[i][0];
            *(imagpart + outputcounter) = fft_result[i][1];

            outputcounter += this->nofmics;
        }


        chunkPosition += this->overlap;

    } // Excuse the formatting, the while ends here.

    fftw_destroy_plan(plan_forward);

    fftw_free(data);
    fftw_free(fft_result);
    free(window);

    return 0;

}