File content as of revision 0:65c41a68b49a:

/*
 * SignalProcessor.cpp
 *
 *  Created on: 
 *      Author: 
 */
 
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "SignalProcessor.h"

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define PI      3.14159265358979323846


void SignalProcessor::CalculateRMSBulk(float *result)
{
    int nChannel,nSample;
    
    for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
        result[nChannel] = 0;
    
    for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
    {
        for(nSample=0;nSample<Settings::get_Samples();nSample++)
        {
            
            unsigned short int v = Capture::GetValue(nSample, nChannel);
            float val = (float)v;
            
            val -= Settings::get_Offset(nChannel);
            val /= Settings::get_Gain(nChannel);
            val *= val;
            result[nChannel] += val;
        }
        result[nChannel] /= (float)Settings::get_Samples();
        result[nChannel] = sqrt(result[nChannel]);
    }
}

float SignalProcessor::CalculateRMS(unsigned short int *buffer,int nChannel)
{
    float result=0;
    int nSample;
    
    for(nSample=0;nSample<Settings::get_Samples();nSample++)
    {
        
        unsigned short int v = buffer[nSample];
        float val = (float)v;
        
        val -= Settings::get_Offset(nChannel);
        val /= Settings::get_Gain(nChannel);
        val *= val;
        result += val;
    }
    result /= (float)Settings::get_Samples();
    result = sqrt(result);
    return result;
}
 
void SignalProcessor::CalculateFFT(unsigned short int *buffer,float *sen,float *cos,float *vm,int sign)
{
    float* fft = ComplexFFT(buffer,1);  //deve desalocar memoria do ptr retornado
    /*
        Mapa do vetor fft.
        O vetor tem 2 vezes o no. de amostras. Cada par de valores (portanto n e n+1), representam, respectivamente 
        COS e SEN.
        Os dois primeiros valores reprensetam a frequencia 0Hz, portanto sao atribuidas ao valor medio.
        Os demais pares de valores representam a fundamental e suas harmonicas,
        sendo que se a fundamental for 60Hz, teremos: 60,120,180,240...
        Para a nossa aplicacao apenas as 12 primeiras harmonicas serao utilizadas (720Hz)
    */
    *vm = fft[0];
    
    for(int i=1;i<Settings::get_MaxHarmonics()+1;i++)
    //for(int i=1;i<257;i++) // para testar com a FFT inversa.
    {
        cos[i-1] = fft[i*2];
        sen[i-1] = fft[i*2+1]*-1;
    }
    
    free(fft);    
}


float* SignalProcessor::ComplexFFT(unsigned short int* data, int sign)
{

    //variables for the fft
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
    float *vector;
    //the complex array is real+complex so the array
    //as a size n = 2* number of complex samples
    //real part is the data[index] and
    //the complex part is the data[index+1]

    //new complex array of size n=2*sample_rate
    //if(vector==0)
    //vector=(float*)malloc(2*SAMPLE_RATE*sizeof(float)); era assim, define estava em Capture.h
    vector=(float*)malloc(2*Settings::get_Samples()*sizeof(float));
    memset(vector,0,2*Settings::get_Samples()*sizeof(float));

    //put the real array in a complex array
    //the complex part is filled with 0's
    //the remaining vector with no data is filled with 0's
    //for(n=0; n<SAMPLE_RATE;n++)era assim, define estava em Capture.h
    for(n=0; n<Settings::get_Samples();n++)
    {
        if(n<Settings::get_Samples())
            vector[2*n]=(float)data[n];
        else
            vector[2*n]=0;
        vector[2*n+1]=0;
    }

    //binary inversion (note that the indexes
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes
    //and the complex part is on the odd-indexes)
    //n=SAMPLE_RATE << 1; //multiply by 2era assim, define estava em Capture.h
    n=Settings::get_Samples() << 1; //multiply by 2
    j=0;
    for (i=0;i<n/2;i+=2) {
        if (j > i) {
            SWAP(vector[j],vector[i]);
            SWAP(vector[j+1],vector[i+1]);
            if((j/2)<(n/4)){
                SWAP(vector[(n-(i+2))],vector[(n-(j+2))]);
                SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]);
            }
        }
        m=n >> 1;
        while (m >= 2 && j >= m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    //end of the bit-reversed order algorithm

    //Danielson-Lanzcos routine
    mmax=2;
    while (n > mmax) {
        istep=mmax << 1;
        theta=sign*(2*PI/mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        for (m=1;m<mmax;m+=2) {
            for (i=m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*vector[j-1]-wi*vector[j];
                tempi=wr*vector[j]+wi*vector[j-1];
                vector[j-1]=vector[i-1]-tempr;
                vector[j]=vector[i]-tempi;
                vector[i-1] += tempr;
                vector[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
    //end of the algorithm

    return vector;
}