MIT/BIH 부정맥 데이터베이스를 이용한 다중스케일 기반 피크검출 알고리즘의 검증

A assessment of multiscale-based peak detection algorithm using MIT/BIH Arrhythmia Database

Abstract

A robust new algorithm for R wave detection named for Multiscale-based Peak Detection(MSPD) is assessed in this paper using MIT/BIH Arrhythmia Database. MSPD algorithm is based on a matrix composed of local maximum and find R peaks using result of standard deviation in the matrix. Furthermore, By reducing needless procedure of proposed algorithm, improve algorithm ability to detect R peak efficiently. And algorithm performance is assessed according to detection rates about various arrhythmia database.

1. 서 론

생체신호처리 분야에서 신호의 피크를 검출하는 과정은 오프라인에서 후처리를 위해서도 매우 중요하므로 여러 가지 새로운 방법들이 끊임없이 개발되고 있다. 일반적으로 사용되고 있는 피크 검출 방법에는 대표적으로 역치값 설정 기술(Adaptive threshold)[3], 미분(Differentiation)[8], 웨이블렛 변환(Wavelet transform)[2], 힐버트 변환(Hilbert transform)[11] 등을 바탕으로 하는 알고리즘들이 있다. 하지만 이 알고리즘들은 신호를 분석하기 전에 상황에 따라서 많은 변수들(역치값[3] 또는 윈도우 길이[4])을 상황에 맞게 수정해야 하는 문제점을 가지고 있다. 알고리즘 적용단계에서 최적화를 위해 변수에 대한 수작업이 동반되는데, 잘 적용되는 경우와 잘 적용되지 않는 경우들이 발생하게 된다. 그리고 피크 검출률은 높은 수준의 결과들을 보이지만 계산양이 많거나, 고주파 잡음에 민감하거나, 요구되는 변수들이 많아 시스템으로 구현하기에는 어려움이 있는 알고리즘들이 대부분이다. 반대로 변수가 거의 필요 없는 알고리즘들은 심전도[5] 같은 비주기적인 스파이크 경향의 신호와 같은 특별한 경우일 때 낮은 검출률로 인해 사용에 제약이 따른다.

따라서 본 논문에서는 심전도 신호 중에서도 각종 심장질환을 진단할 때 중요한 척도가 될 수 있는 심전도 신호의 R-피크를 검출하기 위해, Felix Scholkmann[1] 등이 제안한 방법을 사용했다. 제안된 알고리즘은 기존 알고리즘들과는 달리 상황에 맞게 고려해야 할 변수들이 요구되지 않고, 비교적 간단하다는 장점을 가지고 있다. 추가적으로 심전도 신호에서의 R피크를 검출하는데 불필요한 계산과정을 줄임으로써 계산량을 줄이고 R피크의 특징이 잘 반영되도록 최적화시킨 후, 이의 성능을 비교 평가하고자 한다.

 

2. 본 론

2.1 심전도(Electrocardiogram)

심전도[5]는 심장의 전기적 활동을, 체표면에서 측정한 전기적 신호를 의미한다. 심전도의 한 주기(Cycle)는 각각 P, Q, R, S, T파(wave)로 구성되어 있다. P파는 심방의 탈분극, QRS파는 좌, 우심실의 탈분극, T파는 심실의 재분극에 의해 각각 발생한다. 일반적으로 심전도 신호는 호흡 및 전원 장치로 인하여 신호에 왜곡 현상이 나타난다. 심전도의 주파수대역은 0.05Hz에서 100Hz로 대부분 낮은 대역에 형성되고 있어서, 호흡과 움직임 등에 의한 저주파 잡음에 큰 영향을 받는다. 따라서 필터링을 통해 QRS파를 검출해 내는 것이 중요하다.

2.2 다중스케일 기반 R피크 검출 알고리즘

본 논문에서 제안하는 다중스케일 기반 피크검출 알고리즘은 Felix Scholkmann[1] 등에 의하여 처음으로 제안된 방법으로, 대상 신호에서 국지적 최대값 스칼로그램(scalogram)을 계산하여 피크 성분이 있는 곳의 위치를 확률적으로 파악한 뒤, 최종적으로 피크의 위치는 스칼로그램의 표준편차가 최소가 된다는 점을 이용해 피크를 검출하는 알고리즘이다.

