아두이노(arduino) 심박센서 (heart rate sensor) 심박수 측정 example code

하드 정리 중 예전에 만든 아두이노에서 심박센서를 사용하여 심박수 측정하는 소스 코드가 있어, 이를 블로그 글로 남긴다.

인터넷 쇼핑몰에서 쉽게 구할 수 있는 저렴함 센서가 사용 되었으며, 연결은 아래 그림과 같이 연결 되었다. 
  
심박센서의 하트 부위에 손가락을 올려 놓고, 아날로그 A0값을 읽으면 심박의 움직임에 따라 A0값이 변하는 동작이 간단한 센서이다. 아래 그림의 왼쪽은 일정 시간 간격 마다 A0값을 읽어 시리얼 플로터로 본 심박에 따른 데이터이고, 오른쪽은 평균 필터를 사용해 노이즈를 제거한 후 시리얼 플로터로 본 데이터이다.
 
심박수는 수집된 데이터에서 threshold를 넘는 peak들의 시간 간격을 측정하여 계산한다. 
peak위치를 찾는 방법은 여러 가지 방법이 있으나, 소스에 구현된 방법은 신호의 기울기가 +에서 0또는 -로 바뀌는 위치를 peak으로 선정하도록 하였다.
심장이 t시간동안 N회 뛰는 경우 심박수는 heart rate bpm = (N/t)*60 을 표현 할 수 있다.
 
threshold는 데이터를 보고 임의로 정했으며, 측정 데이터에 따라 자동으로 정하고 싶을 때는 cfar같은 방법을 사용하면 된다. 
소스는 아래와 같다.


#define heartbit_adc_threshold      511
#define meanfilter_size             5
#define heartbit_tick_buffer_size   10
#define measure_buffer_size       3

float meanfilter[meanfilter_size] = {0,};
float measure_buffer[measure_buffer_size] = {0,};
unsigned long heartbit_tick[heartbit_tick_buffer_size] = {0,};
float heartrate_bpm = 0;
float prev_heartrate_bpm = 0;
unsigned long last_heartrate_bpm_update_tick = 0;

void setup() {
  // put your setup code here, to run once:
  Serial.begin(9600);
  for (int i = 0; i < meanfilter_size; i++)
  {
    meanfilter[i] = 0;
  }

  clear_measure_buffer();

  heartrate_bpm = 0;
  prev_heartrate_bpm = 0;
}

float filtering(float val)
{
  int i = 0;
  float mean = 0;
  int cnt = 0;
  float max_val = 0;
  float min_val = 1024;

  for (i = 1; i < meanfilter_size; i++)
  {
    meanfilter[i - 1] = meanfilter[i];
  }

  meanfilter[meanfilter_size - 1] = val;

  for (int i = 0; i < meanfilter_size; i++)
  {
    if (meanfilter[i])
    {
      mean += meanfilter[i];
      cnt++;
      if (max_val < meanfilter[i])
      {
        max_val = meanfilter[i];
      }

      if (min_val > meanfilter[i])
      {
        min_val = meanfilter[i];
      }
    }
  }

  if (max_val)
  {
    mean -= max_val;
    cnt--;
  }

  if (min_val != 1024 && min_val != max_val)
  {
    mean -= min_val;
    cnt--;
  }

  if (cnt > 0)
  {
    return mean / cnt;
  }

  return val;
}

float get_trend(float *signals, int n)
{
  float x = 0, y = 0, a = 0, b = 0;
  float sy = 0.0,
        sxy = 0.0,
        sxx = 0.0;
  int i = 0;

  x = (-n / 2.0 + 0.5);

  while (i < n)
  {
    y = signals[i];
    sy += y;
    sxy += x * y;
    sxx += x * x;

    i++;
    x += 1.0;
  }

  b = sxy / sxx;
  a = sy / n;

  return b;
}

void clear_measure_buffer()
{
  for (int i = 0; i < heartbit_tick_buffer_size; i++)
  {
    heartbit_tick[i] = 0;
  }

}

float prev_trend = 0;
boolean found_high = false;

float measure_heartbitrate(float val, float threshold)
{
  int i = 0;
  float cur_trend = 0;

  found_high = false;

  for (i = 1; i < measure_buffer_size; i++)
  {
    measure_buffer[i - 1] = measure_buffer[i];
  }

  measure_buffer[measure_buffer_size - 1] = val;

  cur_trend = get_trend(measure_buffer, measure_buffer_size);


  if (prev_trend > 0 && cur_trend <= 0 && val > threshold )
  {
    found_high = true;
  }
  prev_trend = cur_trend;


  if (found_high)
  {
    for (i = 1; i < heartbit_tick_buffer_size; i++)
    {
      heartbit_tick[i - 1] = heartbit_tick[i];
    }

    heartbit_tick[heartbit_tick_buffer_size - 1] = millis();

    if (heartbit_tick[0] != 0 && heartbit_tick[heartbit_tick_buffer_size - 1] != 0)
    {
      float timedur = (float)heartbit_tick[heartbit_tick_buffer_size - 1] - (float)heartbit_tick[0];
      float bpm = (((heartbit_tick_buffer_size-1) * 1000) / timedur) * 60;

      if(bpm < 50 || bpm > 120)
      {
        // wrong bpm...
        clear_measure_buffer();
        return 0;
      }

      return bpm;
    }
  }

  return 0;
}

