Python pywt.wavedec() Examples

def _wavelet_coefs(data, wavelet_name='db4'):
    """Compute Discrete Wavelet Transform coefficients.

    Parameters
    ----------
    data : ndarray, shape (n_channels, n_times)

    wavelet_name : str (default: db4)
         Wavelet name (to be used with ``pywt.Wavelet``). The full list of
         Wavelet names are given by: ``[name for family in pywt.families() for
         name in pywt.wavelist(family)]``.

    Returns
    -------
    coefs : list of ndarray
         Coefficients of a DWT (Discrete Wavelet Transform). ``coefs[0]`` is
         the array of approximation coefficient and ``coefs[1:]`` is the list
         of detail coefficients.
    """
    wavelet = pywt.Wavelet(wavelet_name)
    levdec = min(pywt.dwt_max_level(data.shape[-1], wavelet.dec_len), 6)
    coefs = pywt.wavedec(data, wavelet=wavelet, level=levdec)
    return coefs 

def apply(self, data):
        # data[ch][dim0]
        shape = data.shape
        out = np.empty((shape[0], 4 * (self.n * 2 + 1)), dtype=np.float64)

        def set_stats(outi, x, offset):
            outi[offset*4] = np.mean(x)
            outi[offset*4+1] = np.std(x)
            outi[offset*4+2] = np.min(x)
            outi[offset*4+3] = np.max(x)

        for i in range(len(data)):
            outi = out[i]
            new_data = pywt.wavedec(data[i], 'db%d' % self.n, level=self.n*2)
            for i, x in enumerate(new_data):
                set_stats(outi, x, i)

        return out 

def _assert_all_coeffs_equal(coefs1, coefs2):
    # return True only if all coefficients of SWT or DWT match over all levels
    if len(coefs1) != len(coefs2):
        return False
    for (c1, c2) in zip(coefs1, coefs2):
        if isinstance(c1, tuple):
            # for swt, swt2, dwt, dwt2, wavedec, wavedec2
            for a1, a2 in zip(c1, c2):
                assert_array_equal(a1, a2)
        elif isinstance(c1, dict):
            # for swtn, dwtn, wavedecn
            for k, v in c1.items():
                assert_array_equal(v, c2[k])
        else:
            return False
    return True 

def test_coeffs_to_array():
    # single element list returns the first element
    a_coeffs = [np.arange(8).reshape(2, 4), ]
    arr, arr_slices = pywt.coeffs_to_array(a_coeffs)
    assert_allclose(arr, a_coeffs[0])
    assert_allclose(arr, arr[arr_slices[0]])

    assert_raises(ValueError, pywt.coeffs_to_array, [])
    # invalid second element:  array as in wavedec, but not 1D
    assert_raises(ValueError, pywt.coeffs_to_array, [a_coeffs[0], ] * 2)
    # invalid second element:  tuple as in wavedec2, but not a 3-tuple
    assert_raises(ValueError, pywt.coeffs_to_array, [a_coeffs[0],
                                                     (a_coeffs[0], )])
    # coefficients as None is not supported
    assert_raises(ValueError, pywt.coeffs_to_array, [None, ])
    assert_raises(ValueError, pywt.coeffs_to_array, [a_coeffs,
                                                     (None, None, None)])

    # invalid type for second coefficient list element
    assert_raises(ValueError, pywt.coeffs_to_array, [a_coeffs, None])

    # use an invalid key name in the coef dictionary
    coeffs = [np.array([0]), dict(d=np.array([0]), c=np.array([0]))]
    assert_raises(ValueError, pywt.coeffs_to_array, coeffs) 

def test_wavedecn_coeff_reshape_even():
    # verify round trip is correct:
    #   wavedecn - >coeffs_to_array-> array_to_coeffs -> waverecn
    # This is done for wavedec{1, 2, n}
    rng = np.random.RandomState(1234)
    params = {'wavedec': {'d': 1, 'dec': pywt.wavedec, 'rec': pywt.waverec},
              'wavedec2': {'d': 2, 'dec': pywt.wavedec2, 'rec': pywt.waverec2},
              'wavedecn': {'d': 3, 'dec': pywt.wavedecn, 'rec': pywt.waverecn}}
    N = 28
    for f in params:
        x1 = rng.randn(*([N] * params[f]['d']))
        for mode in pywt.Modes.modes:
            for wave in wavelist:
                w = pywt.Wavelet(wave)
                maxlevel = pywt.dwt_max_level(np.min(x1.shape), w.dec_len)
                if maxlevel == 0:
                    continue

                coeffs = params[f]['dec'](x1, w, mode=mode)
                coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs)
                coeffs2 = pywt.array_to_coeffs(coeff_arr, coeff_slices,
                                               output_format=f)
                x1r = params[f]['rec'](coeffs2, w, mode=mode)

                assert_allclose(x1, x1r, rtol=1e-4, atol=1e-4) 