그림 1다중스케일 기반 R피크 검출 알고리즘 흐름도 Fig. 1 Block diagram of proposed peak detection algorithm

2.2.1 호흡성분을 제거하기 위한 고역통과 필터

피실험자의 심전도 신호를 10초간 측정하고 샘플링 주파수 120Hz로 양자화 한 뒤, 심전도의 한 주기가 4개가 포함되어 있도록 약 4초간의 심전도 신호를 추출하였다.

그림 2피실험자의 심전도 신호 Fig. 2 ECG signal

2.2.2 국지적 최대값 스케일러그램(Local Maximum Scalogram)

양의 방향으로 솟아오른 피크의 특징을 이용하기 위하여, 총 길이가 N개인 양자화 된 심전도 신호를 x라 했을 때, xi를 기준으로 양옆 샘플간격이 k개인 심전도 신호 3개 {xi−k, xi, xi+k}를 식 (1)에 따라 비교해 국지적 최대값 mk,i을 계산한다. 비교한 세 값 중 xi가 가장 큰 값이라면 성분 값을 0이라 하고, 그렇지 않으면 임의의 값 r + α로하여 LMS 행렬을 계산한다. 여기서, r은 0과 1 사이의 임의의 값이고, α는 상수 1이다. 예를 들어 i = k+1, ..., N−k에 대해, 행렬의 첫 번째 행(k = 1)은 xi를 기준으로 xi−1과 xi+1을 비교하고, 두 번째 행(k = 2)은 xi를 기준으로 xi−2와 xi+2를 비교하는 방법으로 N/2 번째 행까지 만들어진다. 이렇게 만들어진 행렬을 ‘LMS 행렬’ 이라고 하며 식 (2)처럼 정의한다.

그림 3은 k간격으로 이격된 샘플을 이용하여 계산된 LMS 행렬을 나타낸 그래프이다. 검은색으로 나타난 성분들이 값이 0, 즉 국지적 최대값이다. 따라서 LMS 행렬은 국지적 최대값의 분포에 대한 정보를 담고 있다

그림 3k간격으로 이격된 샘플을 이용한 LMS 행렬 Fig. 3 LMS Matrix using sample of distant k intervals

2.2.3 Gamma 배열과 Lambda 인덱스의 설정

다음 단계는 2.2.2절에서 구한 LMS 행렬에 대하여, 다음의 식 (3)와 같이, 각각의 행(row)을 기준으로 합계를 구한다. 그 결과 만들어진 배열을 ‘Gamma 배열’이 라고 정의한다. Gamma는 국지적 최대값의 분포에 대한 정보를 담고 있는 중요한 배열이다.

식 (1)에 의해, LMS 행렬에서 1차적으로 피크성분이라고 추정되는 부분은 그 성분값들이 0이므로, 식 (3)에 의하여 계산되는 각 행의 총합은 상대적으로 작은 값을 가지는 중요한 성질을 확인 할 수 있다(그림 4).

그림 4LMS배열에서 각각의 행을 기준으로 한 Gamma배열 Fig. 4 Gamma array based on row criteria in LMS matrix

Gamma 배열의 값 중 최소값을 가질 때의 인덱스를 ‘Lambda(λ) 인덱스’ 라고 정의한다. 그림 3에서 볼 수 있듯이, λ = 60인 경우, LMS 행렬의 어떤 한 데이터(열)가 주기적으로 반복되는 피크들 중 하나라면 60번째 행까지 그 열의 성분값 mk,i가 0이라는 것을 뜻함과 동시에, 이 데이터의 양옆 샘플간격이 60개 내에는 자신보다 큰 값이 없다는 의미가 되므로 Lambda 인덱스는 샘플간격과 연관되어 있음을 알 수 있다.

2.2.4 축소된 LMS 행렬(Reduced LMS matrix)

앞에서 구한 LMS 행렬은 행의 범위가 {1, 2, ⋯, N/2}, 열의 범위가 {1, 2, ⋯, N}에 이른다. 여기서 2.2.3절에서 구한 Lambda 인덱스를 적용시켜 행의 범위를 {1, 2, ⋯, λ}로 축소할 수 있다. 이 축소된 행렬을 “축소된 LMS 행렬(RLMS matrix)” 이라고 정의한다.

