On May 10, 6:07=A0pm, Tim Wescott <t...@seemywebsite.com> wrote:
> On 05/10/2011 07:24 AM, Andy365 wrote:
>
> > Hello,
> > I have irregular spaced samples in 3 dimensions: f(xi, yi, zi) that I
> > wish to take the FFT of.
> > Therefore I need to resample the data to a regular grid.
> > Since the data will be Fourier transformed; it will also be good to
> > smooth the data somewhat to avoid Gibbs noise ("ringing")?
>
> Gibbs 'noise' is not a consequence of data that is not smooth. =A0It is a
> consequence of data that is improperly windowed.
>
> > Therefore I am considering using a convolution kernel to resample the
> > data to a regular grid.
> > A convolution kernel could help me with both the resampling and the
> > smoothing in one step?
>
> In theory, yes. =A0In practice, I think you'll run into scaling problems.
> =A0 Resampling _regularly_ sampled data with a convolution kernel would
> give you a meaningful result, that was perhaps off by the ratio of the
> sampling rates involved. =A0Resampling _irregularly_ sampled data that wa=
y
> will give you a transform whose effective gain is low where your samples
> happen to be sparse, and high where they happen to be dense.
>
> > I am not sure which kernel to use, but read a bit about Lanczos
> > resampling on wikipedia and it seems promising.
> > According to wikipedia the kernel in 2 dimensions would be the product
> > of the kernel for each dimension:
> > L(x,y)=3DL(x)L(y),
> > where L(x)=3Dsinc(x)sinc(x/a) for -a<x<a, 0 otherwise.
> > Generalizing this to 3 coordinates I get:
> > L(x,y,z)=3DL(x)L(y)L(z)
>
> Correct, but as mentioned not applicable to your problem.
>
> > How about defining the kernel in radial coordinates instead:
> > L(r)=3Dsinc(r)sinc(r/a) for -a<r<a,
> > where r=3Dsqrt(x^2+y^2+z^2),
> > would that give the same result?
>
> Would it? =A0Does your L(r) equal your L(x, y, z)?
>
>
>
> > Also the formula for L(x) is given in normalized coordinates?
> > If so I think the formula with not normalized coordinates with
> > sampling in x =3D DX will be:
> > L(x)=3Dsinc(x/DX)sinc(x/(DX*a))
> > ?
>
> > According to wikipedia again; it is common to use a=3D2 or 3.
> > But if I want to achieve smoothing I might need to influence more than
> > 2 or 3 nodes in each direction?
> > Will things work with say a=3D10?
>
> As I said, I don't think that this approach will work well, if at all.
>
> You've gotten yourself lost in DSP land, and you're forgetting why
> you're doing the effort in the first place. =A0You want to extract some
> information from some noisy data. =A0Your starting question is "what is
> the best way to separate information from noise". =A0The answer to that
> question, in many cases, is to use filters. =A0But I don't think the
> answer is going to come from any combination of filters designed for
> consistently sampled data, because that's not what you have.
>
> So you need to back off a bit, and ask yourself "what's the best
> estimate of the value of this signal at a point where I don't have
> samples, given that I have samples of surrounding points?". =A0To get tha=
t
> answer you have to answer questions like "is my signal continuous", "how
> continuous is it*", "what is the nature of my noise", etc.
>
> * In continuous time or regularly sampled time, this would translate to
> "is this signal band limited?"
>
> --
>
> Tim Wescott
> Wescott Design Serviceshttp://www.wescottdesign.com
>
> Do you need to implement control loops in software?
> "Applied Control Theory for Embedded Systems" was written for you.
> See details athttp://www.wescottdesign.com/actfes/actfes.html
Thank you Tim for a very thorough reply!
My problem is really a bit more complex than what I have told so far
(sorry).
As an example say that I have n irregularly sampled points in space:
(xi, yi, zi), where each point represent that there is a microscopic
particle at that position.
I introduce the function g(x, y, z) that have the value 1: particle at
that position or 0: no particle at that position. I have n samples of
this function: g(xi, yi, zi) =3D 1.
Now I introduce the particle density: f(x,y,z): the sum of g(xi, yi,
zi) inside a small cell around (x,y,z).
My goal is to find the Fourier transform of this density: F(kx, ky,
kz).
In order to do that I need to resample the particle density to an
uniform grid.
Say I choose the following uniform grid: x =3D i, y =3D j and z =3D k, wher=
e
i, j and k are positive integers.
Now let us consider only one sample: g(1.6, 1.9, 1.9) =3D 1.
So far I have been adding this value (1) to the nearest neighbour
node: f(2, 2, 2): the density in the cell 1.5 <=3D x <=3D 2.5, 1.5 <=3D y <=
=3D
2.5 and 1.5 <=3D z <=3D 2.5.
However the problem with this is that this sample is also quite close
to the node f(1, 2, 2) but this cell does not get any contribution
from the sample :-(.
The result is a somewhat "blocky" density.
Instead I would like a more sophisticated method that divides the
value 1 between all surrounding nodes: f(1, 1, 1), f(1, 1, 2), f(1, 2,
1), f(1, 2, 2), f(2, 1, 1), f(2, 1, 2), f(2, 2, 1) and f(2, 2, 2).
f(2, 2, 2) is closest so this should get the largest share, say 0.6.
f(1, 2, 2) is also quite close so this get say 0.2, which leaves 0.2
to be divided among the remaining 6 nodes.
Therefore I was considering convolving with a kernel, that is looking
up from a window function e.g. Lanczos a weight (0 <=3D weight <=3D1)
based on how far the node is from the sample, and then adding each
weight to that node.
Andreas