def test_wavedecn_coeff_reshape_axes_subset():
    # verify round trip is correct when only a subset of axes are transformed:
    #   wavedecn - >coeffs_to_array-> array_to_coeffs -> waverecn
    # This is done for wavedec{1, 2, n}
    rng = np.random.RandomState(1234)
    mode = 'symmetric'
    w = pywt.Wavelet('db2')
    N = 16
    ndim = 3
    for axes in [(-1, ), (0, ), (1, ), (0, 1), (1, 2), (0, 2), None]:
        x1 = rng.randn(*([N] * ndim))
        coeffs = pywt.wavedecn(x1, w, mode=mode, axes=axes)
        coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs, axes=axes)
        if axes is not None:
            # if axes is not None, it must be provided to coeffs_to_array
            assert_raises(ValueError, pywt.coeffs_to_array, coeffs)

        # mismatched axes size
        assert_raises(ValueError, pywt.coeffs_to_array, coeffs,
                      axes=(0, 1, 2, 3))
        assert_raises(ValueError, pywt.coeffs_to_array, coeffs,
                      axes=())

        coeffs2 = pywt.array_to_coeffs(coeff_arr, coeff_slices)
        x1r = pywt.waverecn(coeffs2, w, mode=mode, axes=axes)

        assert_allclose(x1, x1r, rtol=1e-4, atol=1e-4) 

def WTfilt_1d(sig):
    """
    # 使用小波变换对单导联ECG滤波
    # 参考：Martis R J, Acharya U R, Min L C. ECG beat classification using PCA, LDA, ICA and discrete
    wavelet transform[J].Biomedical Signal Processing and Control, 2013, 8(5): 437-448.
    :param sig: 1-D numpy Array，单导联ECG
    :return: 1-D numpy Array，滤波后信号
    """
    coeffs = pywt.wavedec(sig, 'db6', level=9)
    coeffs[-1] = np.zeros(len(coeffs[-1]))
    coeffs[-2] = np.zeros(len(coeffs[-2]))
    coeffs[0] = np.zeros(len(coeffs[0]))
    sig_filt = pywt.waverec(coeffs, 'db6')
    return sig_filt 

def getWaveletData(eda):
    '''
    This function computes the wavelet coefficients
    INPUT:
        data:           DataFrame, index is a list of timestamps at 8Hz, columns include EDA, filtered_eda
    OUTPUT:
        wave1Second:    DateFrame, index is a list of timestamps at 1Hz, columns include OneSecond_feature1, OneSecond_feature2, OneSecond_feature3
        waveHalfSecond: DateFrame, index is a list of timestamps at 2Hz, columns include HalfSecond_feature1, HalfSecond_feature2
    '''

    # Create wavelet dataframes
    oneSecond =
    halfSecond =

    # Compute wavelets
    cA_n, cD_3, cD_2, cD_1 = pywt.wavedec(eda, 'Haar', level=3) #3 = 1Hz, 2 = 2Hz, 1=4Hz

    # Wavelet 1 second window
    N = int(len(eda)/sampling_rate)
    coeff1 = np.max(abs(np.reshape(cD_1[0:4*N],(N,4))), axis=1)
    coeff2 = np.max(abs(np.reshape(cD_2[0:2*N],(N,2))), axis=1)
    coeff3 = abs(cD_3[0:N])
    wave1Second = pd.DataFrame({'OneSecond_feature1':coeff1,'OneSecond_feature2':coeff2,'OneSecond_feature3':coeff3})
    wave1Second.index = oneSecond[:len(wave1Second)]

    # Wavelet Half second window
    N = int(np.floor((len(data)/8.0)*2))
    coeff1 = np.max(abs(np.reshape(cD_1[0:2*N],(N,2))),axis=1)
    coeff2 = abs(cD_2[0:N])
    waveHalfSecond = pd.DataFrame({'HalfSecond_feature1':coeff1,'HalfSecond_feature2':coeff2})
    waveHalfSecond.index = halfSecond[:len(waveHalfSecond)]

    return wave1Second,waveHalfSecond 

