我正在嘗試估計ECG信號的心率變異性的PSD。為了測試我的代碼，我從fantasia ECG database中提取了R-R間隔。我已經提取了信號，可以訪問here。為了計算PSD，我使用的是如下所示的welch方法：

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)

接下來，計算曲線下面積以估算不同HRV頻帶的功率譜，如下所示

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)

然后我繪制PSD并獲得以下圖表

起初我認為輸出看起來還不錯。但是，當我將我的輸出與Kubios(比分析HRV的軟件)進行比較時，我注意到存在差異。下圖顯示了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

因此，我想知道：

有人可以建議我應該閱讀的文件，以提高我對光譜分析的理解嗎？

我的方法出了什么問題？

如何為welch功能選擇最合適的參數？

雖然這兩個圖以某種方式具有相同的形狀，但數據完全不同。我該如何改進呢？

有沒有更好的方法來解決這個問題？我正在考慮使用Lomb-Scargle估計，但我等著至少讓Welch方法起作用。

總結

以上是生活随笔為你收集整理的python信号分析_Python频谱分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。