用连续小波变换估计功率谱密度

Question

我想用连续小波变换和Morlet小波估计功率谱密度。贝娄，你可以找到我正在使用的函数。对于以下代码是否正确有任何评论或建议？

import pycwt as wavelet


mother_wave_dict = {
    'gaussian': wavelet.DOG(),
    'paul': wavelet.Paul(),
    'mexican_hat': wavelet.MexicanHat()
}

def trace_PSD_wavelet(x, dt, dj,  mother_wave='morlet'):
    """
    Method to calculate the  power spectral density using wavelet method.
    Parameters
    ----------
    x : array-like
        the components of the field to apply wavelet tranform
    dt: int
        the sampling time of the timeseries
    dj: determines how many scales are used to estimate wavelet coeff
    
        (e.g., for dj=1 -> 2**numb_scales 
    mother_wave: str
 
    Returns
    -------
    db_x,db_y,db_zz: array-like
        component coeficients of th wavelet tranform
    freq : list
        Frequency of the corresponding psd points.
    psd : list
        Power Spectral Density of the signal.
    psd : list
        The scales at which wavelet was estimated
    """
    

    if mother_wave in mother_wave_dict.keys():
        mother_morlet = mother_wave_dict[mother_wave]
    else:
        mother_morlet = wavelet.Morlet()
        
    N                       = len(x)

    db_x, _, freqs, _, _, _ = wavelet.cwt(x, dt,  dj, wavelet=mother_morlet)

     
    # Estimate trace powerspectral density
    PSD = (np.nanmean(np.abs(db_x)**2, axis=1))*( 2*dt)
    
    # Also estimate the scales to use later
    scales = ((1/freqs)/dt)#.astype(int)
    
    return db_x, freqs, PSD, scales

Answer 1

如果不知道你所说的“正确”是什么意思，就很难回答这个问题。

据我所知，您的代码允许：

选择一个给定的wavelet

• Compute，所选wavelet

• Compute的输入数据的cwt，整个记录的CWT转换数据的
• 。

我能够对scipy库中的一个电图数据示例运行您的代码，它按预期运行：

from scipy.misc import electrocardiogram
ecg = electrocardiogram()

由于这些数据是360 at采样的，所以我使用dt=1/360：

db_x, freqs, PSD, scales = trace_PSD_wavelet(ecg, 1/360, 1/24, 'morlet')

绘制输出db_x

fig = plt.imshow(np.abs(db_x), extent=[db_x.shape[1],db_x.shape[0],scales[-1],scales[]], aspect='auto')
plt.xlabel('time')
plt.ylabel('scale')

​

​

绘制相应的"PSD"：

​

​

你所谓的"PSD“测量的是CWT中包含的能量--在每个尺度上转换的数据，在整个记录中的平均值，这个例子中的5分钟数据。我不知道你计划如何使用这些信息，但请注意，这不是原始输入时域数据的PSD。

最后，关于Python实现，您可以简化调用默认小波的方式。只需将Morlet小波添加到字典中：

mother_wave_dict = {
    'gaussian': wavelet.DOG(),
    'paul': wavelet.Paul(),
    'mexican_hat': wavelet.MexicanHat()
    'morlet': wavelet.Morlet()
}

然后，您可以避免函数中的if语句，只需调用：

mother_morlet = mother_wave_dict[mother_wave]。