def dwt_features(data, trials, level, sampling_freq, w, n, wavelet):
    """Extract discrete wavelet features

    Args:
        data: 2D (time points, Channels)
        trials: Trials vector
        lLevel: level of DWT decomposition
        sampling_freq: Sampling frequency

    Returns:
        The features matrix (Nbre trials, Nbre features)
    """

    # number of features per trial
    n_features = 9
    # allocate memory to store the features
    X = np.zeros((len(trials), n_features))

    # Extract Features
    for t, trial in enumerate(trials):
        signals = data[trial + fs*4 + (w[0]) : trial + fs*4 + (w[1])]
        coeffs_c3 = pywt.wavedec(data = signals[:,0], wavelet=wavelet, level=level)
        coeffs_c4 = pywt.wavedec(data = signals[:,1], wavelet=wavelet, level=level)
        coeffs_cz = pywt.wavedec(data = signals[:,2], wavelet=wavelet, level=level)

        X[t, :] = np.array([
            np.std(coeffs_c3[n]), np.mean(coeffs_c3[n]**2),
            np.std(coeffs_c4[n]), np.mean(coeffs_c4[n]**2),
            np.std(coeffs_cz[n]), np.mean(coeffs_cz[n]**2),
            np.mean(coeffs_c3[n]),
            np.mean(coeffs_c4[n]),
            np.mean(coeffs_cz[n])])

    return X 

def dwt(raw_eeg_data, level, **kwargs):
    """Multilevel Discrete Wavelet Transform (DWT).

    Compute the DWT for a raw eeg signal on multiple levels.

    Args:
        raw_eeg_data (array_like): input data
        level (int >= 0): decomposition levels
        **kwargs: Additional arguments that will be forwarded to ``pywt.wavedec``

    Returns:
        A 2-element tuple containing

        - **float**: mean value of the first decomposition coefficients
        - **list**: list of mean values for the individual (detail) decomposition coefficients

    """
    wt_coeffs = pywt.wavedec(data = raw_eeg_data, level=level, **kwargs)

    # A7:  0 Hz - 1 Hz
    cAL_mean = np.nanmean(wt_coeffs[0], axis=0)
    details = []

    # For Fs = 128 H
    for i in range(1, level+1):
        # D7:  1 Hz - 2 Hz
        cDL_mean = np.nanmean(wt_coeffs[i], axis=0)
        details.append(cDL_mean)

    return cAL_mean, details 

def wavelet_smoother(x, wavelet="db4", level=1, title=None):
    coeff = pywavelets.wavedec(x, wavelet, mode="per")
    sigma = mad(coeff[-level])
    uthresh = sigma * np.sqrt(2 * np.log(len(x)))
    coeff[1:] = (pywavelets.threshold(i, value=uthresh, mode="soft") for i in coeff[1:])
    return pywavelets.waverec(coeff, wavelet, mode="per") 

def test_compare_downcoef_coeffs():
    rstate = np.random.RandomState(1234)
    r = rstate.randn(16)
    # compare downcoef against wavedec outputs
    for nlevels in [1, 2, 3]:
        for wavelet in pywt.wavelist():
            wavelet = pywt.DiscreteContinuousWavelet(wavelet)
            if isinstance(wavelet, pywt.Wavelet):
                max_level = pywt.dwt_max_level(r.size, wavelet.dec_len)
                if nlevels <= max_level:
                    a = pywt.downcoef('a', r, wavelet, level=nlevels)
                    d = pywt.downcoef('d', r, wavelet, level=nlevels)
                    coeffs = pywt.wavedec(r, wavelet, level=nlevels)
                    assert_allclose(a, coeffs[0])
                    assert_allclose(d, coeffs[1]) 

def test_waverec_shape_mismatch_error():
    c = pywt.wavedec(np.ones(16), 'haar')
    # truncate a detail coefficient to an incorrect shape
    c[3] = c[3][:-1]
    assert_raises(ValueError, pywt.waverec, c, 'haar', axis=1) 

