Marcelo Rebonatto / Mbed 2 deprecated PM_COPIA

Versão atual 13-12-2013.

Dependencies:   EthernetInterface mbed-rtos mbed

Codes/SignalProcessor.cpp@1:238ac24e46dd, 2013-12-13 (annotated)

Committer:
rebonatto
Date:
Fri Dec 13 11:45:06 2013 +0000
Revision:
1:238ac24e46dd
Parent:
0:65c41a68b49a 
Atual 13-12-2013.

Who changed what in which revision?

UserRevisionLine numberNew contents of line
rebonatto 0:65c41a68b49a 1 /*
rebonatto 0:65c41a68b49a 2 * SignalProcessor.cpp
rebonatto 0:65c41a68b49a 3 *
rebonatto 0:65c41a68b49a 4 * Created on:
rebonatto 0:65c41a68b49a 5 * Author:
rebonatto 0:65c41a68b49a 6 */
rebonatto 0:65c41a68b49a 7
rebonatto 0:65c41a68b49a 8 #include <math.h>
rebonatto 0:65c41a68b49a 9 #include <stdlib.h>
rebonatto 0:65c41a68b49a 10 #include <stdio.h>
rebonatto 0:65c41a68b49a 11 #include <string.h>
rebonatto 0:65c41a68b49a 12
rebonatto 0:65c41a68b49a 13 #include "SignalProcessor.h"
rebonatto 0:65c41a68b49a 14
rebonatto 0:65c41a68b49a 15 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
rebonatto 0:65c41a68b49a 16 #define PI 3.14159265358979323846
rebonatto 0:65c41a68b49a 17
rebonatto 0:65c41a68b49a 18
rebonatto 0:65c41a68b49a 19 void SignalProcessor::CalculateRMSBulk(float *result)
rebonatto 0:65c41a68b49a 20 {
rebonatto 0:65c41a68b49a 21 int nChannel,nSample;
rebonatto 0:65c41a68b49a 22
rebonatto 0:65c41a68b49a 23 for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
rebonatto 0:65c41a68b49a 24 result[nChannel] = 0;
rebonatto 0:65c41a68b49a 25
rebonatto 0:65c41a68b49a 26 for(nChannel=0;nChannel<Settings::get_MaxChannels();nChannel++)
rebonatto 0:65c41a68b49a 27 {
rebonatto 0:65c41a68b49a 28 for(nSample=0;nSample<Settings::get_Samples();nSample++)
rebonatto 0:65c41a68b49a 29 {
rebonatto 0:65c41a68b49a 30
rebonatto 0:65c41a68b49a 31 unsigned short int v = Capture::GetValue(nSample, nChannel);
rebonatto 0:65c41a68b49a 32 float val = (float)v;
rebonatto 0:65c41a68b49a 33
rebonatto 0:65c41a68b49a 34 val -= Settings::get_Offset(nChannel);
rebonatto 0:65c41a68b49a 35 val /= Settings::get_Gain(nChannel);
rebonatto 0:65c41a68b49a 36 val *= val;
rebonatto 0:65c41a68b49a 37 result[nChannel] += val;
rebonatto 0:65c41a68b49a 38 }
rebonatto 0:65c41a68b49a 39 result[nChannel] /= (float)Settings::get_Samples();
rebonatto 0:65c41a68b49a 40 result[nChannel] = sqrt(result[nChannel]);
rebonatto 0:65c41a68b49a 41 }
rebonatto 0:65c41a68b49a 42 }
rebonatto 0:65c41a68b49a 43
rebonatto 0:65c41a68b49a 44 float SignalProcessor::CalculateRMS(unsigned short int *buffer,int nChannel)
rebonatto 0:65c41a68b49a 45 {
rebonatto 0:65c41a68b49a 46 float result=0;
rebonatto 0:65c41a68b49a 47 int nSample;
rebonatto 0:65c41a68b49a 48
rebonatto 0:65c41a68b49a 49 for(nSample=0;nSample<Settings::get_Samples();nSample++)
rebonatto 0:65c41a68b49a 50 {
rebonatto 0:65c41a68b49a 51
rebonatto 0:65c41a68b49a 52 unsigned short int v = buffer[nSample];
rebonatto 0:65c41a68b49a 53 float val = (float)v;
rebonatto 0:65c41a68b49a 54
rebonatto 0:65c41a68b49a 55 val -= Settings::get_Offset(nChannel);
rebonatto 0:65c41a68b49a 56 val /= Settings::get_Gain(nChannel);
rebonatto 0:65c41a68b49a 57 val *= val;
rebonatto 0:65c41a68b49a 58 result += val;
rebonatto 0:65c41a68b49a 59 }
rebonatto 0:65c41a68b49a 60 result /= (float)Settings::get_Samples();
rebonatto 0:65c41a68b49a 61 result = sqrt(result);
rebonatto 0:65c41a68b49a 62 return result;
rebonatto 0:65c41a68b49a 63 }
rebonatto 0:65c41a68b49a 64
rebonatto 0:65c41a68b49a 65 void SignalProcessor::CalculateFFT(unsigned short int *buffer,float *sen,float *cos,float *vm,int sign)
rebonatto 0:65c41a68b49a 66 {
rebonatto 0:65c41a68b49a 67 float* fft = ComplexFFT(buffer,1); //deve desalocar memoria do ptr retornado
rebonatto 0:65c41a68b49a 68 /*
rebonatto 0:65c41a68b49a 69 Mapa do vetor fft.
rebonatto 0:65c41a68b49a 70 O vetor tem 2 vezes o no. de amostras. Cada par de valores (portanto n e n+1), representam, respectivamente
rebonatto 0:65c41a68b49a 71 COS e SEN.
rebonatto 0:65c41a68b49a 72 Os dois primeiros valores reprensetam a frequencia 0Hz, portanto sao atribuidas ao valor medio.
rebonatto 0:65c41a68b49a 73 Os demais pares de valores representam a fundamental e suas harmonicas,
rebonatto 0:65c41a68b49a 74 sendo que se a fundamental for 60Hz, teremos: 60,120,180,240...
rebonatto 0:65c41a68b49a 75 Para a nossa aplicacao apenas as 12 primeiras harmonicas serao utilizadas (720Hz)
rebonatto 0:65c41a68b49a 76 */
rebonatto 0:65c41a68b49a 77 *vm = fft[0];
rebonatto 0:65c41a68b49a 78
rebonatto 0:65c41a68b49a 79 for(int i=1;i<Settings::get_MaxHarmonics()+1;i++)
rebonatto 0:65c41a68b49a 80 //for(int i=1;i<257;i++) // para testar com a FFT inversa.
rebonatto 0:65c41a68b49a 81 {
rebonatto 0:65c41a68b49a 82 cos[i-1] = fft[i*2];
rebonatto 0:65c41a68b49a 83 sen[i-1] = fft[i*2+1]*-1;
rebonatto 0:65c41a68b49a 84 }
rebonatto 0:65c41a68b49a 85
rebonatto 0:65c41a68b49a 86 free(fft);
rebonatto 0:65c41a68b49a 87 }
rebonatto 0:65c41a68b49a 88
rebonatto 0:65c41a68b49a 89
rebonatto 0:65c41a68b49a 90 float* SignalProcessor::ComplexFFT(unsigned short int* data, int sign)
rebonatto 0:65c41a68b49a 91 {
rebonatto 0:65c41a68b49a 92
rebonatto 0:65c41a68b49a 93 //variables for the fft
rebonatto 0:65c41a68b49a 94 unsigned long n,mmax,m,j,istep,i;
rebonatto 0:65c41a68b49a 95 double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
rebonatto 0:65c41a68b49a 96 float *vector;
rebonatto 0:65c41a68b49a 97 //the complex array is real+complex so the array
rebonatto 0:65c41a68b49a 98 //as a size n = 2* number of complex samples
rebonatto 0:65c41a68b49a 99 //real part is the data[index] and
rebonatto 0:65c41a68b49a 100 //the complex part is the data[index+1]
rebonatto 0:65c41a68b49a 101
rebonatto 0:65c41a68b49a 102 //new complex array of size n=2*sample_rate
rebonatto 0:65c41a68b49a 103 //if(vector==0)
rebonatto 0:65c41a68b49a 104 //vector=(float*)malloc(2*SAMPLE_RATE*sizeof(float)); era assim, define estava em Capture.h
rebonatto 0:65c41a68b49a 105 vector=(float*)malloc(2*Settings::get_Samples()*sizeof(float));
rebonatto 0:65c41a68b49a 106 memset(vector,0,2*Settings::get_Samples()*sizeof(float));
rebonatto 0:65c41a68b49a 107
rebonatto 0:65c41a68b49a 108 //put the real array in a complex array
rebonatto 0:65c41a68b49a 109 //the complex part is filled with 0's
rebonatto 0:65c41a68b49a 110 //the remaining vector with no data is filled with 0's
rebonatto 0:65c41a68b49a 111 //for(n=0; n<SAMPLE_RATE;n++)era assim, define estava em Capture.h
rebonatto 0:65c41a68b49a 112 for(n=0; n<Settings::get_Samples();n++)
rebonatto 0:65c41a68b49a 113 {
rebonatto 0:65c41a68b49a 114 if(n<Settings::get_Samples())
rebonatto 0:65c41a68b49a 115 vector[2*n]=(float)data[n];
rebonatto 0:65c41a68b49a 116 else
rebonatto 0:65c41a68b49a 117 vector[2*n]=0;
rebonatto 0:65c41a68b49a 118 vector[2*n+1]=0;
rebonatto 0:65c41a68b49a 119 }
rebonatto 0:65c41a68b49a 120
rebonatto 0:65c41a68b49a 121 //binary inversion (note that the indexes
rebonatto 0:65c41a68b49a 122 //start from 0 witch means that the
rebonatto 0:65c41a68b49a 123 //real part of the complex is on the even-indexes
rebonatto 0:65c41a68b49a 124 //and the complex part is on the odd-indexes)
rebonatto 0:65c41a68b49a 125 //n=SAMPLE_RATE << 1; //multiply by 2era assim, define estava em Capture.h
rebonatto 0:65c41a68b49a 126 n=Settings::get_Samples() << 1; //multiply by 2
rebonatto 0:65c41a68b49a 127 j=0;
rebonatto 0:65c41a68b49a 128 for (i=0;i<n/2;i+=2) {
rebonatto 0:65c41a68b49a 129 if (j > i) {
rebonatto 0:65c41a68b49a 130 SWAP(vector[j],vector[i]);
rebonatto 0:65c41a68b49a 131 SWAP(vector[j+1],vector[i+1]);
rebonatto 0:65c41a68b49a 132 if((j/2)<(n/4)){
rebonatto 0:65c41a68b49a 133 SWAP(vector[(n-(i+2))],vector[(n-(j+2))]);
rebonatto 0:65c41a68b49a 134 SWAP(vector[(n-(i+2))+1],vector[(n-(j+2))+1]);
rebonatto 0:65c41a68b49a 135 }
rebonatto 0:65c41a68b49a 136 }
rebonatto 0:65c41a68b49a 137 m=n >> 1;
rebonatto 0:65c41a68b49a 138 while (m >= 2 && j >= m) {
rebonatto 0:65c41a68b49a 139 j -= m;
rebonatto 0:65c41a68b49a 140 m >>= 1;
rebonatto 0:65c41a68b49a 141 }
rebonatto 0:65c41a68b49a 142 j += m;
rebonatto 0:65c41a68b49a 143 }
rebonatto 0:65c41a68b49a 144 //end of the bit-reversed order algorithm
rebonatto 0:65c41a68b49a 145
rebonatto 0:65c41a68b49a 146 //Danielson-Lanzcos routine
rebonatto 0:65c41a68b49a 147 mmax=2;
rebonatto 0:65c41a68b49a 148 while (n > mmax) {
rebonatto 0:65c41a68b49a 149 istep=mmax << 1;
rebonatto 0:65c41a68b49a 150 theta=sign*(2*PI/mmax);
rebonatto 0:65c41a68b49a 151 wtemp=sin(0.5*theta);
rebonatto 0:65c41a68b49a 152 wpr = -2.0*wtemp*wtemp;
rebonatto 0:65c41a68b49a 153 wpi=sin(theta);
rebonatto 0:65c41a68b49a 154 wr=1.0;
rebonatto 0:65c41a68b49a 155 wi=0.0;
rebonatto 0:65c41a68b49a 156 for (m=1;m<mmax;m+=2) {
rebonatto 0:65c41a68b49a 157 for (i=m;i<=n;i+=istep) {
rebonatto 0:65c41a68b49a 158 j=i+mmax;
rebonatto 0:65c41a68b49a 159 tempr=wr*vector[j-1]-wi*vector[j];
rebonatto 0:65c41a68b49a 160 tempi=wr*vector[j]+wi*vector[j-1];
rebonatto 0:65c41a68b49a 161 vector[j-1]=vector[i-1]-tempr;
rebonatto 0:65c41a68b49a 162 vector[j]=vector[i]-tempi;
rebonatto 0:65c41a68b49a 163 vector[i-1] += tempr;
rebonatto 0:65c41a68b49a 164 vector[i] += tempi;
rebonatto 0:65c41a68b49a 165 }
rebonatto 0:65c41a68b49a 166 wr=(wtemp=wr)*wpr-wi*wpi+wr;
rebonatto 0:65c41a68b49a 167 wi=wi*wpr+wtemp*wpi+wi;
rebonatto 0:65c41a68b49a 168 }
rebonatto 0:65c41a68b49a 169 mmax=istep;
rebonatto 0:65c41a68b49a 170 }
rebonatto 0:65c41a68b49a 171 //end of the algorithm
rebonatto 0:65c41a68b49a 172
rebonatto 0:65c41a68b49a 173 return vector;
rebonatto 0:65c41a68b49a 174 }