그림 4의 왼쪽 그림은 축소되기 이전 LMS 행렬, 오른쪽 그림은 축소된 이후의 LMS 행렬을 나타낸다. LMS 행렬을 축소하게 되면 오른쪽 그림의 빨간색 네모박스와 같은 특징이 드러난다. 각 열의 요소값이 모두 0을 나타내고 있다. 이 특징을 이용해 다음 절에서 피크의 위치를 검출해낼 수 있다.

그림 5기존 LMS 행렬(좌),축소된 LMS 행렬(우) Fig. 5 LMS Matrix(L), Rescaled LMS Matrix(R)

2.2.5 RLMS 행렬 열(column) 표준편차의 계산과 R피크의 위치

다음 식 (5)와 같이 축소된 LMS 행렬의 열(column)을 기준으로 표준편차를 구함으로써 R피크의 위치를 추정할 수 있다.

예를 들어, 100번 째 열의 데이터가 찾고자 하는 심전도의 R피크일 때, LMS 행렬에서 1행 ～ λ행까지 모두 mk,100 = 0이므로 표준편차 σ100 또한 0이다. 그 결과 100번째 열의 데이터를 R피크라고 인식하게 된다(그림 6).

그림 6축소된 LMS 행렬에서 열(column) 기준의 표준편차 Fig. 6 Standard deviation of Rescaled LMS Matrix

그림 7은 축소된 LMS 행렬에서 열(column)을 기준으로 계산한 표준편차 σi 중에 σi = 0 일 때의 i값만 배열로 구성하여 ECG 원본 신호와 매칭 시킨 결과를 나타낸다.

그림 7표준편차 σi = 0와 신전도 신호 R피크와의 관계 Fig. 7 Detection result of R peak

정확한 R피크를 검출하기 위해서는 Lambda 인덱스 값이 어느 범위에 있는지가 중요한 요소로 작용한다. 만약 정확한 R피크를 검출하기 위한 Lambda 인덱스 값이 60인데, 60이하의 값이 나오면 실제 피크인 위치에서 표준편차가 0이 나오지 않아 R피크를 놓치는 경우가 발생한다. 반대로 Lambda 인덱스 값이 60이상이면 R피크는 물론이고 P파나 T파의 위치까지도 표준편차가 0이 나오게 되므로 원하지 않은 피크까지 검출되는 문제가 발생할 수 있는 가능성이 있다.

2.3 알고리즘 최적화

이런 과정을 통해 생체 신호의 피크 위치를 검출해낼 수 있다. 그러나 일반적인 신호에서의 피크성분을 검출하기 위하여 버퍼로 사용하는 메모리 공간이 과다하게 많으며, 심전도 신호의 경우, 심전도 신호의 특성이 반영이 되지 않는 경우가 많으므로 알고리즘을 심전도에 그대로 적용하기에는 무리가 있다. 따라서 심전도 신호에서의 R피크 검출을 위해서는 계산량을 줄이고 심전도의 특성이 잘 반영되기 위한 최적화가 필요하다.

2.3.1 P파, T파를 제외시키기 위한 Gamma 범위 축소

심전도의 QRS파뿐만 아니라 P파와 T파 또한 주기적인 성향을 가지고 있으므로, 제안한 MSPD 알고리즘은 P와 T파 또한 피크로 인식할 가능성이 있다. 따라서 P파와 T파가 검출되지 않는 범위로 mk를 수정할 필요가 있다. 이를 위하여, 실험적으로 mk를 다음의 식 (6)과 같이 제한하였다. mk의 범위를 축소시킨 후에 Gamma 배열과 Lambda 인덱스를 찾는 과정을 진행한다.

그림 8Gamma에 나타난 심전도특성 Fig. 8 ECG Characteristics in Gamma

정상적인 심전도의 경우 P파의 특성은 mk = 40, T파의 특성이 mk = 75 부근에서 나타난 모습을 볼 수 있다. MIT/BIH Arrhythmia Database의 심전도 데이터는 fs = 360Hz 이므로, mk의 범위는 102 ～ 240로 축소된다. 따라서 제안한 알고리즘은 일반적인 R피크일 경우, P와 T 파는 검출되지 않고 R피크만 검출한다.

