Python频谱分析 [英] Python Spectrum Analysis
问题描述
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的输出进行比较时，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方法起作用.
此处的问题是您没有正确处理信号采样.在欢迎会议中，您将考虑采样频率为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)
您的信号因此没有得到定期采样.因此，您需要考虑这一点，否则您将无法正确估计PSD，这将给您带来错误的估计.
第一个更正应该是正确使用welsch中的参数fs.此参数指示给定信号的采样频率.设为4表示您的时间向量应为规则的[0, .25, .5, .75, .1, ....].更好的估计值是dt或len(t)/(t.max()-t.min())的中位数，所以在4/3左右.
这样可以为某些常数提供更好的PSD估计和正确的阶数，但与Kubios值相比仍然有所不同.
要正确估计PSD，您应该使用非均匀DFT . 可以在此处找到实现这种转换的软件包.该软件包的文档非常含糊，但是您需要使用adjoint方法来获得傅立叶变换而不会出现缩放问题:
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估计.您可以尝试将焊接技术与该nfft结合使用，方法是将窗口信号的估计值平均为它不依赖于FFT，而是依赖于PSD的任何估计.
I am trying to estimate the PSD of the heart rate variability of an ECG signal. To test my code,I have extracted the R-R interval from from the fantasia ECG database. I have extracted the signal can be accessed here. To calculate the PSD, I am using the welch method as shown below:
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)
Next,the area under the curve is calculated to estimate the power spectrum of the different HRV bands as shown below
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)
I then plot the PSD and get the following plot
At first I tough the output looks okay. However, when I compare my output with that of Kubios, which is a software than analyze HRV, I noticed that there are differences. The following chart shows the expected value for the PSD as calculated by Kubios Namely, the two plots are visually different and their values are way different. To confirm this, a print out of my data clearly shows that my calculation are wrong
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
I am thus, wondering:
- Can someone suggest a document I should read to improve my understanding on spectra analysis?
- What's wrong with my approach?
- How do I choose the most suitable parameters for the welch function?
- While the two plots somehow have the same shape, the data is completely different. How can I improve this?
- Is there a better approach to solve this? I am thinking about using the the Lomb-Scargle estimate but I am waiting to get at least the Welch method to work.
The problem here is that you do not handle correctly the sampeling of your signal. In your welsch call, you consider a regularly sampled signal with sample frequency 4Hz. If you look at the time vector 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)
Your signal is thus not regularly sampled. You thus need to take that into acount, else you do not estimate correctly the PSD and this gives you bad estimates.
A first correction should be to use correctly the parameter fs in welsch. This parameter indicates the sampling frequency of the given signal. Putting it ot 4 means that your time vector should be a regular [0, .25, .5, .75, .1, ....]. A better estimate would be either the median of dt or len(t)/(t.max()-t.min()), so around 4/3.
This gives better PSD estimate and correct order for some of the constant but it is still different compared to Kubios values.
To get correct estimate of the PSD, you should use a non uniform DFT. A package that implement such transform can be found here. The documentation is quite cryptic for this package but you need to use the adjoint method to get the Fourier Transform without scaling issue:
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))
Here the estimates do not seems to be in line with the one from Kubios. It is possible the estimation is probably of because you do a psd estimate on the whole signal. You can try to use the welch technique combined with this nfft by averaging estimates on windowed signals as it do not rely on FFT but on any estimation of the PSD.
这篇关于Python频谱分析的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