void loop() {
  // put your main code here, to run repeatedly:

  float cur_measure_val = 0;
  unsigned long cur_time_ms = millis();

  float heartbeat_level = filtering(analogRead(A0));

  cur_measure_val = measure_heartbitrate(heartbeat_level, heartbit_adc_threshold);
  if (cur_measure_val)
  {
    heartrate_bpm = cur_measure_val;
    last_heartrate_bpm_update_tick = cur_time_ms;
  }
  else
  {
    if (cur_time_ms - last_heartrate_bpm_update_tick >= 10 * 1000)
    {
      heartrate_bpm = 0;
      last_heartrate_bpm_update_tick = cur_time_ms;
      clear_measure_buffer();
    }
  }

  /*
  */
#if 0
  Serial.print(heartbeat_level);
  Serial.print(" ");
  //Serial.print(heartrate_bpm);
  //Serial.print(" ");

  if (found_high)
  {
    Serial.print(heartbeat_level);
  }
  else
  {
    Serial.print(heartbit_adc_threshold);

  }

  Serial.print("\n");
  #else
    if (prev_heartrate_bpm != heartrate_bpm)
  {
    if (heartrate_bpm)
    {
      Serial.print("HeartBeat Rate :");
      Serial.print(heartrate_bpm);
      Serial.println(" bpm");
    }
    else
    {
      Serial.println("HeartBeat Rate : ..... ");
    }
    
    prev_heartrate_bpm = heartrate_bpm;
  }
#endif

  /**/
  cur_time_ms = millis() - cur_time_ms;
  if (cur_time_ms < 10)
  {
    delay(10 - cur_time_ms);
  }
}


아래 그림은 실제 심박수를 측정한 결과다.


댓글

댓글 쓰기

이 블로그의 인기 게시물

쉽게 설명한 파티클 필터(particle filter) 동작 원리와 예제

