Python pywt.dwt_max_level函数代码示例

def test_default_level():
    # default level is the maximum permissible for the transformed axes
    data = np.ones((128, 32, 4))
    wavelet = ('db8', 'db1')
    for dec_func in [pywt.wavedec2, pywt.wavedecn]:
        for axes in [(0, 1), (2, 1), (0, 2)]:
            c = dec_func(data, wavelet, axes=axes)
            max_lev = np.min([pywt.dwt_max_level(data.shape[ax], wav)
                              for ax, wav in zip(axes, wavelet)])
            assert_equal(len(c[1:]), max_lev)

    for ax in [0, 1]:
        c = pywt.wavedecn(data, wavelet[ax], axes=(ax, ))
        assert_equal(len(c[1:]),
                     pywt.dwt_max_level(data.shape[ax], wavelet[ax]))

def my_waverec(data, wavelet, mode='sym', level=None):
    """
    Multilevel 1D Discrete Wavelet Transform of data.
    Returns coefficients list - [cAn, cDn, cDn-1, ..., cD2, cD1]

    data    - input data
    wavelet - wavelet to use (Wavelet object or name string)
    mode    - signal extension mode, see MODES
    level   - decomposition level. If level is None then it will be
              calculated using `dwt_max_level` function.
    """

    if not isinstance(wavelet, pywt.Wavelet):
        wavelet = pywt.Wavelet(wavelet)

    if level is None:
        level = pywt.dwt_max_level(len(data), wavelet.dec_len)
    elif level < 0:
        raise ValueError("Level value of %d is too low . Minimum level is 0." % level)

    a = data[0:2]
    d = data[2:4]
    for i in xrange(int(level)):
        offs = pow(2,i+2)
        a = pywt.idwt(a, d, wavelet, mode)
        d = data[offs:2*offs]

    return a

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 my_wavedec(data, wavelet, mode='sym', level=None):
    """
    Multilevel 1D Discrete Wavelet Transform of data.
    Returns coefficients list - [cAn, cDn, cDn-1, ..., cD2, cD1]

    data    - input data
    wavelet - wavelet to use (Wavelet object or name string)
    mode    - signal extension mode, see MODES
    level   - decomposition level. If level is None then it will be
              calculated using `dwt_max_level` function.
    """

    if not isinstance(wavelet, pywt.Wavelet):
        wavelet = pywt.Wavelet(wavelet)

    if level is None:
        level = pywt.dwt_max_level(len(data), wavelet.dec_len)
    elif level < 0:
        raise ValueError("Level value of %d is too low . Minimum level is 0." % level)

    coeffs_list = []

    a = data
    for i in xrange(level):
        a, d = pywt.dwt(a, wavelet, mode)
        d = list(d)
        #d.reverse()
        coeffs_list.append(d)

    a = list(a)
    #a.reverse()
    coeffs_list.append(a)
    #coeffs_list.reverse()

    return coeffs_list

	def filterData(self, icurr, Fs):
		"""
			Denoise an ionic current time-series and store it in self.eventData

			:Parameters:
				- `icurr` :	ionic current in pA
				- `Fs` :	original sampling frequency in Hz
		"""
		# self.eventData=icurr
		self.Fs=Fs

		# Set up the wavelet
		w=pywt.Wavelet(self.waveletType)

		# Calculate the maximum wavelet level for the data length
		self.maxWaveletLevel=pywt.dwt_max_level(len(icurr), filter_len=w.dec_len)

		# Perform a wavelet decomposition to the specified level
		wcoeff = pywt.wavedec(icurr, w, mode='sym', level=self.waveletLevel)

		# Perform a simple threshold by setting all the detailed coefficients
		# up to level n-1 to zero
		thresh=np.std(wcoeff[-1])*self._thselect(wcoeff, self.waveletThresholdSubType)
		thrfunc=self.thrtypedict[self.waveletThresholdType]

		# print thresh, np.std(wcoeff[-1])
		wcoeff[1:] = [ thrfunc(wc, thresh) for wc in wcoeff[1:] ]
		# for i in range(1, self.waveletLevel):
		# 	wcoeff[-i]=np.zeros(len(wcoeff[-i]))

		# Reconstruct the signal with the thresholded wavelet coefficients
		self.eventData = pywt.waverec(wcoeff, self.waveletType, mode='sym')