2.3.2 구간 overlapping

그림 9Overlapping 흐름도 Fig. 9 Overlapping Block diagram

또한 제안한 알고리즘은 메모리 사용의 효율성을 향상시키기 위하여, 전체 심전도 데이터를 한 번에 처리하지 않고, window의 크기를 1000개로 고정시켜 심전도 샘플구간을 나눠서 피크를 검출했다. 전체 데이터가 N개 일 때, 구간을 1～1000, 501～1500, 1001～2000, ..., (N-999)～N으로 나누고, 겹치는 구간(501～1000)구간은 앞 구간(1～1000)에서 검출된 피크와 뒷 구간(501～1500)에서 검출된 피크를 OR연산을 통해 결정해준다. 그림 10은 데이터 샘플의 중첩과정의 유무에 따른 R피크 검출 결과를 나타낸다.

그림 10Overlap 하지 않았을 때(위), Overlap 했을 때(아래) Fig. 10 No overlap(up), overlap(down)

데이터의 전체길이가 길어 한 번에 알고리즘을 처리하기에는 다소 무리가 있으므로 구간을 나누어 알고리즘에 적용시킨다. 그리고 이에 따라 발생하는 문제는 원도우 크기에 따라 구간이 나눠지는 경계에 근접해 있는 R피크가 검출되지 않은 결과를 그림 10의 위 그림에서 볼 수 있으나, 중첩 처리를 적용함으로써 R피크 검출의 정확도를 높일 수 있다

2.3.3 QRS파 검출을 위한 대역통과 필터

QRS파의 주파수 대역만 통과시키기 위해 Butterworth 8～20Hz 대역의 대역통과 1차 필터를 적용시켰다. QRS파를 검출하기 위해 제안된 차단주파수 대역은 다음과 같다. Thakor[6]와 Chen[7] 5～15Hz, Pan and Tompkins[8] 5～11Hz, Cuiwei[9] 8～58.8Hz, Sahambi[10] 3～40Hz, Benitez[11] 8～20Hz, Moraes[12] 9～30Hz, Mahmoodabadi[13] 2～40Hz. 논문에 제시된 대역을 알고리즘에 적용시켜본 결과 8～20Hz의 경우가 가장 효율적이었다. 결과는 표 4에서 확인할 수 있다.

표 4BPF 주파수 대역별 예측도와 민감도 결과 Table 4 +P and SE results on BPF frequency ranges

 

3. 결 과

심전도의 R피크 검출 성능을 평가하기 위하여 MIT/BIH Arrhythmia Database에서 제공하는 약 1시간 길이의 48개 2채널 심전도 신호 중, 첫 번째 채널을 사용하였다. 피험자 1명당 데이터 개수는 650000개이다. MIT/BIH Arrhythmia Database는 Annotation R피크 정보 또한 제공해 주기 때문에, 알고리즘이 찾아낸 R피크가 정확한지 아닌지를 확인할 수 있다. 피크검출 알고리즘의 성능을 판단하기 위해 검출 결과를 표 1에 제시된 조건(Condition)과 검출결과(Test)에 따라 True Positive, False Positive, False Negative, True Negative로 나눈 후, 그 결과를 표 2에 제시된 평가기준에 적용시킨다. 이 결과를 통해 알고리즘이 얼마나 효과적으로 R피크를 검출했는지 알 수 있다. TP(True Positive)는 Annotation에 제시된 샘플번호를 기준으로 샘플 개수 40개 이내로 피크가 검출되었을 경우 TP라고 설정하였다. 마찬가지로 FP, FN, TN 또한 표 2의 정의에 따라 검출결과를 설정했다.

표 1분할표[15] Table 1 Contingency table

표 2분할표를 이용한 평가기준 Table 2 Appraisal standard using contingency table

