Re: masked arrays and nd_image

On Wed, 2004-08-11 at 06:24, Todd Miller wrote:

> I think the key here may be the "filled()" method which lets you convert
> a masked array into a NumArray with the masked values filled into some
> fill value, say 0.  I'm not sure what the post convolution mask value
> should be.

I hope I'm not jumping in where I don't belong, but here at SFO/CSUN
we've had quite a bit of experience with convolutions and correlations
of time series (serieses?) with missing data points.

I'm sure you all know this, but:  For the scaling to be correct, you
have to not only mask out the values you don't want, but normalize the
sum to reflect the fact that different numbers of values will appear in
the sum.  Our MATLAB code to convolve x and y, with bad points marked by
NaNs, is:

for i = 1 : xlen-ylen+1
   j(i)=i;
   x1=x(i:i+ylen-1);
   a=x1.*y;
   b=a(find(~isnan(a)));
   if isempty(b)
      z(i)= NaN;
   else
       z(i)=sum(b)/length(b)
   end
end

I'd be happy to know how to code up the equivalent in numarray.  In the
above, note that x1 is x padded on both ends with ylen-1 NaNs.

Unfortunately, and again I'm sure everyone knows this, you can't use
FFTs to speed up convolutions/correlations if you have missing data
points, so you have to use discrete techniques.  The numerical analysis
literature refers to this problem as Fourier analysis of unequally
spaced data.  The only publications and algorithms I could find went the
wrong way:  given an unequally spaced set of points in Fourier space,
find the most likely reconstruction in real space.

-- 
Stephen Walton <stephen.walton@xxxxxxxx>
Dept. of Physics & Astronomy, Cal State Northridge

signature.asc
Description: This is a digitally signed message part