On 16 Nov 2006 13:09:03 -0800, "sturlamolden" <[EMAIL PROTECTED]> wrote:
...SNIP... >To compare Matlab with NumPy we can e.g. use the D4 discrete wavelet >transform. I have here coded it in Matlab and Python/NumPy using Tim >Swelden's lifting scheme. > >First the Matlab version (D4_Transform.m): > >function x = D4_Transform(x) > % D4 Wavelet transform in Matlab > % (C) Sturla Molden > C1 = 1.7320508075688772; > C2 = 0.4330127018922193; > C3 = -0.066987298107780702; > C4 = 0.51763809020504137; > C5 = 1.9318516525781364; > s1 = zeros(ceil(size(x)/2)); > d1 = zeros(ceil(size(x)/2)); > d2 = zeros(ceil(size(x)/2)); > odd = x(2:2:end); > even = x(1:2:end-1); > d1(:) = odd - C2*even; > s1(1) = even(1) + C2*d1(1) + C3*d1(end); > s1(2:end) = even(2:end) + C2*d1(2:end) + C3*d1(1:end-1); > d2(1) = d1(1) + s1(end); > d2(2:end) = d1(2:end) + s1(1:end-1); > x(1:2:end-1) = C4*s1; > x(2:2:end) = C5*d2; > if (length(x) >2) > x(1:2:end-1) = D4_Transform(x(1:2:end-1)); > end > > >Then the Python version (D4.py): > >import numpy >import time > >def D4_Transform(x, s1=None, d1=None, d2=None): > """ > D4 Wavelet transform in NumPy > (C) Sturla Molden > """ > C1 = 1.7320508075688772 > C2 = 0.4330127018922193 > C3 = -0.066987298107780702 > C4 = 0.51763809020504137 > C5 = 1.9318516525781364 > if d1 == None: > d1 = numpy.zeros(x.size/2) > s1 = numpy.zeros(x.size/2) > d2 = numpy.zeros(x.size/2) > odd = x[1::2] > even = x[:-1:2] > d1[:] = odd[:] - C1*even[:] > s1[0] = even[0] + C2*d1[0] + C3*d1[-1] > s1[1:] = even[1:] + C2*d1[1:] + C3*d1[:-1] > d2[0] = d1[0] + s1[-1] > d2[1:] = d1[1:] + s1[:-1] > even[:] = C4 * s1[:] > odd[:] = C5 * d2[:] > if x.size > 2: > >D4_Transform(even,s1[0:even.size/2],d1[0:even.size/2],d2[0:even.size/2]) > >if __name__ == "__main__": > x = numpy.random.rand(2**23) > t0 = time.clock() > D4_Transform(x) > t = time.clock() > print "Elapsed time is %.6f seconds" % (t-t0) > ...SNIP... I couldn't help but notice that the following two lines from the above code are doing different things, but I think they are supposed to do the same things (the first is the Matlab version, the 2nd is the Python version): > d1(:) = odd - C2*even; > d1[:] = odd[:] - C1*even[:] I don't understand what this code is doing, so I don't know which is wrong. I don't even know if this is the actual code you are using. Just thought I'd mention it. Bill -- http://mail.python.org/mailman/listinfo/python-list
