Re: How to rebin complex array? [message #39011] |
Wed, 14 April 2004 13:01 |
Yunxiang Zhang
Messages: 19 Registered: October 2003
|
Junior Member |
|
|
Folks,
I did some test to rebin a complex array. Thanks for all you guys' help!:)
speed test result: a for loop is not as worse as I would imagine in
this case.
"#" > "for loop" > "special rebin" > "component rebin"
for loop 0.25066495
component rebin 0.53236008
special rebin 0.37491202
# 0.16867399
I also found using fltarr is better than using replicate
fltarr > replicate
replicate 0.039183140
fltarr 0.027969122
;here's my testing code
aa=bytscl(reform(ran2_normal(-738345,dim=5000),2,2500))
a=reform(complex(aa[0,*],aa[1,*]))
na=n_elements(a) & k=na
;using for loop
t0=systime(1)
b1=complexarr(na,k)
for i=0,k-1 do begin
b1[*,i]=a
endfor
print, 'for loop', systime(1)-t0
;using component rebin
t0=systime(1)
ra=real_part(a)
ia=imaginary(a)
rb=rebin(ra,[na,k])
ib=rebin(ia,[na,k])
b2=complex(rb,ib)
print,'component rebin', systime(1)-t0
;i call it special rebin, please refer to Peter's post
t0=systime(1)
d=double(a,0,k)
r=rebin(temporary(d),na,k,/sample)
b3=complex(temporary(r),0,na,k)
print,'special rebin', systime(1)-t0
;using #
t0=systime(1)
b4=a#(1+fltarr(k))
print,'#', systime(1)-t0
t0=systime(1)
for i=0,1000 do aaa=replicate(1.,10000)
print,'replicate', systime(1)-t0
t0=systime(1)
for i=0,1000 do bbb=1+fltarr(10000)
print,'fltarr', systime(1)-t0
end
|
|
|
Re: How to rebin complex array? [message #39013 is a reply to message #39011] |
Wed, 14 April 2004 09:31  |
R.G. Stockwell
Messages: 363 Registered: July 1999
|
Senior Member |
|
|
"James Kuyper" <kuyper@saicmodis.com> wrote in message news:407D40F3.F045EC70@saicmodis.com...
> Timm Weitkamp wrote:
> ...
...
>> By the way, there is a good reason for REBIN not to work on complex
>> arrays, and that is the fact that there is not really any proper way of
>> interpolating between two complex values.
>
> Why in the world not?
In the 2D complex space, it is ambiguous as how to interpolate.
Breifly, one could interpolate real and imaginary seperately, as
suggested in this thread. Or one could average the lengths of
the vector, and the angle of the vector.
For any rotating vector time series, component-wise interpolation
will almost always underestimate the amplitude of the interpolated value.
Consider two subsequent values, [1,0] and [0,1].
If this is a unit vector rotating, then the "correct" interpolated
value would be [1/root 2,1/root 2] = [.707,.707]. The componently
interpolated value is [.5,.5].
So it really depends on what the type of time series one has.
It may be that a "norm preserving interpolation" may be required.
(average the lengths of the vector rather than just average components).
Cheers,
bob
|
|
|
Re: How to rebin complex array? [message #39016 is a reply to message #39013] |
Wed, 14 April 2004 06:47  |
James Kuyper
Messages: 425 Registered: March 2000
|
Senior Member |
|
|
Timm Weitkamp wrote:
...
> u = 1 + FLTARR(k) ; create unit vector with k elements
A unit vector is a vector such that total(u*u) eq 1.0" (or for complex
numbers, total(u*conj(u)) eq 1.0). What you've described is a vector
whose elements are all 1, which has no special name.
...
> By the way, there is a good reason for REBIN not to work on complex
> arrays, and that is the fact that there is not really any proper way of
> interpolating between two complex values.
Why in the world not? All of the normal interpolation formulas work
perfectly well when applied seperately to the real and imaginary
components of a table of complex values. The formulas just involve
multiplcation, division, addition, and subtraction, all of which are
perfectly well defined for complex numbers, and they produce reasonable
results (well, not quite - in some higher-order interpolation algorithms
it's necessary to replace c^2 with c*conj(c) in some places).
|
|
|
Re: How to rebin complex array? [message #39019 is a reply to message #39016] |
Wed, 14 April 2004 00:39  |
Timm Weitkamp
Messages: 66 Registered: August 2002
|
Member |
|
|
On 13.04.04 at 15:57 -0700, Yunxiang Zhang wrote:
> I wanna do the following thing as fast as possible.
>
> ;creat a complex array
> a=complexarr(na)
> ;replicate the array k times
> b=rebin(a,na,k)
>
> But rebin won't accept complex array. :( What should I do?
This here is probably what you want:
a=complexarr(na) ; create a complex array
u = 1 + FLTARR(k) ; create unit vector with k elements
b = a # u ; matrix multiplication of a and u
By the way, there is a good reason for REBIN not to work on complex
arrays, and that is the fact that there is not really any proper way of
interpolating between two complex values. (You could argue that
next-neighbor sampling should still be possible.)
> Another similiar question is what is the best solution to replicate a 2d
> image eg dist(512,512) k times to get a k-frame "still movie" as a 3d array
> movie(512,512,k) without using any loops?
As previous posters have said, there's no reason why REBIN shouldn't work
here, given that DIST yields a real-valued result. But *if* you did have a
complex array that you wanted to replicate in the manner described above,
or if you want to avoid REBIN for some other reason, and not use loops,
then I'd probably do this:
s = SIZE(image)
movie = MAKE_ARRAY(s[1], s[2], TYPE=s[3])
movie[*] = image[*] # (1.0 + FLTARR(k))
Cheers,
Timm
--
Timm Weitkamp <http://people.web.psi.ch/weitkamp>
|
|
|
Re: How to rebin complex array? [message #39021 is a reply to message #39019] |
Tue, 13 April 2004 17:47  |
Peter Mason
Messages: 145 Registered: June 1996
|
Senior Member |
|
|
Yunxiang Zhang wrote:
> Folks,
>
> I wanna do the following thing as fast as possible.
>
> ;creat a complex array
> a=complexarr(na)
> ;replicate the array k times
> b=rebin(a,na,k)
>
> But rebin won't accept complex array. :( What should I do?
Further to what Mark has written, if this *really* is what you want to do
and if you can live with REBIN's /SAMPLE option (which is faster than its
default bilinear interpolation anyway), then you might try something like
the following:
d=double( a, 0, na )
r=rebin( temporary(d), na, k, /SAMPLE )
b=complex( temporary(r), 0, na, k )
Of course, this'll only work for single-precision complex "a". (8 bytes
per value, same as "non-complex" double precision.)
Peter Mason
|
|
|
Re: How to rebin complex array? [message #39022 is a reply to message #39021] |
Tue, 13 April 2004 17:10  |
Mark Hadfield
Messages: 783 Registered: May 1995
|
Senior Member |
|
|
Yunxiang Zhang wrote:
>
> I wanna do the following thing as fast as possible.
>
> ;creat a complex array
> a=complexarr(na)
> ;replicate the array k times
> b=rebin(a,na,k)
>
> But rebin won't accept complex array. :( What should I do?
Extract the real & imaginary parts, rebin each one & recombine them?
> Another similiar question is what is the best solution to replicate a 2d
> image eg dist(512,512) k times to get a k-frame "still movie" as a 3d array
> movie(512,512,k) without using any loops?
You can use REBIN in this case.
You might want to look at JD's dimension juggling tutorial, on David
Fanning's site:
http://www.dfanning.com/tips/rebin_magic.html
--
Mark Hadfield "Ka puwaha te tai nei, Hoea tatou"
m.hadfield@niwa.co.nz
National Institute for Water and Atmospheric Research (NIWA)
|
|
|