Muhammad Khairul Nizar Bin Khairuddin
/
HRandRR
HR and RR measurement
main.cpp
- Committer:
- nizarpang
- Date:
- 2017-06-12
- Revision:
- 1:16fb5b4d1024
- Parent:
- 0:99c4524cab70
File content as of revision 1:16fb5b4d1024:
#include "mbed.h" #include "math.h" #define PI 3.14159265358979 //円周率 #define N 1024 //FFT点数 #define Step 0.015 //Step time #define M 512 #define S 1 #include "TextLCD_SB1602E.h" LocalFileSystem local("local"); DigitalOut myled1(LED1); DigitalOut myled2(LED2); DigitalOut myled3(LED3); DigitalOut myled4(LED4); AnalogIn myinput(p15); DigitalOut rst1(p11); I2C i2c( p28, p27 ); // sda, scl TextLCD_SB1602E lcd( &i2c ); InterruptIn takuto_sw(p18); Ticker timer; FILE *fp1;//File Pointer1 FILE *fp2;//File Pointer2 FILE *fp3; int i,j,cnt; int e; float x0[N]; float xIm[N]; float q0[N]; float qIm[N]; float s0[N]; float sIm[N]; char fname1[8]; char fname2[8]; char fname3[8]; float r; float k; //int k; int d; int num; float T; float max; /*least square*/ void sai1(float s0[],float keisu[]) { //int i; float a0,a1; float A00,A01,A02,A11,A12; A00=A01=A02=A11=A12=0.0; for (i=0;i<N;i++) { A00+=1.0; A01+=(int)i*Step; A02+=s0[i]; A11+=(int)i*Step*(int)i*Step; A12+=(int)i*Step*s0[i]; } a0=(A00*A12-A01*A02)/(A00*A11-A01*A01); a1=(A02*A11-A01*A12)/(A00*A11-A01*A01); keisu[0]=a0; keisu[1]=a1; } /*FFT*/ void fft(float x0[],float xIm[]){ int ip,m,me,me1,n,nv2,p,k; double uRe,uIm,vRe,vIm,wRe,wIm,tRe,tIm; p=10; n= 1 << p; nv2=n/2; //ビット逆順 j=0; for(i=0; i<n-1; i++){ if(j>i){ tRe=x0[j]; x0[j]=x0[i]; x0[i]=tRe; tIm=xIm[j]; xIm[j]=xIm[i]; xIm[i]=tIm; } k=nv2; while(j >= k){ j -= k; k /=2; } j+=k; } //main for(m=1;m<=p;m++){ me=1<<m; me1=me/2; uRe=1.0; uIm=0.0; wRe=cos(PI/me1); wIm=-sin(PI/me1); for(j=0; j<me1;j++){ for(i=j;i<n;i+=me){ ip=i+me1; tRe=x0[ip]*uRe-xIm[ip]*uIm; tIm=x0[ip]*uIm+xIm[ip]*uRe; x0[ip]=x0[i]-tRe; xIm[ip]=xIm[i]-tIm; x0[i]+=tRe; xIm[i]+=tIm; } vRe=uRe*wRe-uIm*wIm; vIm=uRe*wIm+uIm*wRe; uRe=vRe; uIm=vIm; } } } /*IFFT*/ void ifft(float q0[],float qIm[]){ int p,n; p=10; n=1<<p; for(i=0;i<n;i++) qIm[i]= -1*qIm[i]; fft(q0,qIm); for(i=0;i<N;i++){ q0[i]/=n; qIm[i]/=(-n); } } void aifft(float s0[],float sIm[]){ int p,n; p=10; n=1<<p; for(i=0;i<n;i++) sIm[i]= -1*sIm[i]; fft(s0,sIm); for(i=0;i<N;i++){ s0[i]/=n; sIm[i]/=(-n); } } /*Timer Fanction*/ void tflip() { float Ain; Ain=0.0; Ain=(float)myinput.read()*3.3; x0[i]=Ain; if(i<N) i++; myled3=1; myled4=1; myled2=1; } /*極大検知*/ void result(float q0[]){ d=0; e=0; for(i=1; i < N; i++){ if(q0[i-1] < q0[i] && q0[i+1] < q0[i]){ x0[d]=i; d=d+1; } } for(i=0; i<d; i++){ // fprintf( fp2, "%f\r\n",x0[i] ); } for(i=0; i<d; i++){ if(k < (x0[i+1]-x0[i])*0.015){ e=e+1; } } //fprintf( fp3, "%d\r\n",e ); r=0; for(i=0; i<e; i++){ r=x0[i+1]-x0[i]+r; } r=r*Step; // fprintf( fp3, "%f\r\n",r ); r=r/e; //fprintf( fp3, "%f\r\n",r ); r=60/r; //fprintf( fp3, "%f\r\n",r ); } void flip() { float keisu[2]; //double max1[N]; rst1=1; /*create txt File(for derivation of pulse rate)*/ sprintf(fname1,"/local/WP%d.txt",cnt); //Filename of Rest_PulseRate if ( NULL == (fp1 = fopen( fname1, "w" )) ) error( "" ); /*create txt File(for derivation of pulse rate)*/ sprintf(fname2,"/local/WS%d.txt",cnt); //Filename of Rest_PulseRate if ( NULL == (fp2 = fopen( fname2, "w" )) ) error( "" ); /*create txt File(for derivation of pulse rate)*/ sprintf(fname3,"/local/WT%d.txt",cnt); //Filename of Rest_PulseRate if ( NULL == (fp3 = fopen( fname3, "w" )) ) error( "" ); wait(2); lcd.clear(); /*functions formatting*/ for ( i = 0; i < N; i++ ) { x0[i]=0.0; xIm[i]=0.0; } /*measuring start and import Ain to x[N]*/ i=0; j=0; timer.attach(&tflip,Step); wait(Step*N); timer.detach(); /*writing to a text to x0*/ for( i = 0; i < N; i++ ){ fprintf( fp1, "%f\r\n",x0[i] ); myled2=0; } /*FFT*/ fft(x0,xIm); for(i=0;i<N;i++){ //P[i]=sqrt(x0[i]*x0[i]+xIm[i]*xIm[i]); // fprintf( fp2, "%f\r\n",P[i]); // fprintf( fp3, "%f\r\n",xIm[i]); } myled4=0; /*呼吸数抽出*/ for(i=0;i<N;i++){ q0[i]=0; qIm[i]=0; } for(i=0;i<8;i++){ q0[i]=x0[i]; qIm[i]=xIm[i]; } for(i=1017;i<1025;i++){ q0[i]=x0[i]; qIm[i]=xIm[i]; } /*心音抽出*/ for(i=0;i<N;i++){ s0[i]=0; sIm[i]=0; } for(i=0;i<1;i++){ s0[i]=x0[i]; sIm[i]=xIm[i]; } for(i=14;i<31;i++){ s0[i]=x0[i]; sIm[i]=xIm[i]; } for(i=994;i<1011;i++){ s0[i]=x0[i]; sIm[i]=xIm[i]; } /*IFFT(呼吸)*/ ifft(q0,qIm); for(i=0;i<N;i++){ fprintf( fp2, "%f\r\n",q0[i]); } k=1.8; result(q0); //fprintf( fp3, "%f\r\n",r ); myled4 = 0; /*least square*/ /* keisu[0]=keisu[1]=0.0; sai1(q0,keisu); for(i=0;i<N;i++){ q0[i]=q0[i]-(int)i*Step*keisu[0]-keisu[1]; // fprintf( fp3, "%f\r\n",q0[i]); } */ /*極大値検知*/ /* j=0; for(i=1;i<(N-1);i++){ if(q0[i+1]<q0[i] && q0[i]>q0[i-1]){ max1[j]=i; j=j+1; } } fprintf( fp2, "%d\r\n",j ); */ /*自己相関関数呼吸*/ /* // coll(s0); // float r; for(j=0;j <= M;j++){ r=0; for(i=1;i<=M;i++){ r=r+q0[i]*q0[i+j]; } t[j] = r; } if(t[0]){ for(j=0; j <= M; j++){ t[j]=t[j]/t[0]; } } for(i=0;i<M;i++){ fprintf( fp2, "%f\r\n",t[i] ); } */ /*極大検知*/ /* T=max=bpm=0.0; num=intbpm=0; for( i=(0.5/(float)Step); i<(5.5/(float)Step); i++ ){ if(t[i]>max){ max=t[i]; num=i; } } T=(int)num*Step; bpm=60.0/T; */ lcd.clear(); lcd.printf( 0, "RR %f\r", r ); /*極大検知*/ // result(q0); //fprintf( fp2, "%d\r\n",d ); /*IFFT(心音)*/ aifft(s0,sIm); for(i=0;i<N;i++){ //fprintf( fp3, "%f\r\n",s0[i]); } /*least square*/ keisu[0]=keisu[1]=0.0; sai1(s0,keisu); for(i=0;i<N;i++){ s0[i]=s0[i]-(int)i*Step*keisu[0]-keisu[1]; // fprintf( fp3, "%f\r\n",s0[i]); } /*自己相関関数心拍*/ // coll(s0); // float r; for(j=0;j <= M;j++){ r=0; for(i=1;i<=M;i++){ r=r+s0[i]*s0[i+j]; } x0[j] = r; } if(x0[0]){ for(j=0; j <= M; j++){ x0[j]=x0[j]/x0[0]; } } for(i=0;i<M;i++){ fprintf( fp3, "%f\r\n",x0[i] ); } /*極大検知*/ T=max=0.0; num=0; for( i=(0.5/(float)Step); i<(1.4/(float)Step); i++ ){ if(x0[i]>max){ max=x0[i]; num=i; } } max=0.0; T=(int)num*Step; max=60.0/T; //fprintf( fp3, "%f\r\n",bpm ); //k=0.5; //result(x0); //fprintf( fp2, "%f\r\n",r ); //lcd.clear(); //lcd.printf( 0, "RR %f\r", d ); //result(t); //fprintf( fp2, "%d\r\n",e ); lcd.printf( 1, "HR %f\r", max ); fclose( fp1 ); fclose( fp2 ); fclose( fp3 ); myled3 = 0; } int main() { cnt=0; lcd.printf( 0, "PUSH RED\n"); takuto_sw.fall(&flip); // attach the address of the flip function to the falling edge NVIC_SetPriority( EINT3_IRQn, 255);//Interrupt Priority setting to "InterruptIn" while(1) { // wait around, interrupts will interrupt this! myled1 = 1; wait(0.1); myled1 = 0; wait(2.9); } }