표 2를 R피크 검출 알고리즘에 적용해 설명하면, 예측도는 알고리즘에 의해서 피크라고 검출된 것 중 실제 피크가 얼마나 들어있는지에 대한 비율을 의미한다. 민감도는 실제 피크 중 알고리즘에 의해 몇 개나 검출되었는지에 대한 비율이다. 특이성은 피크가 아닌 것 중 알고리즘에 의해 피크가 아니라고 검출된 것의 비율이다. 정확도는 전체 샘플에 대해서 알고리즘이 정확하게 피크인 것과 아닌 것을 구분했는지에 대한 확률이다.

표 3은 본 논문에서 제안한 방법으로 피크 검출 알고리즘을 최적화 시킨 후, MIT/BIH Arrhythmia Database에 적용하여 얻어진 R피크 검출 결과를 나타낸다.

표 3MIT/BIH arrhythmia database에서의 R피크 검출 결과 Table 3 Detection results on MIT/BIH arrhythmia database

표 4는 주파수 대역별 예측도와 민감도 결과만 정리한 것이다. 이를 보아 알 수 있는 것은 대역통과 필터를 사용하지 않은 경우에 예측도(+P) 88.82%와 민감도(SE) 94.32%로, 가장 낮은 검출율을 보였다. 그리고 필터를 사용한 경우, 8～20Hz의 대역을 가진 필터에서 예측도(+P) 96.17%와 민감도(SE) 99.54%로 가장 높은 검출 결과를 얻을 수 있었다.

그러나 본 논문에서 제안한 알고리즘이 모든 경우에 대하여, 최적화된 결과를 보여주는 것은 아니다. 다음의 그림 11은 정확한 피크검출이 어려운 심전도 신호의 예를 보여주고 있다.

그림 11정확한 R피크를 검출하기 어려운 경우 Fig. 11 Bad case to detect R peak

일반적으로 심전도 신호가 구분이 명확한 P-QRS-T 파의 형태를 가지고 있으면 검출 정확도가 매우 높게 나타나나, 그림 11처럼 비정상적인 파형을 가지고 있을 때는 검출율이 상대적으로 낮게 나오는 경향이 있다.

기존의 문헌에서 제안된 대표적인 R피크 검출 알고리즘과 본 논문에서 제안한 알고리즘의 R피크 검출결과를 비교해 보았을 때, 예측도(+P)는 약 3%정도 차이가 났으며 민감도(SE)는 큰 차이가 없는 것을 확인 할 수 있었다.

표 5다른 R피크 검출 알고리즘의 +P와 SE Table 5 Contingency table

 

4. 결 론

본 논문에서는 국지적 최대값 스케일러그램을 구해 피크의 위치를 검출해내는 새로운 R피크 검출 방법을 MIT/BIH 부정맥 데이터에 적용시켜 보았다. 그 결과, 부정맥 같은 불규칙한 피크가 갑자기 나타나거나, 베이스라인이 흔들려도 유사한 크기의 피크만 주기적으로 나온다면 무리 없이 검출할 수 있다. 차단주파수가 0.1Hz인 고역통과 필터를 추가해 호흡에 의한 잡음을 없애주고, QRS파 검출을 위한 8～20Hz의 대역통과 필터를 추가함으로써 알고리즘을 최적화시킬 수 있다. 실제로 부정맥 데이터를 적용시켜 보았을 때, +P와 SE 가 각각 96.17%, 99.53% 로 최적화과정을 거치기전보다 더 정확한 피크 검출이 가능했다.

R피크를 찾는 다른 기존의 방법들과 비교해보았을 때, 제시된 알고리즘은 다음과 같은 특징을 가진다. 첫 번째, 신호를 분석하기 전에 사용자에 의해 선택되어져야 하는 변수가 필요 없다. 두 번째, 주기적인 성향을 가진 신호에서 정확한 피크를 검출할 수 있다. 이러한 이유로 제시된 알고리즘은 R피크 검출에 대해 고효율 알고리즘이다.

제안한 알고리즘은 심전도뿐만 아니라, 다양한 전기적 생체 신호에 위와 같은 알고리즘을 적용시켜, 효과적인 검출을할 수 있을 것이다. 그리고 알고리즘을 더욱 단축시킨다면 실시간 R피크 검출이 가능할 것이다. 이렇게 다양한 생체신호를 수집하여, 정확한 검출이 가능한 자동화 시스템을 개발한다면 u-health 산업을 크게 발전시킬 수 있을 것이다.