def test_waverec_axis_error():
    c = pywt.wavedec(np.ones(4), 'haar')
    # out of range axis not allowed
    assert_raises(ValueError, pywt.waverec, c, 'haar', axis=1) 

def test_wavedec_axis_error():
    data = np.ones(4)
    # out of range axis not allowed
    assert_raises(ValueError, pywt.wavedec, data, 'haar', axis=1) 

def test_waverec_axes_subsets():
    rstate = np.random.RandomState(0)
    data = rstate.standard_normal((8, 8, 8))
    # test all combinations of 1 out of 3 axes transformed
    for axis in [0, 1, 2]:
        coefs = pywt.wavedec(data, 'haar', axis=axis)
        rec = pywt.waverec(coefs, 'haar', axis=axis)
        assert_allclose(rec, data, atol=1e-14) 

def compute_wavelet_descriptor(beat, family, level):
    wave_family = pywt.Wavelet(family)
    coeffs = pywt.wavedec(beat, wave_family, level=level)
    return coeffs[0]

# Compute my descriptor based on amplitudes of several intervals 

def test_waverec_all_wavelets_modes():
    # test 2D case using all wavelets and modes
    rstate = np.random.RandomState(1234)
    r = rstate.randn(80)
    for wavelet in wavelist:
        for mode in pywt.Modes.modes:
            coeffs = pywt.wavedec(r, wavelet, mode=mode)
            assert_allclose(pywt.waverec(coeffs, wavelet, mode=mode),
                            r, rtol=tol_single, atol=tol_single)

####
# 2d multilevel dwt function tests
#### 

def test_waverec_complex():
    x = np.array([3, 7, 1, 1, -2, 5, 4, 6])
    x = x + 1j
    coeffs = pywt.wavedec(x, 'db1')
    assert_allclose(pywt.waverec(coeffs, 'db1'), x, rtol=1e-12) 

def test_waverec_odd_length():
    x = [3, 7, 1, 1, -2, 5]
    coeffs = pywt.wavedec(x, 'db1')
    assert_allclose(pywt.waverec(coeffs, 'db1'), x, rtol=1e-12) 

def test_waverec_none():
    x = [3, 7, 1, 1, -2, 5, 4, 6]
    coeffs = pywt.wavedec(x, 'db1')

    # set some coefficients to None
    coeffs[2] = None
    coeffs[0] = None
    assert_(pywt.waverec(coeffs, 'db1').size, len(x)) 

def test_waverec_accuracies():
    rstate = np.random.RandomState(1234)
    x0 = rstate.randn(8)
    for dt, tol in dtypes_and_tolerances:
        x = x0.astype(dt)
        if np.iscomplexobj(x):
            x += 1j*rstate.randn(8).astype(x.real.dtype)
        coeffs = pywt.wavedec(x, 'db1')
        assert_allclose(pywt.waverec(coeffs, 'db1'), x, atol=tol, rtol=tol) 

def test_wavedec():
    x = [3, 7, 1, 1, -2, 5, 4, 6]
    db1 = pywt.Wavelet('db1')
    cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
    assert_almost_equal(cA3, [8.83883476])
    assert_almost_equal(cD3, [-0.35355339])
    assert_allclose(cD2, [4., -3.5])
    assert_allclose(cD1, [-2.82842712, 0, -4.94974747, -1.41421356])
    assert_(pywt.dwt_max_level(len(x), db1) == 3) 

def denoise(X,wave0):
	wavelet=wave0
	if len(X)>=8:
		level0= 1
		if np.floor(np.log(len(X)))>7:
		   level0= np.floor(np.log(len(X))/2.0) 
		thres = 2*np.sqrt(2*np.log(len(X))/len(X))*np.std(X)
		thres = 0.0
		WaveletCoeffs = pywt.wavedec(X, wavelet, level=level0)
		NewWaveletCoeffs = map (lambda x: pywt.threshold(x, thres, mode='hard'),WaveletCoeffs)
		newWave2 = pywt.waverec( NewWaveletCoeffs, wavelet)
		return newWave2
	else:
		logging.warning("the series is too short!")
		return X 

