HR and RR measurement

Dependencies:   mbed

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);
    }
}