Boss,you are absolutely right and I 'm sorry for the typos.It seems
the origin post also has some error according to the following
responses.I try to update the three versions as:
1.Numpy version:
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[:-1] = even[:-1] + C2*d1[:-1] + C3*d1[1:]
s1[-1] = even[-1] + C2*d1[-1] + C3*d1[0]
d2[0] = d1[0] + s1[-1]
d2[1:] = d1[1:] + s1[:-1]
even[:] = C5 * s1[:]
odd[:] = C4 * d2[:]
if x.size > 2:
D4_Transform(even,s1[0:even.size/2],d1[0:even.size/2],d2[0:even.size/2])
2.Matlab version:
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 - C1*even;
s1(1:end-1) = even(1:end-1) + C2*d1(1:end-1) + C3*d1(2:end);
s1(end) = even(end) + C2*d1(end) + C3*d1(1);
d2(1) = d1(1) + s1(end);
d2(2:end) = d1(2:end) + s1(1:end-1);
x(1:2:end-1) = C5*s1;
x(2:2:end) = C4*d2;
if (length(x) >2)
x(1:2:end-1) = D4_Transform(x(1:2:end-1));
end
3. J Version
D4_Transform=:3 :0
NB. D4 Wavelet transform in Matlab
NB. (C) Sturla Molden
C1 =. 1.7320508075688772
C2 =. 0.4330127018922193
C3 =. -0.066987298107780702
C4 =. 0.51763809020504137
C5 =. 1.9318516525781364
n=.>.-:#y
r =. (#y)$0
evenkey=.2*i.n
oddkey=.<<<evenkey
odd =.oddkey{y
even =.evenkey{y
d1 =. odd - C1*even
s1 =.even+(C2*d1)+C3*1|.d1
d2=.d1+_1|.s1
r=.(C5*s1)evenkey}r
r=.(C4*d2)oddkey}r
if.2<#y do.
r=.(D4_Transform evenkey{r) evenkey}r
end.
r
)
And using the code in
http://groups.google.com/group/comp.lang.python/msg/a71a5bf00a88ad57
to verify.
For input 1,5,3,2,4,8,5,2,the outputs are all 10.606602 2.6736336
_0.7589746 _0.29349422 1.272112 1.3541544 _4.8391016 _0.90586666.
The new timings are: numpy -- 1.71 seconds,matlab -- 25 seconds,J -- 20s.
2006/12/23, Mike Powell <[EMAIL PROTECTED]>:
The use of the constant C1 looks odd. It is used in the numpy
version, but not in the Matlab or J versions.
Mike Powell
On 22-Dec-06, at 6:01 AM, bill lam wrote:
> Xu Zuoqian wrote:
>> And the timing is:
>> d=.?(2^23)$0
>> ts 'D4_Transform d'
>> 150.59277 7.718007e8
>
> Before talking about efficiency could you give the correct result of
> D4_Transform d
> for validation, and also use
> d=. ?.(2^23)$0
> so that experiment is repeatable.
>
> --
> regards,
> bill
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
