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+1));

    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.53836 - (0.46164 * cos(2 * PI * ((double)i / (((double)windowLength) * 1.0))))) ;
    }
}

void synthesistransform::SetSourceSize(int nofsource)
{

    this->nofmics = nofsource;

}


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


    complex<double>* data, * ifft_output;
    fftw_plan       plan_backward;
    int             i;

    data = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
    ifft_output = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);

    int outputcounter = 0;

    plan_backward = fftw_plan_dft_1d(this->windowsize, reinterpret_cast<fftw_complex*>(data), reinterpret_cast<fftw_complex*>(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;

                reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex);
                reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex); 

            }

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

                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][0]; // *(realpart + readIndex);
                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][1]; // *(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] = (double)(ifft_output[i].real());/* / (double)window[i]*/; // / (double)((this->windowsize / 2.0)); // / 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;

}

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


    complex<double>* data, * ifft_output;
    fftw_plan       plan_backward;
    int             i;

    data = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
    ifft_output = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);

    int outputcounter = 0;

    plan_backward = fftw_plan_dft_1d(this->windowsize, reinterpret_cast<fftw_complex*>(data), reinterpret_cast<fftw_complex*>(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;
    int zeropadding = 0;

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

       /* for (int ki = 0; ki < (this->signallength * this->nofmics); ki++)
        {
            istft_out[ki] = 0.0;
        }*/

        //float pressurelev = 0.00001 / pressurelevel[channelindex];

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

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

                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][0]; // *(realpart + readIndex);
                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][1]; // *(imagpart + readIndex);
            }

            for (i = zeropadding; i < this->windowsize / 2; i++)
            {
                readIndex = chunkPosition  * this->nofmics + i * this->nofmics + channelindex;
                // If there is not an analysis use directly
                reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex); // *window[i]; // *(float)((this->windowsize));
                reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex); // *window[i]; // *(float)((this->windowsize));
                
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
              //reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex);
              //reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex);
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////

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

                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // *(realpart + readIndex)* window[(this->windowsize / 2) - i];
                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; //  *(imagpart + readIndex)* (-1) * window[this->windowsize - i];

            }
            // 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  ; i++)
            {
                // If there is not an analysis use directly
                istft_out[outputcounter] = (double)(ifft_output[i].real()); // / windowSize; // The i term should be non-existent.
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
                // Without analysis use this one:
              //  istft_out[outputcounter] = ((double)(ifft_output[i].real()) / window[i]); // / (double)((this->windowsize));
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
                outputcounter += this->nofmics;
            }

            chunkPosition += this->overlap;
            outputcounter -= this->overlap * this->nofmics;

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

        bStop = 0;
    }

    fftw_destroy_plan(plan_backward);

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

}