파이썬 스펙트럼 분석

Question

ECG 신호의 심박수 가변성의 PSD를 추정하려고합니다. 코드를 테스트하기 위해 fantasia ECG database에서 R-R 간격을 추출했습니다. 신호를 추출한 후 here에 액세스 할 수 있습니다. 아래에 도시 된 바와 같이, I는 웰치 방식을 사용하고있는 PSD를 계산하기 :

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.signal import welch 
ibi_signal = np.loadtxt('fantasia-f1y01-RR.txt') 
t = np.array(ibi_signal[:, 0]) # time index in seconds 
ibi = np.array(ibi_signal[:, 1]) # the IBI in seconds 
# Convert the IBI in milliseconds 
ibi = ibi * 1000 
# Calculate the welch estimate 
Fxx, Pxx = welch(ibi, fs=4.0, window='hanning', nperseg=256, noverlap=128) 

ulf = 0.003 
vlf = 0.04 
lf = 0.15 
hf = 0.4 
Fs = 250 
# find the indexes corresponding to the VLF, LF, and HF bands 
ulf_freq_band = (Fxx <= ulf) 
vlf_freq_band = (Fxx >= ulf) & (Fxx <= vlf) 
lf_freq_band = (Fxx >= vlf) & (Fxx <= lf) 
hf_freq_band = (Fxx >= lf) & (Fxx <= hf) 
tp_freq_band = (Fxx >= 0) & (Fxx <= hf) 
# Calculate the area under the given frequency band 
dy = 1.0/Fs 
ULF = np.trapz(y=abs(Pxx[ulf_freq_band]), x=None, dx=dy) 
VLF = np.trapz(y=abs(Pxx[vlf_freq_band]), x=None, dx=dy) 
LF = np.trapz(y=abs(Pxx[lf_freq_band]), x=None, dx=dy) 
HF = np.trapz(y=abs(Pxx[hf_freq_band]), x=None, dx=dy) 
TP = np.trapz(y=abs(Pxx[tp_freq_band]), x=None, dx=dy) 
LF_HF = float(LF)/HF 
HF_LF = float(HF)/LF 
HF_NU = float(HF)/(TP - VLF) 
LF_NU = float(LF)/(TP - VLF) 

아래에 도시 된 바와 같이, 곡선 아래의 면적은 다른 HRV 밴드의 파워 스펙트럼을 추정 산출

그 다음 PSD를 플롯하고 다음 플롯을 얻습니다.

처음에는 출력이 괜찮아 보입니다. 그러나 HRV를 분석하는 소프트웨어 인 Kubios의 결과와 내 결과를 비교할 때 차이가 있음을 발견했습니다. 다음 차트는 Kubios 에 의해 계산 된 PSD의 예상 값을 보여줍니다. 즉, 두 플롯은 시각적으로 다르고 그 값은 다릅니다. 이를 확인하려면 내 데이터 중 인쇄 명확하게 내 계산이 잘못 보여줍니다

ULF 0.0 
    VLF 13.7412277853 
    LF  45.3602063444 
    HF  147.371442221 
    TP  239.521363002 
    LF_HF 0.307795090152 
    HF_LF 3.2489147228 
    HF_NU 0.652721029154 
    LF_NU 0.200904328012 

궁금 따라서 이니

• 누군가가 내가 스펙트럼에 대한 나의 이해를 개선하기 위해 읽어야 할 문서를 제안 할 수 분석?
• 내 접근 방식에 문제가 있습니까?
• 웰치 기능에 가장 적합한 매개 변수는 어떻게 선택합니까?
• 두 플롯이 어떻게 든 같은 모양을 가지지 만 데이터는 완전히 다릅니다. 어떻게 개선 할 수 있습니까?
• 해결 방법이 더 있습니까? 나는 Lomb-Scargle 견적을 사용하려고 생각하고 있습니다 만 적어도 Welch 방법을 사용하려면 기다리고 있습니다.

Answer 1

여기에서 문제는 신호의 신호가 올바르게 처리되지 않는다는 것입니다. 귀하의 welsch 호출에서, 당신은 샘플 주파수 4Hz로 정기적으로 샘플링 된 신호를 고려하십시오. 시간 벡터를 보면 t

In [1]: dt = t[1:]-t[:-1] 

In [2]: dt.mean(), np.median(dt) 
Out[2]: 0.76693059125964014, 0.75600000000000023 

In [3]: dt.min(), dt.max() 
Out[3]: (0.61599999999998545, 1.0880000000000081) 

신호가 정기적으로 샘플링되지 않습니다. 따라서 이것을 acount로 가져 가야합니다. 그렇지 않으면 PSD를 정확하게 추정하지 못하고 나쁜 예측을하게됩니다.

첫 번째 수정은 welsch에서 fs 매개 변수를 올바르게 사용해야합니다. 이 파라미터는 주어진 신호의 샘플링 주파수를 나타냅니다. 당신의 시간 벡터가 보통 [0, .25, .5, .75, .1, ....]이어야한다는 것을 의미합니다. 더 나은 추정치는 dt 또는 len(t)/(t.max()-t.min())의 중간 값이므로 4/3 주위에 있습니다. 이것은 상수의 일부에 대해 더 나은 PSD 추정 및 올바른 순서를 제공하지만 Kubios 값과 비교하여 여전히 다릅니다.

정확한 PSD 추정치를 얻으려면 non uniform DFT을 사용해야합니다. 이러한 변환을 구현하는 패키지는 here입니다.문서는이 패키지에 대한 매우 애매하지만 당신은 푸리에 확장 문제없이 변환 얻기 위해 수반 행렬 방법을 사용해야합니다

N = 128 # Number of frequency you will get 
M = len(t) # Number of irregular samples you have 
plan = NFFT(N, M) 

# Give the sample times and precompute variable for 
# the NFFT algorithm. 
plan.x = t 
plan.precompute() 

# put your signal in `plan.f` and use the `plan.adjoint` 
# to compute the Fourier transform of your signal 
plan.f = ibi 
Fxx = plan.adjoint() 
plt.plot(abs(Fxx)) 

가 여기에 추정 Kubios에서 하나 줄 것으로 보인다 없습니다. 전체 신호에 대한 PSD 추정을하기 때문에 추정이 가능할 수 있습니다. FFT에 의존하지 않고 PSD를 추정 할 때 창 신호의 평균을 계산하여이 nfft와 결합 된 welch technique을 사용할 수 있습니다.