Python pywt.cwt方法代碼示例

# 需要導入模塊: import pywt [as 別名]
# 或者: from pywt import cwt [as 別名]
def MyPywtCWT(data):
    pass
    ''' --原始數據信息初始化--'''

    # Fs = 500000 #采樣頻率：500 000 Hz ; 采樣周期：2 us
    ''' --尺度計算-- '''
    wavename = 'gaus1'
    totalscal = 256          #尺度序列的長度
    
    #Fc = 2000;             #小波中心頻率(Hz)（“主波峰之間的差值”=2000Hz）
    #Fc = pywt.central_frequency(wavename, precision=8) 
    Fc = pywt.central_frequency(wavename)               
                   
    C = 2*Fc*totalscal      # C為常數,用於計算尺度序列. C = 2*Fc/totalscal
    scal= C/np.arange(1,totalscal+1)    #尺度序列,範圍（2*Fc,inf）

    #--連續小波變換--
    coef,freqs = pywt.cwt(data, scal, wavename)
    coef = np.abs(coef)
    return coef,freqs 

# 需要導入模塊: import pywt [as 別名]
# 或者: from pywt import cwt [as 別名]
def MyScipyCwt(data,MyWidths):
    ''' 將int型data轉為float型sig '''
    sig = np.ones(len(data),np.float)  #產生空的float型sig
    for i in range(0,len(data)): 
        sig[i] = float(data[i])

    # widths = np.arange(1, 31)
    widths = np.arange(1, MyWidths+1)
    ''' 
    signal.cwt(sig, signal.ricker, widths) 
        - CWT: Continuous wavelet transform 連續小波變換
        - signal.ricker返回一個Ricker小波，也被稱為“墨西哥帽子小波”   
    '''
    cwtmatr =  signal.cwt(sig, signal.ricker, widths) 
    # cwtmatr = np.abs(cwtmatr)
    # plt.imshow(cwtmatr, extent=[-1, 1, 31, 1], cmap='PRGn', aspect='auto', vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
    return cwtmatr 

# 需要導入模塊: import pywt [as 別名]
# 或者: from pywt import cwt [as 別名]
def child_wav(wavelet, scale):
    """Returns an array of complex values with the child wavelet used at the
    given ``scale``.

    The ``wavelet`` argument can be either a string like 'cmor1-1.5' or
    a ``pywt.ContinuousWavelet`` instance
    """

    wavelet = _wavelet_instance(wavelet)

    # the following code has been extracted from pywt.cwt() 1.0.2
    precision = 10
    int_psi, x = pywt.integrate_wavelet(wavelet, precision=precision)
    step = x[1] - x[0]
    j = np.floor(
            np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step))
    if np.max(j) >= np.size(int_psi):
                j = np.delete(j, np.where((j >= np.size(int_psi)))[0])

    return int_psi[j.astype(np.int)] 

# 需要導入模塊: import pywt [as 別名]
# 或者: from pywt import cwt [as 別名]
def _peaks_delineator(ecg, rpeaks, sampling_rate=1000):
    # Try loading pywt
    try:
        import pywt
    except ImportError:
        raise ImportError(
            "NeuroKit error: ecg_delineator(): the 'PyWavelets' module is required for this method to run. ",
            "Please install it first (`pip install PyWavelets`).",
        )
    # first derivative of the Gaissian signal
    scales = np.array([1, 2, 4, 8, 16])
    cwtmatr, __ = pywt.cwt(ecg, scales, "gaus1", sampling_period=1.0 / sampling_rate)

    qrs_duration = 0.1

    search_boundary = int(0.9 * qrs_duration * sampling_rate / 2)
    significant_peaks_groups = []
    for i in range(len(rpeaks) - 1):
        # search for T peaks and P peaks from R peaks
        start = rpeaks[i] + search_boundary
        end = rpeaks[i + 1] - search_boundary
        search_window = cwtmatr[4, start:end]
        height = 0.25 * np.sqrt(np.mean(np.square(search_window)))
        peaks_tp, heights_tp = scipy.signal.find_peaks(np.abs(search_window), height=height)
        peaks_tp = peaks_tp + rpeaks[i] + search_boundary
        # set threshold for heights of peaks to find significant peaks in wavelet
        threshold = 0.125 * max(search_window)
        significant_index = []
        significant_index = [j for j in range(len(peaks_tp)) if heights_tp["peak_heights"][j] > threshold]

        significant_peaks_tp = []
        for index in significant_index:
            significant_peaks_tp.append(peaks_tp[index])
        significant_peaks_groups.append(_find_tppeaks(ecg, significant_peaks_tp, sampling_rate=sampling_rate))

    tpeaks, ppeaks = zip(*[(g[0], g[-1]) for g in significant_peaks_groups])

    tpeaks = np.array(tpeaks, dtype="int")
    ppeaks = np.array(ppeaks, dtype="int")
    return tpeaks, ppeaks 

# 需要導入模塊: import pywt [as 別名]
# 或者: from pywt import cwt [as 別名]
def _find_tppeaks(ecg, keep_tp, sampling_rate=1000):
    # Try loading pywt
    try:
        import pywt
    except ImportError:
        raise ImportError(
            "NeuroKit error: ecg_delineator(): the 'PyWavelets' module is required for this method to run. ",
            "Please install it first (`pip install PyWavelets`).",
        )
    # first derivative of the Gaissian signal
    scales = np.array([1, 2, 4, 8, 16])
    cwtmatr, __ = pywt.cwt(ecg, scales, "gaus1", sampling_period=1.0 / sampling_rate)
    max_search_duration = 0.05
    tppeaks = []
    for index_cur, index_next in zip(keep_tp[:-1], keep_tp[1:]):
        # limit 1
        correct_sign = cwtmatr[4, :][index_cur] < 0 and cwtmatr[4, :][index_next] > 0  # pylint: disable=R1716
        #    near = (index_next - index_cur) < max_wv_peak_dist #limit 2
        #    if near and correct_sign:
        if correct_sign:
            index_zero_cr = signal_zerocrossings(cwtmatr[4, :][index_cur:index_next])[0] + index_cur
            nb_idx = int(max_search_duration * sampling_rate)
            index_max = np.argmax(ecg[index_zero_cr - nb_idx : index_zero_cr + nb_idx]) + (index_zero_cr - nb_idx)
            tppeaks.append(index_max)
    return tppeaks


# =============================================================================
#                              PEAK METHOD
# =============================================================================