functions.py 文件源码

python
阅读 24 收藏 0 点赞 0 评论 0

项目:decode 作者: deshima-dev 项目源码 文件源码
def psd(data, dt, ndivide=1, window=hanning, overlap_half=False):
    """Calculate power spectrum density of data.

    Args:
        data (np.ndarray): Input data.
        dt (float): Time between each data.
        ndivide (int): Do averaging (split data into ndivide, get psd of each, and average them).
        ax (matplotlib.axes): Axis you want to plot on.
        doplot (bool): Plot how averaging works.
        overlap_half (bool): Split data to half-overlapped regions.

    Returns:
        vk (np.ndarray): Frequency.
        psd (np.ndarray): PSD
    """
    logger = getLogger('decode.utils.ndarray.psd')

    if overlap_half:
        step = int(len(data) / (ndivide + 1))
        size = step * 2
    else:
        step = int(len(data) / ndivide)
        size = step

    if bin(len(data)).count('1') != 1:
        logger.warning('warning: length of data is not power of 2: {}'.format(len(data)))
    size = int(len(data) / ndivide)
    if bin(size).count('1') != 1.:
        if overlap_half:
            logger.warning('warning: ((length of data) / (ndivide+1)) * 2 is not power of 2: {}'.format(size))
        else:
            logger.warning('warning: (length of data) / ndivide is not power of 2: {}'.format(size))
    psd = np.zeros(size)
    T   = (size - 1) * dt
    vs  = 1 / dt
    vk_ = fftfreq(size, dt)
    vk  = vk_[np.where(vk_ >= 0)]

    for i in range(ndivide):
        d = data[i * step:i * step + size]
        if window is None:
            w    = np.ones(size)
            corr = 1.0
        else:
            w    = window(size)
            corr = np.mean(w**2)
        psd = psd + 2 * (np.abs(fft(d * w)))**2 / size * dt / corr

    return vk, psd[:len(vk)] / ndivide
评论列表
文章目录


问题


面经


文章

微信
公众号

扫码关注公众号