def c_dists(Y,use_swt=True,level_weights=False):
	w = pywt.Wavelet('sym2')
	if use_swt:
		L = pywt.swt_max_level(Y.shape[0])
		C = [pywt.swt(Y[:,i],w,level=L) for i in range(Y.shape[1])]
		C = [[list(reshape(l[0],-1)) + list(reshape(l[1],-1)) for l in c] for c in C]
	else:
		L = pywt.dwt_max_level(Y.shape[0],w)
		C = [pywt.wavedec(Y[:,i],w,level=L) for i in range(Y.shape[1])]
	if level_weights:
		if use_swt:
			raise NameError('No level weights with SWT')
		Wc = [1. for x in range(1,L+1)]
		D = zeros((len(C),len(C)))
		for i in range(len(C)):
			for j in range(i+1,len(C)):
				d = sum([distance.cosine(C[i][x],C[j][x])*Wc[x] for x in range(L)])/sum(Wc)
				D[i,j] = d
				D[j,i] = d
		return D
	else:
		Cn = []
		for c in C:
			cn = []
			for l in c:
				cn += list(l)
			Cn.append(cn)
		return abs(pdist(Cn,'cosine'))

def wp(S, costf, wavelet="db4", mode=pywt.MODES.ppd, level=None):
    '''
    Returns the 1D discrete wavelet packet transformation, with the best basis according
    to the given cost function, for the given 1D input signal.
    @param S:         Input signal.
                      Both single and double precision floating-point data types are supported
                      and the output type depends on the input type. If the input data is not
                      in one of these types it will be converted to the default double precision
                      data format before performing computations.
    @param costf:      The (single parameter) cost function that must be used while
                      searching for the best basis.
    @param wavelet:   Wavelet to use in the transform. 
                      This must be a name of the wavelet from the wavelist() list.
    @param mode:      Signal extension mode to deal with the border distortion problem.
                      The default mode is periodic-padding.
    @param level:     Number of decomposition steps to perform. If the level is None, then the
                      full decomposition up to the level computed with dwt_max_level() function for
                      the given data and wavelet lengths is performed.
    @return:          A list containing the nodes of the 1D discrete wavelet packet transformation,
                      with the best basis according to the given cost function, for the given input signal. 
    '''
    if (level == None):
        level = pywt.dwt_max_level(S.shape[0], pywt.Wavelet(wavelet))
    
    #Data collection step
    Nodes = collect(S, wavelet=wavelet, mode=mode, level=level)
    #Dynamic programming upstream traversal
    mark(Nodes, costf)
    #node.print_nodes(Nodes)
    #Dynamic programming downstream traversal
    Result = []
    traverse(Nodes[0][0], Nodes, Result)
    traverse(Nodes[0][1], Nodes, Result)
    return sorted(Result, cmp=node.compare_low_level_first, reverse=False)

def test_wavelet_denoising_levels():
    rstate = np.random.RandomState(1234)
    ndim = 2
    N = 256
    wavelet = 'db1'
    # Generate a very simple test image
    img = 0.2*np.ones((N, )*ndim)
    img[[slice(5, 13), ] * ndim] = 0.8

    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
                                             wavelet_levels=1)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    psnr_denoised_1 = compare_psnr(img, denoised_1)

    # multi-level case should outperform single level case
    assert_(psnr_denoised > psnr_denoised_1 > psnr_noisy)

    # invalid number of wavelet levels results in a ValueError
    max_level = pywt.dwt_max_level(np.min(img.shape),
                                   pywt.Wavelet(wavelet).dec_len)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=max_level+1)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=-1)

def wavelet_trans(l, keep_proportion = 0.5):
    #TODO binning per decomposition level?
    #TODO explain why db4
    if not (l is None or l == []):
        coeffs = pywt.wavedec(l, 'db4', level=pywt.dwt_max_level(len(l),pywt.Wavelet('db4')))
        return merge(*coeffs[0:int(len(coeffs)*keep_proportion)])
    else:
        return l

    def full_coeff_len(datalen, filtlen, mode):
        max_level = wt.dwt_max_level(datalen, filtlen)
        total_len = 0

        for i in xrange(max_level):
            datalen = wt.dwt_coeff_len(datalen, filtlen, mode)
            total_len += datalen

        return total_len + datalen

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 startup():
    global w, N, L0, mode, level

    print pywt.families()
    print pywt.wavelist('db')

    w = pywt.Wavelet('db6')
    mode = pywt.MODES.per

    print w
    print "vanishing_moments_psi:", w.vanishing_moments_psi
    print "vanishing_moments_phi:", w.vanishing_moments_phi

    N = 2**9
    print "max level = ", pywt.dwt_max_level(N, w.dec_len)
    L0 = numpy.zeros((N,N), 'double')
    if True:
        for i in xrange(0,N):
            L0[i][i-1], L0[i][i], L0[i-1][i] = (1., -2., 1.)
    else:
        for i in xrange(1,N):
            L0[i][i-1], L0[i][i], L0[i-1][i] = (1., -2., 1.)
        L0[0][0] = -2.

