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 #include "limites.h"
00015 
00016 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
00017 #define PI 3.14159265358979323846F
00018 // 3.141597653564793332212487132
00019 //              3.14159265358979323846
00020                 
00021 
00022 /* Elementos vm2, under, over adicionados em 20/05/2014 por Rebonatto */
00023 /* vm2 eh o calculo do valor medio, under eh a cotagem dos valores do AD com 0 */
00024 /* over e a contagem do AD com 4095 */
00025 /* over e under sao para verificar se o processo de ajuste dos dados esta Ok e vm2 para conferir o vm da fft */ 
00026 void SignalProcessor::CalculateRMSBulk(float *result, float *vm2, int *under, int *over)
00027 {
00028     int nChannel,nSample;    
00029     
00030     for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
00031         result[nChannel] = vm2[nChannel] = under[nChannel] = over[nChannel] = 0;
00032     
00033     for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
00034     {        
00035         for(nSample=0;nSample<Settings::get_Samples();nSample++)
00036         {            
00037             unsigned short int v = Capture::GetValue(nSample, nChannel);
00038             /* novos calculos  */    
00039             vm2[nChannel] += v;
00040             if (v <= 20)
00041                 under[nChannel] = under[nChannel] + 1;
00042             if (v >= 4075)    
00043                 over[nChannel] = over[nChannel] + 1;
00044             float val = (float)v;
00045             
00046             val -= Settings::get_Offset(nChannel);
00047             val /= Settings::get_Gain(nChannel);
00048             val *= val;
00049             result[nChannel] += val;
00050         }
00051         result[nChannel] /= (float)Settings::get_Samples();
00052         result[nChannel] = sqrt(result[nChannel]);
00053         
00054         /* novos calculos */
00055         vm2[nChannel] /= (float)Settings::get_Samples();        
00056     }
00057 }
00058 
00059 float SignalProcessor::CalculateRMS(unsigned short int *buffer,int nChannel)
00060 {
00061     float result=0;
00062     int nSample;
00063     
00064     for(nSample=0;nSample<Settings::get_Samples();nSample++)
00065     {
00066         
00067         unsigned short int v = buffer[nSample];
00068         float val = (float)v;
00069                                                // cada ponto   
00070         val -= Settings::get_Offset(nChannel); // diminui o offset
00071         val /= Settings::get_Gain(nChannel);   // divide pelo ganhp
00072         val *= val;                            // eleva ao quadrado
00073         result += val;                         // soma  
00074     }
00075     result /= (float)Settings::get_Samples();  // divide pelo numero de amostras (256)
00076     result = sqrt(result);
00077     return result;
00078 }
00079  
00080 void SignalProcessor::CalculateFFT(unsigned short int *buffer,float *sen,float *cos,float *vm,int sign, int ch)
00081 {
00082     int i;
00083     //float value[256];
00084     /*
00085     printf("Tamanho float %lu\n", sizeof(float));        
00086     printf("Tamanho double %lu\n", sizeof(double));        
00087     printf("Tamanho unsigned short int %lu\n", sizeof(unsigned short int));   
00088     printf("Tamanho unsigned long %lu\n", sizeof(unsigned long));   
00089     printf("Tamanho unsigned long long %lu\n", sizeof(unsigned long long));       
00090     */   
00091     
00092     /*
00093     for(int i=0; i < Settings::get_Samples(); i++)
00094         printf("%d*",buffer[i]);
00095     printf("\n");
00096     */
00097     
00098     //printf("[0] %d %d %d %d\n", buffer[0], buffer[100], buffer[200], buffer[255]);
00099     /*
00100     for(i=0; i<Settings::get_Samples();i++)        
00101             value[i]= (float) ( (buffer[i] - Settings::get_Offset(ch)) / Settings::get_Gain(ch) );    
00102     */
00103     
00104     float* fft = ComplexFFT(buffer,1, ch);  //deve desalocar memoria do ptr retornado
00105     
00106     
00107     /*
00108         Mapa do vetor fft.
00109         O vetor tem 2 vezes o no. de amostras. Cada par de valores (portanto n e n+1), representam, respectivamente 
00110         COS e SEN.
00111         Os dois primeiros valores reprensetam a frequencia 0Hz, portanto sao atribuidas ao valor medio.
00112         Os demais pares de valores representam a fundamental e suas harmonicas,
00113         sendo que se a fundamental for 60Hz, teremos: 60,120,180,240...
00114         Para a nossa aplicacao apenas as 12 primeiras harmonicas serao utilizadas (720Hz)
00115     */
00116     
00117     //*vm = DFT(value, sen, cos);
00118     *vm = fft[0];
00119     
00120     for(int i=1;i<Settings::get_MaxHarmonics()+1;i++)
00121     {
00122         cos[i-1] = fft[i*2];
00123         sen[i-1] = fft[i*2+1];
00124     }
00125     
00126     /*
00127     for(int i=0;i<Settings::get_MaxHarmonics();i++)
00128     {
00129         printf("[%dHz]\tsen %.4f\tcos %.4f\n", (i+1)*60, sen[i], cos[i]);
00130     }
00131     */
00132     free(fft);    
00133     //printf("[3] %d %d %d %d\n", buffer[0], buffer[100], buffer[200], buffer[255]);
00134 }
00135 
00136 
00137 float* SignalProcessor::ComplexFFT(unsigned short int* data, int sign, int ch)
00138 {
00139 
00140     //variables for the fft 
00141     unsigned long n,mmax,m,j,istep,i;
00142     //double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
00143     float wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
00144     float *vector;
00145     //the complex array is real+complex so the array
00146     //as a size n = 2* number of complex samples
00147     //real part is the data[index] and
00148     //the complex part is the data[index+1]
00149 
00150     //new complex array of size n=2*sample_rate
00151     //if(vector==0)
00152     //vector=(float*)malloc(2*SAMPLE_RATE*sizeof(float)); era assim, define estava em Capture.h
00153     
00154     printf("Antes malloc\n");
00155     vector=(float*)malloc(2*Settings::get_Samples()*sizeof(float));
00156     memset(vector,0,2*Settings::get_Samples()*sizeof(float));
00157     printf("DEpois malloc\n");
00158     
00159     
00160     DisplayRAMBanks();
00161     
00162     //put the real array in a complex array
00163     //the complex part is filled with 0's
00164     //the remaining vector with no data is filled with 0's
00165     //for(n=0; n<SAMPLE_RATE;n++)era assim, define estava em Capture.h
00166             
00167     for(n=0; n<Settings::get_Samples();n++)
00168     {
00169         if(n<Settings::get_Samples()){
00170             //vector[2*n]= (float) ( (data[n] - Settings::get_Offset(ch)) / Settings::get_Gain(ch) );            
00171             vector[2*n]= (float) data[n] ;            
00172             //  printf("%.4f$", vector[2*n]);
00173             }
00174         else
00175             vector[2*n]=0;
00176         vector[2*n+1]=0;
00177     }
00178 
00179     //printf("\n");    
00180     
00181     //printf("[1] %d %d %d %d\n", data[0], data[100], data[200], data[255]);
00182     
00183     //binary inversion (note that the indexes
00184     //start from 0 witch means that the
00185     //real part of the complex is on the even-indexes
00186     //and the complex part is on the odd-indexes)
00187     //n=SAMPLE_RATE << 1; //multiply by 2era assim, define estava em Capture.h
00188     n=Settings::get_Samples() << 1; //multiply by 2
00189     j=0;
00190     for (i=0;i<n/2;i+=2) {
00191         if (j > i) {
00192             SWAP(vector[j],vector[i]);
00193             SWAP(vector[j+1],vector[i+1]);
00194             if((j/2)<(n/4)){
00195                 SWAP(vector[(n-(i+2))],vector[(n-(j+2))]);
00196                 SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]);
00197             }
00198         }
00199         m=n >> 1;
00200         while (m >= 2 && j >= m) {
00201             j -= m;
00202             m >>= 1;
00203         }
00204         j += m;
00205     }
00206     //end of the bit-reversed order algorithm
00207 
00208     //Danielson-Lanzcos routine
00209     mmax=2;
00210     while (n > mmax) {
00211         istep=mmax << 1;
00212         theta=sign*(2*PI/mmax);
00213         wtemp=sin(0.5*theta);
00214         wpr = -2.0*wtemp*wtemp;
00215         wpi=sin(theta);
00216         wr=1.0;
00217         wi=0.0;
00218         for (m=1;m<mmax;m+=2) {
00219             for (i=m;i<=n;i+=istep) {
00220                 j=i+mmax;
00221                 tempr=wr*vector[j-1]-wi*vector[j];
00222                 tempi=wr*vector[j]+wi*vector[j-1];
00223                 vector[j-1]=vector[i-1]-tempr;
00224                 vector[j]=vector[i]-tempi;
00225                 vector[i-1] += tempr;
00226                 vector[i] += tempi;
00227             }
00228             wr=(wtemp=wr)*wpr-wi*wpi+wr;
00229             wi=wi*wpr+wtemp*wpi+wi;
00230         }
00231         mmax=istep;
00232     }
00233     //end of the algorithm
00234 
00235     /*
00236     // Ajustes a FFT
00237     for(i = 0; i < Settings::get_Samples()*2; i++ ){          
00238         vector[i] = (float) ((2 * vector[i]) / Settings::get_Samples() );
00239         /*
00240         if (i % 2 == 1)
00241             vector[i] = vector[i] * -1;
00242         
00243     }
00244     */
00245     //printf("[2] %d %d %d %d\n", data[0], data[100], data[200], data[255]);
00246     
00247     return vector;
00248 }
00249 
00250 /*
00251 float SignalProcessor::DFT(float *data, float *seno, float *coss){
00252     int i, j;
00253     
00254     for(i=0; i < Settings::get_MaxHarmonics()+1; i++)
00255         seno[i] = coss[i] = 0;
00256     
00257     for(i=0; i < Settings::get_Samples(); i++){
00258         for(j = 0; j < Settings::get_MaxHarmonics()+1; j++ ){          
00259             coss[j] += (data[i] * (cos( (2 * PI * i * j) / Settings::get_Samples() ) ) ) ;
00260             seno[j] += (data[i] * (sin( (2 * PI * i * j) / Settings::get_Samples() ) ) ) ;
00261         }
00262     }
00263        
00264     for(j = 1; j < Settings::get_MaxHarmonics()+1; j++ ){          
00265         coss[j] = 2 * coss[j] / Settings::get_Samples();
00266         seno[j] = 2 * seno[j] / Settings::get_Samples() ;
00267     }
00268     return (float) (coss[0] / Settings::get_Samples()) + (seno[0] / Settings::get_Samples());
00269 }
00270 */