Python pywt.waverec() Examples

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_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 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 test_first_stage(self):
        """ Test of perfect reconstruction of first stage wavelets. """
        for wavelet in available_first_stage_filters():
            # Using waverec and wavedec instead of dwt and idwt because parameters
            # don't need as much parsing.
            self.assertTrue(
                np.allclose(
                    self.array, pywt.waverec(pywt.wavedec(self.array, wavelet), wavelet)
                )
            )


##############################################################################
###           DUAL-TREE COMPLEX WAVELET TRANSFORM
############################################################################## 

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_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_db2():
    # test for fix to issue gh-293
    rstate = np.random.RandomState(0)
    data = rstate.standard_normal((16, 16))
    for axis in [0, 1]:
        coefs = pywt.wavedec(data, 'db2', axis=axis)
        rec = pywt.waverec(coefs, 'db2', axis=axis)
        assert_allclose(rec, data, atol=1e-14) 

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 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 test_multilevel_dtypes_1d():
    # only checks that the result is of the expected type
    wavelet = pywt.Wavelet('haar')
    for dt_in, dt_out in zip(dtypes_in, dtypes_out):
        # wavedec, waverec
        x = np.ones(8, dtype=dt_in)
        errmsg = "wrong dtype returned for {0} input".format(dt_in)

        coeffs = pywt.wavedec(x, wavelet, level=2)
        for c in coeffs:
            assert_(c.dtype == dt_out, "wavedec: " + errmsg)
        x_roundtrip = pywt.waverec(coeffs, wavelet)
        assert_(x_roundtrip.dtype == dt_out, "waverec: " + errmsg) 

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_waverec_invalid_inputs():
    # input must be list or tuple
    assert_raises(ValueError, pywt.waverec, np.ones(8), 'haar')

    # input list cannot be empty
    assert_raises(ValueError, pywt.waverec, [], 'haar') 

def plot_signal_decomp(data, w, title):
    """Decompose and plot a signal S.
    S = An + Dn + Dn-1 + ... + D1
    """
    w = pywt.Wavelet(w)
    a = data
    ca = []
    cd = []
    for i in range(5):
        (a, d) = pywt.dwt(a, w, mode)
        ca.append(a)
        cd.append(d)

    rec_a = []
    rec_d = []

    for i, coeff in enumerate(ca):
        coeff_list = [coeff, None] + [None] * i
        rec_a.append(pywt.waverec(coeff_list, w))

    for i, coeff in enumerate(cd):
        coeff_list = [None, coeff] + [None] * i
        rec_d.append(pywt.waverec(coeff_list, w))

    fig = plt.figure()
    ax_main = fig.add_subplot(len(rec_a) + 1, 1, 1)
    ax_main.set_title(title)
    ax_main.plot(data)
    ax_main.set_xlim(0, len(data) - 1)

    for i, y in enumerate(rec_a):
        ax = fig.add_subplot(len(rec_a) + 1, 2, 3 + i * 2)
        ax.plot(y, 'r')
        ax.set_xlim(0, len(y) - 1)
        ax.set_ylabel("A%d" % (i + 1))

    for i, y in enumerate(rec_d):
        ax = fig.add_subplot(len(rec_d) + 1, 2, 4 + i * 2)
        ax.plot(y, 'g')
        ax.set_xlim(0, len(y) - 1)
        ax.set_ylabel("D%d" % (i + 1)) 

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 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:预测的点位数