Daubechies D4小波变换

原文算法

import sys
from math import sqrt
h0=(1+sqrt(3))/(4*sqrt(2))
h1=(3+sqrt(3))/(4*sqrt(2))
h2=(3-sqrt(3))/(4*sqrt(2))
h3=(1-sqrt(3))/(4*sqrt(2))

g0=h3
g1=-h2
g2=h1
g3=-h0


Ih0 = h2
Ih1 = g2
Ih2 = h0
Ih3 = g0

Ig0 = h3
Ig1 = g3
Ig2 = h1
Ig3 = g1

def transform(a,n):

	if n>=4:
		half=n>>1
		tmp=[]
		j=0
		i=0
		while j<n-3:
			tmp.insert(i,a[j]*h0+a[j+1]*h1+a[j+2]*h2+a[j+3]*h3)
			tmp.insert(i+half,a[j]*g0+a[j+1]*g1+a[j+2]*g2+a[j+3]*g3)
			i+=1
			j+=2
		tmp.insert(i,a[n-2]*h0+a[n-1]*h1+a[0]*h2+a[1]*h3)
		tmp.insert(i+half,a[n-2]*g0+a[n-1]*g1+a[0]*g2+a[1]*g3)

		for i in range(n):
			a[i]=tmp[i]
			print('a[%d]=:%4.2f' % (i,a[i]))
		return a

def iTransform(a,n):

	if n>=4:
		i=0
		j=2
		half=n>>1
		halfPls1=half+1
		tmp=[]

		tmp.append(a[half-1]*Ih0+a[n-1]*Ih1+a[0]*Ih2+a[half]*Ih3)
		tmp.append(a[half-1]*Ig0+a[n-1]*Ig1+a[0]*Ig2+a[half]*Ig3)

		for i in range(half-1):
			tmp.insert(j,a[i]*Ih0+a[i+half]*Ih1+a[i+1]*Ih2+a[i+halfPls1]*Ih3)
			j+=1
			tmp.insert(j,a[i]*Ig0+a[i+half]*Ig1+a[i+1]*Ig2+a[i+halfPls1]*Ig3)
			j+=1
		for i in range(n):
			a[i]=tmp[i]
			print('a[%d]=%4.2f' %(i,a[i]))
		return a


def main():
	a=[1265,2498,1634,2353,1828,2214,
	1904,2168,1949,2130,1976,2116,2003,2091,
	2002,2081,2008,2070]
	b=transform(a,18)
	t=input('input thread value:')

	for i in range(len(b)/2,len(b)):
		if b[i] < t:
			b[i]=0
	print('--------------------------')
	for bb in b:
		print(bb)

	iTransform(b,18)

if __name__ == '__main__':
	main()