이미지
파티클 필터(particle filter)는 칼만 필터(kalman filter) 와 마찬가지로 노이즈가 있는 환경에서 측정된 데이터를 필터를 사용해 실제 위치를 추정하는 도구다. 파티클 필터(particle filter)는 보통 가우시안 분포가 아닌 측정 데이터를 다루기 위해 사용된다고 한다. 물론 가우시안 분포의 데이터에서 사용하지 말라는 건 아니다.  본 포스트에서는 파티클 필터(particle filter)의 어려운 수학적인 내용은 제하고, 쉽게 예를 들며 필터의 동작원리에 대해 알아보았다. 파티클 필터(particle filter)에 대해 검색해 보면 아래와 같은 그림을 많이 보게 된다. 아래 그림은 파티클 필터의 estimation cycle을 도식화한 것이다. <’Real-Time Tracking of Multiple Moving Objects Using Particle Filters and Probabilistic Data Association’ original scientific paper 中> 위 그림의 검은색 원들은 particle을 의미한다. 이 particle에는 보통 위치 데이터와 weight가 포함된다. 검은색 원이 큰 것이 있고 작은 것이 있는 이유는 weight가 크고 작음을 의마한다. typedef struct _particle_t {    int x;    int y;    float weight; }particle_t; 노이즈가 있는 환경에서 레이저 센서나 레이더등을 이용하여 물체의 위치를 측정할 때를 가정해 순서대로 이 필터가 동작하는 과정에 대해 알아 보겠다. (순서는 위 그림과 조금 다르다) 1. 초기 상태 측정 범위안에 랜덤 혹은 일정한 간경으로 particle을 뿌려 놓는다. 이 떄 particle의 weight는 모두 0으로 한다.   for(i
자세한 내용 보기

바로 프로젝트 적용 가능한 FIR Filter (low/high/band pass filter )를 c나 python으로 만들기

이미지
 band pass filter 코드가 필요해 인터넷을 검색을 해봤지만, 쓸 만한 코드와 자료가 별로 없어 그냥 만들어 쓰기로 했다. 필요한 공식은 위키피디아에 나와있다. https://en.wikipedia.org/wiki/Finite_impulse_response 1. filter 크기 FIR filter를 만들 때 filter의 크기(tpas의 크기)를 결정할 때 도움될 수식이 있다.   N=A㏈Fs / 22*Δf A㏈ : 감쇄될 db Fs : Sampling Frequency Δf : frequency bandwidth 출처 :  https://www.allaboutcircuits.com/technical-articles/design-of-fir-filters-design-octave-matlab/ def estimatefilterlen(fl,fh,fs,db):     N = int(np.round(db * fs / (22 * (fh-fl)))-1)     return N 2. Low Pass Filter [ python code ] def lowpassfilter(f,fs,N):     taps = np.zeros(N)          fc = f / fs     omega = 2*np.pi*fc     middle = int(N/2)     for i in range(N):         if i == middle:             taps[i] = 2*fc         else:             taps[i] = np.sin(omega*(i-middle))/(np.pi*(i-middle))              return taps [ c code ] #ifndef PI #define PI 3.14159265358979323846 #endif
자세한 내용 보기

Eigen 라이브러리 matrix 사용법

이미지
본 글은 eigen 라이브러리의 matrix 간단한 사용법에 관한 글을 싣고 있다. Eigen 라이브러리는 http://eigen.tuxfamily.org/ 에서 다운받아 사용할 수 있다. 별도의 라이브러리 빌드 없이 헤더만 불러와 사용할 수 있어 편리하다. 개인적으론 opencv의 Mat보다 Eigen의 Matrix가 사용하기 편한 것 같다. 라이브러리 경로 project 설정 추가 다운받은 라이브러리 파일은 압축해제 후 project 설정에서 경로만 추가하면 된다. 라이브러리 헤더  #include <Eigen/Dense> namespace  아래와 같이 Eigen namespace를 사용하도록 설정하던가, 아니면, 호출 때 마다 Eigen::을 사용하면 된다.  using namespace Eigen; matrix 데이터 형식  Matrix 타입의 맨 뒤는 matrix내 데이터의 형식을 의미한다. EIGEN_MAKE_TYPEDEFS_ALL_SIZES(int,                  i) EIGEN_MAKE_TYPEDEFS_ALL_SIZES(float,                f) EIGEN_MAKE_TYPEDEFS_ALL_SIZES(double,               d) EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<float>,  cf) EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<double>, cd) matrix 타입 Eigen의 matrix 타입은 고정 사이즈와 가변 사이즈 두가지가 있다. 아래는 두 타입을 비교한 코드다. Dynamic-Size의 matrix타입은 사용시 matrix의 row수와 column수를 정해줘야 한다.  Fixed-Size  Dynamic-Size  Eigen::Matrix3d m; m << 1, 2, 3, 4, 5, 6, 7, 8, 9; m = m + Eigen::Matrix3d::Constant
자세한 내용 보기

base64 인코딩 디코딩 예제 c 소스

이미지
본 글은 base64 encoder/decoder예제 c 소스 코드를 싣고 있다. BASE64 BASE64는 일련의 데이터를 64종류" ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/ "의 문자만을 이용하여 표현하도록 하는 부호화 방식이다. Base64는 인터넷 이메일 프로토콜이나 http url등에서 사용된다.  부호화 과정은 아래와 같다. 문자 : Hello Hex 코드 : 0x48 0x65 0x6c 0x6c 0x6f 8bits bin 코드: 01001000 01100101 01101100 01101100 01101111 8비트 단위의 문자를 6비트 단위로 분리한다. 6bits bin 코드 : 010010  000110  010101  101100  011011  000110  1111  00  ‘=’ 6비트 2진 코드 값은 BASE64 코드 테이블의 인덱스이고 각 인덱스의 문자로 base64 부호화 한다. 마지막 모자라는 비트는 0으로 채우며, 4글자 단위로 부호화 하기 때문에 4글자가 안되면 ‘=’기호를 사용해 4글자를 맞춘다. Base64 코드 : aGVsbG8= Base64 Encoder source code Base64 인코더 소스는 아래와 같다. 위의 부호화 예대로 text데이터를 받아 6bit 단위로 나누어 base64 코드 테이블에 맞춰 부호화를 진행한다. static const char base64code_ascii[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"; int text2base64_size(const char *text) { int len = strlen(te
자세한 내용 보기

python winsound를 이용한 윈도우 환경에서 소리 재생 예제

 본 포스트에는 윈도우 환경에서 python의 winsound라이브러리를 이용한 사운드 재생의 예제를 싣고 있다. import winsound 특정 주파수 톤 재생 Beep를 사용하여 특정 주파수의 소리를 정해진 시간동안 재생할 수 있다. winsound.Beep(frequency, duration) 파라미터 frequency : Hz 단위며, 범위는 37Hz~32,767Hz duration : milliseconds 단위이다. example code import winsoind winsound.Beep(370,2000) wave 재생 지정된 파일, 윈도우 레지스트리에 등록된 사운드 혹은 오디오 데이터(bytes-like object)를 재생한다. 만약 sound 리소스가 없는 경우 시스템 디폴트 사운드를 재생한다.  winsound.PlaySound(sound, flags) 파라미터 sound : 파일 이름, 시스템 사운드 별명, 오디오 데이터 flags : winsound.SND_FILENAME : sound 파라미터는 wave file 이름 이며, SND_ALIAS와 함께 사용하면 안된다.  winsound.SND_ALIAS : sound 파라미터는 시스템 사운드 별명이며, SND_FILENAME와 함께 사용하면 안된다. ex) 'SystemAsterisk', 'SystemExclamation', 'SystemExit', 'SystemHand', 'SystemQuestion' winsound.SND_LOOP : 반복 재생한다. 테스트해보니 SND_ASY
자세한 내용 보기