def denoise(X,wave0):
    wavelet=wave0
    if len(X)>=8:
        level0= 1
        if np.floor(np.log(len(X)))>7:
           level0= np.floor(np.log(len(X))/2.0) 
        thres = 2*np.sqrt(2*np.log(len(X))/len(X))*np.std(X)
        thres = 0.0
        WaveletCoeffs = pywt.wavedec(X, wavelet, level=level0)
        NewWaveletCoeffs = map (lambda x: pywt.threshold(x, thres, mode='hard'),WaveletCoeffs)
        newWave2 = pywt.waverec( NewWaveletCoeffs, wavelet)
        return newWave2
    else:
        logging.warning( "the series is too short")
        return X


#compute the liquidity index 

def getWaveletData(values, waveletName, level, threadMethodName):
    mode = 'sym'
    #小波系数分解
    data = pywt.wavedec(values, waveletName, mode, level)
    #cA4, cD4, cD3, cD2, cD1 = pywt.wavedec(values, waveletName, mode, level)
    coeffs = [] #小波重构系数
    #阈值处理
    if threadMethodName == 'sqtwolog':
        #print len(cA4), len(cD4), len(cD3), len(cD2), len(cD1),len(values)
        for i in np.arange(0, len(data)):
            if i > 0:
                data[i] = softThreshold(threshold_sqtwolog(len(data[i])), data[i])
            coeffs.append(data[i])
    #小波重构
    zValues = pywt.waverec(coeffs, waveletName, mode)
        
#     if level == 4:
#         #cA4, cD4, cD3, cD2, cD1 = pywt.wavedec(values, waveletName, mode, level)
#         coeffs = []
#         #阈值处理
#         if threadMethodName == 'sqtwolog':
#             #print len(cA4), len(cD4), len(cD3), len(cD2), len(cD1),len(values)
# #             cD4 = softThreshold(threshold_sqtwolog(len(cD4)), cD4)
# #             cD3 = softThreshold(threshold_sqtwolog(len(cD3)), cD3)
# #             cD2 = softThreshold(threshold_sqtwolog(len(cD2)), cD2)
# #             cD1 = softThreshold(threshold_sqtwolog(len(cD1)), cD1)
#         #小波重构
#         #coeffs = [cA4, cD4, cD3, cD2, cD1]
#         zValues = pywt.waverec(coeffs, waveletName, mode)
#     elif level ==2:
#         cA2, cD2, cD1 = pywt.wavedec(values, waveletName, mode, level)
#         #阈值处理
#         if threadMethodName == 'sqtwolog':
#             print len(cA2), len(cD2), len(cD1),len(values)
#             cD2 = softThreshold(threshold_sqtwolog(len(cD2)), cD2)
#             cD1 = softThreshold(threshold_sqtwolog(len(cD1)), cD1)
#         #小波重构
#         coeffs = [cA2, cD2, cD1]
#         zValues = pywt.waverec(coeffs, waveletName, mode) 
    return zValues

#小波包分解
# forecastCount:预测的点位数 

def get_data_label_mitdb( list_patient, mit_db ):
    labels = np.array([], dtype=np.int32)
    data = np.array([], dtype=float)

    for p in list_patient:
        index = mit_db.patients.index(str(p))
        for b in range(0, len(mit_db.classes[index]), 1):
            RR = [mit_db.temporal_features[index].pre_R[b], mit_db.temporal_features[index].post_R[b], mit_db.temporal_features[index].local_R[b], mit_db.temporal_features[index].global_R[b]]
            signal = mit_db.signals[index][b]
            beat_type = mit_db.classes[index][b]
            # Name class by numbers np.int32
            #['N', 'L', 'R', 'e', 'j', 'A', 'a', 'J', 'S', 'V', 'E', 'F', 'P', '/', 'f', 'u']
            for i in range(0,5,1):
                if beat_type in superclass[i]:
                    class_n = i
                    break #exit loop
            
            labels = np.append(labels, class_n)
            
            #Display raw and wave signal
            if not compute_wavelets:
                features = signal

            else: # db of order 8 
                db8 = pywt.Wavelet('db8')
                coeffs = pywt.wavedec(signal, db8, level=4)
                features = coeffs[1]

            #TODO add RR interval
            if compute_RR_interval_feature:
                features = np.append(features, RR)

            if len(data) == 0:
                data = features
            else:
                data = np.vstack((data, features))           
        
                #plt.subplot(211)
                #plt.plot(signal)
                #plt.subplot(212)
                #plt.plot(coeffs[1])
                #plt.show()
    return (data, labels)