#L0[0][N-1], L0[0][0], L0[N-1][0] = (1, -2, 1)

#L0 = numpy.eye(N)

#for i in xrange(0,N):
        #L0[i] = [1,2,3,4,5,6,7,8]

    numpy.core.arrayprint.set_printoptions(threshold=N*N+1, linewidth=100000)

#print L0

#coeffs = pywt.wavedec2(L0, w) #, level=pywt.dwt_max_level(N, w.dec_len))
#print coeffs
#print pywt.waverec2(coeffs, w)
#
#coeffs = pywt.wavedec2(L0, w) #, level=pywt.dwt_max_level(N, w.dec_len))
#print coeffs

    print "max level = ", pywt.dwt_max_level(N, w.dec_len)

def wavelet_levels(Y):
	w = pywt.Wavelet('sym2')
	levels = pywt.dwt_max_level(Y.shape[0],w)
	w0 = pywt.wavedec(Y[:,0],w,level=levels)[1:]
	L = [np.empty((Y.shape[1],len(x))) for x in w0]
	for i in range(Y.shape[1]):
		wd = pywt.wavedec(Y[:,i],w)[1:]
		for j,x in enumerate(wd):
			L[j][i,:] = x
	return L,[Y.shape[0]/len(x) for x in w0]

def haar_decomp(img_file_name):
    img = Image.open(img_file_name).convert('L')
    x,y = img.size
    coeffs = pywt.wavedec2(img, 'haar', level=pywt.dwt_max_level(max(x,y), pywt.Wavelet('haar')))
    params = coeffs[0]
    final_arr = []
    for i in xrange(1,len(coeffs)):
	    for j in coeffs[i]:
		    final_arr = np.concatenate((final_arr, np.ndarray.flatten(j)))
    return params, coeffs[1:], final_arr

def wavelet_avg(Y,X=None,reshape=False,plot_xy=False,zero_thresh=1*10**-7,Names=None,Title=None):
	if not X:
		X = range(Y.shape[0])
	w = pywt.Wavelet('sym2')
	L = pywt.dwt_max_level(Y.shape[0],w)
	#L = pywt.swt_max_level(Y.shape[0])
	Zavg = zeros((L,len(X)))
	Zavgabs = zeros((L,len(X)))
	if reshape:
		e = [0,L,0,L]
	else:
		e = [X[0],X[-1],0,L]
	sp = 1
	if not Names:
		Names = ["Series "+str(y+1) for y in range(Y.shape[1])]
	if Title:
		pylab.suptitle(Title)
	if plot_xy:
		tp = str(Y.shape[1]+3)
		pylab.subplot(int(''.join([tp,'1',str(sp)])))
		P = pylab.plot(X,Y)
		pylab.legend(P,Names,prop={'size':6})
		pylab.xlim([0,Y.shape[0]-1])
		pylab.title("Relative Abundance Time Series")
		sp += 1
	else:
		tp = str(Y.shape[1]+2)
	for y in range(Y.shape[1]):
		C = pywt.wavedec(Y[:,y],w,level=L)
		#C = pywt.swt(Y[:,y],w,level=L)
		#C = ['dummy'] + [c[1] for c in C]
		Z = wavelet_matrix(C,L,len(X))
		pylab.subplot(int(''.join([tp,'1',str(sp)])))
		pylab.imshow(abs(Z),extent=e)
		pylab.title(Names[y])
		sp += 1
		Zavg += Z
		Zavgabs += abs(Z)
	Zavg /= Y.shape[1]
	Zavg[Zavg<zero_thresh] = 0
	pylab.subplot(int(''.join([tp,'1',str(sp)])))
	pylab.imshow(abs(Zavg),extent=e)
	pylab.title("Decomposition Avg")
	sp += 1
	Zavgabs /= Y.shape[1]
	Zavgabs[Zavgabs<zero_thresh] = 0
	pylab.subplot(int(''.join([tp,'1',str(sp)])))
	pylab.imshow(abs(Zavgabs),extent=e)
	pylab.title("abs(Decomposition) Sum")
	sp += 1
	pylab.show()
	return Zavg