trouble with r peak detection using discrete wavelt transform

31 views Asked by At

I got peak indexes on detailed coefficient5. However, I need to get indexes on signal. How do I need to transform index?

The individual recordings are each 10 hours in duration, and contain two ECG signals each sampled at 250 samples per second with 12-bit resolution over a range of ±10 millivolts. ambulatory ECG recorders with a typical recording bandwidth of approximately 0.1 Hz to 40 Hz.

import numpy as np
import pywt
import wfdb
import os
import matplotlib.pyplot as plt

def soft_threshold(signal,std):
    
    # threshold =  std *np.sqrt(2 * np.log(len(signal)))  #OK  
    threshold =0.8*max(signal)
    
    return np.sign(signal) * np.maximum(np.abs(signal) - threshold, 0)

def denoising(signal):
    (cA10,cD10,cD9,cD8,cD7,cD6,cD5,cD4,cD3,cD2,cD1)= pywt.wavedec(signal, 'db10', level = 10,mode = 'symmetric') #분해해 bior1.3
    std =  np.median(np.abs(cD1)) / 0.6745
    # 기저선 변동 제거를 위해 cA10,cD10 없앰
    # 고주파 오염 노이즈원 및 전력간섭 제거를 위 cD1 ~ cD4 제거
    # np.zeros_like는 변수의 사이즈 만큼의 0로 가득찬 numpy array를 배출, np.zeros는 array 사이즈만 넣지만 like은 array 넣음 
    for coeffs in [ cA10 ,cD10,cD4,cD3,cD2,cD1]:
        #print(coeffs)
        coeffs  = np.zeros_like(coeffs)

    denoised_coeffs = [cA10]+  [cD10,cD9,cD8,cD7,cD6,cD5,cD4,cD3,cD2,cD1]

    denoised_signal = pywt.waverec(denoised_coeffs, 'db10') #디노이징 된 신호 재건하기
    
    # cD5 thresholding
    denoised_coeffs[6] = soft_threshold(denoised_coeffs[6],std)
    
    return denoised_signal,denoised_coeffs[6]

def detect_r_peaks(modified_cD5):
   

    # qrs 값 증폭하고 양수로 만들기 위해 제곱
    cD5_s = np.square(modified_cD5)
    
    # cD5를 10단계로 분해
    (A10,D10,D9,D8,D7,D6,D5,D4,D3,D2,D1) = pywt.wavedec(cD5_s, 'db10', level=10, mode = 'symmetric')
    
    # D1 ~ D5 제거
    for coeffs_d5 in [D5,D4,D3,D2,D1]:
        coeffs_d5 = np.zeros_like(coeffs_d5) #
    
    # reference window에서 R 피크 찾기
    Rf_window = pywt.waverec((A10,D10,D9,D8,D7,D6,D5,D4,D3,D2,D1), 'db10')

    # 임계값 설정
    T_window  = Rf_window.mean() /4 # /는 소수점 반환, //는 round down 하고 정수 반환
    
    # reference window에서 임계값보다 높은 값 가지면 r피크 영역
    # r_peaks = []

    # for j in range(1, len(Rf_window) - 1):
    #     if Rf_window[j] > T_window:
    #         r_peaks.append(j)

    ##
    r_peaks = []
    for j in range(1, len(Rf_window) - 1):
        if Rf_window[j] > T_window and Rf_window[j] > Rf_window[j-1] and Rf_window[j] > Rf_window[j+1]:
            r_peaks.append(j)
      
    return r_peaks,Rf_window

#sampling_rate = 250
ecg_20s = ecg[0: 5000]
#remove dc
ecg_non_dc = ecg_20s - np.mean(ecg_20s)
#denoising
ecg_denoised,cD5 = denoising(ecg_non_dc)
  
r_peaks,Rf_window = detect_r_peaks(cD5)
    
#plot
labels = ['denoised signal','D5 Coefficient after Soft Thresholding T=300','D5 Coefficient after Squaring and D1-D5 Band Cancellation','ECG Signal after R-peak Point Detection']
n =4
#for (i,j) in zip(range(n),[ecg_denoised,cD5,Rf_window,r_peaks]):
#    plt.subplot(n,1,i+1)
#    plt.plot(j) #5초
#    plt.title(labels[i])
#    plt.tight_layout()
    
plt.plot(Rf_window)
plt.plot(r_peaks[:99],Rf_window[r_peaks[:99]],'ro')

enter image description here

What I want to get is the exact point of R peaks.

0

There are 0 answers