synthesistransform.cpp

#include "synthesistransform.h"




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

    window = (double*)malloc(sizeof(double) * this->windowsize);

    hamming(this->windowsize, window);
}

// Create a constructor for analysis transform
synthesistransform::synthesistransform(int windowsize, int signallength, int hopsize, int nofmics)
{
    this->signallength = signallength;
    this->windowsize = windowsize;
    this->overlap = hopsize;
    this->nofmics = nofmics;

    window = (double*)malloc(sizeof(double) * this->windowsize);

    hamming(this->windowsize, window);
}

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

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

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


void* synthesistransform::IFFT(double* istft_out, double* realpart, double* imagpart) {


    fftw_complex* data, * ifft_output;
    fftw_plan       plan_backward;
    int             i;

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

    int outputcounter = 0;

    plan_backward = fftw_plan_dft_1d(this->windowsize, data, ifft_output, FFTW_BACKWARD, FFTW_ESTIMATE);

    int readIndex;

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

    // Process each chunk of the signal
    int channelindex = 0;

    for (channelindex = 0; channelindex < this->nofmics; channelindex++)
    {
        int outputcounter = channelindex;
        int chunkPosition = 0;

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

            for (i = 0; i < this->windowsize / 2; i++)
            {
                readIndex = (chunkPosition + i) * this->nofmics + channelindex;

                data[i][0] = *(realpart + readIndex);
                data[i][1] = *(imagpart + readIndex);

            }

            for (i = this->windowsize / 2; i < this->windowsize; i++)
            {
                readIndex = (chunkPosition + (i - this->windowsize / 2)) * this->nofmics + channelindex;

                data[i][0] = 0.0; // *(realpart + readIndex);
                data[i][1] = 0.0; // *(imagpart + readIndex);

            }
            // Perform the FFT on our chunk
            fftw_execute(plan_backward);

            // 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.

            //outputcounter = chunkPosition * this->nofmics + channelindex;

            for (i = 0; i < this->windowsize/2; i++)
            {
                istft_out[outputcounter] = (ifft_output[i][0] / window[i]) / (this->windowsize / 2); // / windowSize; // The i term should be non-existent.

                outputcounter += this->nofmics;
            }

            chunkPosition += this->overlap;

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

        bStop = 0;
    }

    fftw_destroy_plan(plan_backward);

    fftw_free(data);
    fftw_free(ifft_output);
    return 0;

}