SignalProcessor.cpp

00001 /*
00002  * SignalProcessor.cpp
00003  *
00004  *  Created on: 
00005  *      Author: 
00006  */
00007  
00008 #include <math.h>
00009 #include <stdlib.h>
00010 #include <stdio.h>
00011 #include <string.h>
00012 
00013 #include "SignalProcessor.h"
00014 
00015 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
00016 #define PI      3.14159265358979323846
00017 
00018 
00019 void SignalProcessor::CalculateRMSBulk(float *result)
00020 {
00021     int nChannel,nSample;
00022     
00023     for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
00024         result[nChannel] = 0;
00025     
00026     for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
00027     {
00028         for(nSample=0;nSample<Settings::get_Samples();nSample++)
00029         {
00030             
00031             unsigned short int v = Capture::GetValue(nSample, nChannel);
00032             float val = (float)v;
00033             
00034             val -= Settings::get_Offset(nChannel);
00035             val /= Settings::get_Gain(nChannel);
00036             val *= val;
00037             result[nChannel] += val;
00038         }
00039         result[nChannel] /= (float)Settings::get_Samples();
00040         result[nChannel] = sqrt(result[nChannel]);
00041     }
00042 }
00043 
00044 float SignalProcessor::CalculateRMS(unsigned short int *buffer,int nChannel)
00045 {
00046     float result=0;
00047     int nSample;
00048     
00049     for(nSample=0;nSample<Settings::get_Samples();nSample++)
00050     {
00051         
00052         unsigned short int v = buffer[nSample];
00053         float val = (float)v;
00054         
00055         val -= Settings::get_Offset(nChannel);
00056         val /= Settings::get_Gain(nChannel);
00057         val *= val;
00058         result += val;
00059     }
00060     result /= (float)Settings::get_Samples();
00061     result = sqrt(result);
00062     return result;
00063 }
00064  
00065 void SignalProcessor::CalculateFFT(unsigned short int *buffer,float *sen,float *cos,float *vm,int sign)
00066 {
00067     float* fft = ComplexFFT(buffer,1);  //deve desalocar memoria do ptr retornado
00068     /*
00069         Mapa do vetor fft.
00070         O vetor tem 2 vezes o no. de amostras. Cada par de valores (portanto n e n+1), representam, respectivamente 
00071         COS e SEN.
00072         Os dois primeiros valores reprensetam a frequencia 0Hz, portanto sao atribuidas ao valor medio.
00073         Os demais pares de valores representam a fundamental e suas harmonicas,
00074         sendo que se a fundamental for 60Hz, teremos: 60,120,180,240...
00075         Para a nossa aplicacao apenas as 12 primeiras harmonicas serao utilizadas (720Hz)
00076     */
00077     *vm = fft[0];
00078     
00079     for(int i=1;i<Settings::get_MaxHarmonics()+1;i++)
00080     //for(int i=1;i<257;i++) // para testar com a FFT inversa.
00081     {
00082         cos[i-1] = fft[i*2];
00083         sen[i-1] = fft[i*2+1]*-1;
00084     }
00085     
00086     free(fft);    
00087 }
00088 
00089 
00090 float* SignalProcessor::ComplexFFT(unsigned short int* data, int sign)
00091 {
00092 
00093     //variables for the fft
00094     unsigned long n,mmax,m,j,istep,i;
00095     double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
00096     float *vector;
00097     //the complex array is real+complex so the array
00098     //as a size n = 2* number of complex samples
00099     //real part is the data[index] and
00100     //the complex part is the data[index+1]
00101 
00102     //new complex array of size n=2*sample_rate
00103     //if(vector==0)
00104     //vector=(float*)malloc(2*SAMPLE_RATE*sizeof(float)); era assim, define estava em Capture.h
00105     vector=(float*)malloc(2*Settings::get_Samples()*sizeof(float));
00106     memset(vector,0,2*Settings::get_Samples()*sizeof(float));
00107 
00108     //put the real array in a complex array
00109     //the complex part is filled with 0's
00110     //the remaining vector with no data is filled with 0's
00111     //for(n=0; n<SAMPLE_RATE;n++)era assim, define estava em Capture.h
00112     for(n=0; n<Settings::get_Samples();n++)
00113     {
00114         if(n<Settings::get_Samples())
00115             vector[2*n]=(float)data[n];
00116         else
00117             vector[2*n]=0;
00118         vector[2*n+1]=0;
00119     }
00120 
00121     //binary inversion (note that the indexes
00122     //start from 0 witch means that the
00123     //real part of the complex is on the even-indexes
00124     //and the complex part is on the odd-indexes)
00125     //n=SAMPLE_RATE << 1; //multiply by 2era assim, define estava em Capture.h
00126     n=Settings::get_Samples() << 1; //multiply by 2
00127     j=0;
00128     for (i=0;i<n/2;i+=2) {
00129         if (j > i) {
00130             SWAP(vector[j],vector[i]);
00131             SWAP(vector[j+1],vector[i+1]);
00132             if((j/2)<(n/4)){
00133                 SWAP(vector[(n-(i+2))],vector[(n-(j+2))]);
00134                 SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]);
00135             }
00136         }
00137         m=n >> 1;
00138         while (m >= 2 && j >= m) {
00139             j -= m;
00140             m >>= 1;
00141         }
00142         j += m;
00143     }
00144     //end of the bit-reversed order algorithm
00145 
00146     //Danielson-Lanzcos routine
00147     mmax=2;
00148     while (n > mmax) {
00149         istep=mmax << 1;
00150         theta=sign*(2*PI/mmax);
00151         wtemp=sin(0.5*theta);
00152         wpr = -2.0*wtemp*wtemp;
00153         wpi=sin(theta);
00154         wr=1.0;
00155         wi=0.0;
00156         for (m=1;m<mmax;m+=2) {
00157             for (i=m;i<=n;i+=istep) {
00158                 j=i+mmax;
00159                 tempr=wr*vector[j-1]-wi*vector[j];
00160                 tempi=wr*vector[j]+wi*vector[j-1];
00161                 vector[j-1]=vector[i-1]-tempr;
00162                 vector[j]=vector[i]-tempi;
00163                 vector[i-1] += tempr;
00164                 vector[i] += tempi;
00165             }
00166             wr=(wtemp=wr)*wpr-wi*wpi+wr;
00167             wi=wi*wpr+wtemp*wpi+wi;
00168         }
00169         mmax=istep;
00170     }
00171     //end of the algorithm
00172 
00173     return vector;
00174 }