Thanks. I've re-written the function to minimize the amount of copying
(i.e. slicing) that is required. But now, I'm befuddled as to how to
parallelize this function using Julia. Any suggestions?
Alex
function expensive_hat(S::Array{Complex{Float64},2}, mx::Array{Int64,2},
my::Array{Int64,2})
samples = 64
A = zeros(size(mx));
B = zeros(size(mx));
C = zeros(size(mx));
for i = 1:samples
Si = S[:,i];
Sx = Si[mx];
Sy = Si[my];
Sxy = Si[mx+my];
Sxyc = conj(Sxy);
A += abs2(Sy .* Sx);
B += abs2(sqrt(Sxyc .* Sxy));
C += Sxyc .* Sy .* Sx;
end
ans = (A .* B ./ samples ./ samples, C./samples)
return ans
end
data = rand(24000,64);
limit = 2000;
ix = int64([1:limit/2]);
iy = ix[1:end/2];
mg = zeros(Int64,length(iy),length(ix));
mx = broadcast(+,ix',mg);
my = broadcast(+,iy,mg);
S = rfft(data,1)./24000;
@time (AB, C) = expensive_hat(S,mx,my);