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;

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

    hamming(this->windowsize, window);
}

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

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

    hamming(this->windowsize, window);
}

// Create a hamming window of windowLength samples in buffer
void analysistransform::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)))));
    }
}

/*
Signal in : x[0] to x[signalLength - 1]
nofmics : the number of the microphones
realpart : fft_real[0] to fft_real[N/2 -1]
imagpart : imagpart[0] to imagpart[N/2 -1]
*/
void* analysistransform::FFT(double* signal, double* realpart, double* imagpart) {


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

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

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

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

    int readIndex;

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

    int channelindex = 0;

    // Process each chunk of the signal

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

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

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


                readIndex = (chunkPosition + i) * this->nofmics + channelindex;

                if (readIndex < this->signallength * this->nofmics) {

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

                }
                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.
        bStop = 0;

    }
    fftw_destroy_plan(plan_forward);

    fftw_free(data);
    fftw_free(fft_result);

    return 0;

}

/*
Signal in : x[0] to x[signalLength - 1]
nofmics : the number of the microphones
realpart : fft_real[0] to fft_real[N/2 -1]
imagpart : imagpart[0] to imagpart[N/2 -1]
*/
void* analysistransform::FFT(double* signal, float* realpart, float* imagpart) {


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

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

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

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

    int readIndex;

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

    int channelindex = 0;

    // Process each chunk of the signal

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

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

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


                readIndex = (chunkPosition + i) * this->nofmics + channelindex;

                if (readIndex < this->signallength * this->nofmics) {

                    // 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++) {
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
               *(realpart + outputcounter) = (float)fft_result[i][0] / (float)((this->windowsize));
              *(imagpart + outputcounter) = (float)fft_result[i][1] / (float)((this->windowsize));
//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
              //  *(realpart + outputcounter) = ((float)fft_result[i][0] * window[i]) / (float)((this->windowsize) / 2.0);
              //  *(imagpart + outputcounter) = ((float)fft_result[i][1] * window[i]) / (float)((this->windowsize) / 2.0);
                outputcounter += this->nofmics;
            }


            chunkPosition += this->overlap;

        } // Excuse the formatting, the while ends here.
        bStop = 0;

    }
    fftw_destroy_plan(plan_forward);

    fftw_free(data);
    fftw_free(fft_result);

    return 0;

}