About alternatives to Matlab

[Bill Maxwell]

On 16 Nov 2006 13:09:03 -0800, "sturlamolden" <sturlamolden at yahoo.no>
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