DSPRelated.com
Forums

how to decide a good sampling rate for sampling a function without obvious frequency?

Started by lucy August 15, 2004
Hi all,

I am having trouble with my sampling problem:

Basically I want to discretize a 2D Gaussian function f=gaussian(x, y)

gaussian(x, y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2));

In my experiments using Matlab, I am using square grids to deiscretize the
above function.

Suppose I have a grid -- [-N..N, -N..N] where N is the number of samples in
one axis. So all together there are (2N+1)^2 samples.

Suppose the stepsize for discretization is dx.

The support of the discretized function will be [-N*dx..N*dx, -N*dy..N*dy].
Since Gaussian function has an infinite support itself, I have to make N*dx
and N*dy sufficiently large in order to capture most of the energy of the
function f(x, y).

I remember for Gaussian distribution, it should be at least 3*sigma.

So I decide that the following conditions must be satisfied:

N*dx >= 4*sigma_x;
N*dy >= 4*sigma_y;

(sigma_x and sigma_y are known).

If the stepsize is large, the discretization will not be accurate. So dx and
dy should be small compared with sigma_x and sigma_y...

For example, I can choose dx=sigma_x / 100, dy= sigma_y /100,

but then N should be >= 400.

But then

[-N..N, -N..N] grids will have (2N+1)^2 = 801^2=641601 samples...

This is too large and will be very bad for subsequent operations...

After many experiments, I want to ask how to determine a good stepsize
taking into consideration of the storage size, while maintain good accuracy?

I know Nyquist sampling theory, but after many thinking, I found nowhere
here I can apply that theory...

Please help me and give me a hand!

Thank you very much

-Lucy










lucy wrote:

> Hi all, > > I am having trouble with my sampling problem: > > Basically I want to discretize a 2D Gaussian function f=gaussian(x, y) > > gaussian(x, y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); > > In my experiments using Matlab, I am using square grids to deiscretize the > above function. > > Suppose I have a grid -- [-N..N, -N..N] where N is the number of samples in > one axis. So all together there are (2N+1)^2 samples. > > Suppose the stepsize for discretization is dx. > > The support of the discretized function will be [-N*dx..N*dx, -N*dy..N*dy]. > Since Gaussian function has an infinite support itself, I have to make N*dx > and N*dy sufficiently large in order to capture most of the energy of the > function f(x, y). > > I remember for Gaussian distribution, it should be at least 3*sigma. > > So I decide that the following conditions must be satisfied: > > N*dx >= 4*sigma_x; > N*dy >= 4*sigma_y; > > (sigma_x and sigma_y are known). > > If the stepsize is large, the discretization will not be accurate. So dx and > dy should be small compared with sigma_x and sigma_y... > > For example, I can choose dx=sigma_x / 100, dy= sigma_y /100, > > but then N should be >= 400. > > But then > > [-N..N, -N..N] grids will have (2N+1)^2 = 801^2=641601 samples... > > This is too large and will be very bad for subsequent operations... > > After many experiments, I want to ask how to determine a good stepsize > taking into consideration of the storage size, while maintain good accuracy? > > I know Nyquist sampling theory, but after many thinking, I found nowhere > here I can apply that theory... > > Please help me and give me a hand! > > Thank you very much > > -Lucy >
Use the Nyquist sampling theory on the frequency _components_ of the signal. Decide how much inaccuracy you can stand, and at what frequencies. Now look at the power spectral density of the signal you want to sample. Due to aliasing any signal energy above the Nyquist frequency will be "folded back" into your real signal, corrupting it. Because you know how much inaccuracy you can stand you can locate the Nyquist frequency to give you no more than your desired accuracy limit. Disclaimer: That made perfect sense coming out, but on re-reading it I think it would be a good growth enhancer in the garden. None the less I'm going to hit the send button... -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
"Tim Wescott" <tim@wescottnospamdesign.com> wrote in message
news:10i0aaprl2va25b@corp.supernews.com...
> lucy wrote: > > > Hi all, > > > > I am having trouble with my sampling problem: > > > > Basically I want to discretize a 2D Gaussian function f=gaussian(x, y) > > > > gaussian(x,
y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2));
> > > > In my experiments using Matlab, I am using square grids to deiscretize
the
> > above function. > > > > Suppose I have a grid -- [-N..N, -N..N] where N is the number of samples
in
> > one axis. So all together there are (2N+1)^2 samples. > > > > Suppose the stepsize for discretization is dx. > > > > The support of the discretized function will be
[-N*dx..N*dx, -N*dy..N*dy].
> > Since Gaussian function has an infinite support itself, I have to make
N*dx
> > and N*dy sufficiently large in order to capture most of the energy of
the
> > function f(x, y). > > > > I remember for Gaussian distribution, it should be at least 3*sigma. > > > > So I decide that the following conditions must be satisfied: > > > > N*dx >= 4*sigma_x; > > N*dy >= 4*sigma_y; > > > > (sigma_x and sigma_y are known). > > > > If the stepsize is large, the discretization will not be accurate. So dx
and
> > dy should be small compared with sigma_x and sigma_y... > > > > For example, I can choose dx=sigma_x / 100, dy= sigma_y /100, > > > > but then N should be >= 400. > > > > But then > > > > [-N..N, -N..N] grids will have (2N+1)^2 = 801^2=641601 samples... > > > > This is too large and will be very bad for subsequent operations... > > > > After many experiments, I want to ask how to determine a good stepsize > > taking into consideration of the storage size, while maintain good
accuracy?
> > > > I know Nyquist sampling theory, but after many thinking, I found nowhere > > here I can apply that theory... > > > > Please help me and give me a hand! > > > > Thank you very much > > > > -Lucy > > > > Use the Nyquist sampling theory on the frequency _components_ of the > signal. Decide how much inaccuracy you can stand, and at what > frequencies. Now look at the power spectral density of the signal you > want to sample. Due to aliasing any signal energy above the Nyquist > frequency will be "folded back" into your real signal, corrupting it. > Because you know how much inaccuracy you can stand you can locate the > Nyquist frequency to give you no more than your desired accuracy limit. > > Disclaimer: That made perfect sense coming out, but on re-reading it I > think it would be a good growth enhancer in the garden. None the less > I'm going to hit the send button... > > -- > > Tim Wescott > Wescott Design Services > http://www.wescottdesign.com
Hi Tim, Thank you very much for your help. Yet your hint reads very abstract. For the 2D gaussian function: gaussian(x, y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); I've computed its FT to be GAUSSIAN(u, v)=exp(-0.5*(u^2/sigmau^2+v^2/sigmav^2)); where sigma_u=1/(2*pi*sigma_x), sigma_v=1/(2*pi*sigma_y)... In my numerical calculation, I have chosen sigma_x = 1 and sigma_y... So sigma_u=1/(2*pi), sigma_v=1/(2*pi)... Since the frequency components are still in Gaussian shape, if I want to 99% of the freuqency components, I should keep the frequency ranging between [-4*sigma_u .. 4*sigma_u, -4*sigma_v .. 4*sigma_v] If I pick 4*sigma_u = 4*sigma_v to be my Nyquist frequency, and double it, my sampling rate should be 8*sigma_u=8*sigma_v = 8/(2*pi)=4/pi=1.27, so I pick 2. So am I correct that I should only need to sample the gaussian function twice? Thank you very much, -Lucy
"lucy" <losemind@yahoo.com> wrote in message news:<cfp1n4$bp4$1@news.Stanford.EDU>...
<snip>
>Lucy is trying to fit a round peg into a square hole
> gaussian(x, y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); > > In my experiments using Matlab, I am using square grids to deiscretize the > above function.
This is your problem right here. Why are you using square grids? What's so special about them, especially since what you're discretising is naturally a cylindrical polar shape? As well as that, you need greater resolution near the centre and less at the edges, so your cylindrical polar grid needs to have radial spacing that is a function of the radius. Perhaps you could tell us what your application is? What will you do with the Gaussian, once you've discretised it?

lucy wrote:

>big snip<
> > So am I correct that I should only need to sample the gaussian function > twice? >
Well, yes you should be able to reconstruct an approximation of a gaussian with 2 sinusoids - 4 would be better - more even better. I think as you decrease the sample spacing geometrically the resolution will increase exponentially so you will not need a very large number of samples. Just make your sampling algorithm scalable. Start with a small managable size and if you need better resolution later you can do that. -jim -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
lucy wrote:

> "Tim Wescott" <tim@wescottnospamdesign.com> wrote in message > news:10i0aaprl2va25b@corp.supernews.com... > >>lucy wrote: >> >> >>>Hi all, >>> >>>I am having trouble with my sampling problem: >>> >>>Basically I want to discretize a 2D Gaussian function f=gaussian(x, y) >>> >>>gaussian(x, > > y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); > >>>In my experiments using Matlab, I am using square grids to deiscretize > > the > >>>above function. >>>
-- snip --
>>> >>>-Lucy >>> >> >>Use the Nyquist sampling theory on the frequency _components_ of the >>signal. Decide how much inaccuracy you can stand, and at what >>frequencies. Now look at the power spectral density of the signal you >>want to sample. Due to aliasing any signal energy above the Nyquist >>frequency will be "folded back" into your real signal, corrupting it. >>Because you know how much inaccuracy you can stand you can locate the >>Nyquist frequency to give you no more than your desired accuracy limit. >> >>Disclaimer: That made perfect sense coming out, but on re-reading it I >>think it would be a good growth enhancer in the garden. None the less >>I'm going to hit the send button... >> >>-- >> >>Tim Wescott >>Wescott Design Services >>http://www.wescottdesign.com > > > > Hi Tim, > > Thank you very much for your help. Yet your hint reads very abstract.
I know; I just didn't have time last night to go over it and figure out why, or how to make it more clear. Some time spent with graph paper and a pencil will be required to really understand the subject, no matter how well I explain it.
> > For the 2D gaussian function: > > gaussian(x, y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); > > I've computed its FT to be > > GAUSSIAN(u, v)=exp(-0.5*(u^2/sigmau^2+v^2/sigmav^2)); > > where sigma_u=1/(2*pi*sigma_x), sigma_v=1/(2*pi*sigma_y)... > > In my numerical calculation, I have chosen sigma_x = 1 and sigma_y... > > So sigma_u=1/(2*pi), sigma_v=1/(2*pi)... > > Since the frequency components are still in Gaussian shape, if I want to 99% > of the freuqency components, I should keep the frequency ranging between > [-4*sigma_u .. 4*sigma_u, -4*sigma_v .. 4*sigma_v] > > If I pick 4*sigma_u = 4*sigma_v to be my Nyquist frequency, > > and double it, my sampling rate should be 8*sigma_u=8*sigma_v = > 8/(2*pi)=4/pi=1.27, so I pick 2. > > So am I correct that I should only need to sample the gaussian function > twice?
This gets back to how much accuracy you need. If you're truncating the function at 3 sigma then it certainly indicates that sampling it twice will be adequate. On the other hand that seems a terribly small sample to me. Are you sure that your 4 sigma/99% assumption is true in 2-D, and are you sure that truncating at 3 sigma is appropriate? If your intuition is telling you that two samples is too small then perhaps the error assumptions that you started with are too generous, and you need to tighten them up? Getting 99% of all the energy sounds good, but it only puts your Nyquist energy down 20 decibels -- a sampling system wouldn't be considered "high quality" unless it were getting the aliasing down to at least 40 to 60dB or more, which implies getting 99.9% or even 99.99% of your energy. If I did this calculation from the precept I stated above and came up with these numbers I would do one of two things: One, I would try it with the indicated number of samples, and 2x that, and perhaps 10x that. Two, I would go over my math to make sure that I really knew what I was talking about.
> > Thank you very much, > > -Lucy > >
-- Tim Wescott Wescott Design Services http://www.wescottdesign.com
"jim" <"N0sp"@m.sjedging@mwt.net> wrote in message
news:4120a820_1@corp.newsgroups.com...
> > > lucy wrote: > > >big snip< > > > > > So am I correct that I should only need to sample the gaussian function > > twice? > > > > Well, yes you should be able to reconstruct an approximation of a > gaussian with 2 sinusoids - 4 would be better - more even better. I > think as you decrease the sample spacing geometrically the resolution > will increase exponentially so you will not need a very large number of > samples. Just make your sampling algorithm scalable. Start with a small > managable size and if you need better resolution later you can do that. > > -jim > > > -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- > http://www.newsfeeds.com - The #1 Newsgroup Service in the World! > -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
What kind of reconstruction filter do you use if you use only 2 samples of the Gaussian filter to approximate it? I did not do interpolation in resconstruction. Maybe I should do that?
"Derek Goring" <nztideman@yahoo.co.nz> wrote in message
news:d94b83b8.0408152329.2b13a3fb@posting.google.com...
> "lucy" <losemind@yahoo.com> wrote in message
news:<cfp1n4$bp4$1@news.Stanford.EDU>...
> <snip> > >Lucy is trying to fit a round peg into a square hole > > > gaussian(x,
y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2));
> > > > In my experiments using Matlab, I am using square grids to deiscretize
the
> > above function. > > This is your problem right here. > Why are you using square grids? > What's so special about them, especially since what you're > discretising is naturally a cylindrical polar shape? > As well as that, you need greater resolution near the centre and less > at the edges, so your cylindrical polar grid needs to have radial > spacing that is a function of the radius. > > Perhaps you could tell us what your application is? > What will you do with the Gaussian, once you've discretised it?
I need to use these sampled grids for 2D convolution... if I use a non-uniform sampling, how can I handle the subsequent 2D convolution operation? I am intested in learning that if the non-uniform sampling can handle convolution after the discretization? Thank you very much! -Lucy

lucy wrote:

> What kind of reconstruction filter do you use if you use only 2 samples of > the Gaussian filter to approximate it? >
The frequency response (mag) of the filter [1 1] is a raised cosine that's a reasonable approximation of a gaussian. Adding more coefficients will give more accuracy. You don't need very many to have a fairly accurate representation. -jim -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
"lucy" <losemind@yahoo.com> wrote in message news:<cfqusi$mif$1@news.Stanford.EDU>...
> "Derek Goring" <nztideman@yahoo.co.nz> wrote in message > news:d94b83b8.0408152329.2b13a3fb@posting.google.com... > > "lucy" <losemind@yahoo.com> wrote in message > news:<cfp1n4$bp4$1@news.Stanford.EDU>... > > <snip> > > >Lucy is trying to fit a round peg into a square hole > > > > gaussian(x, > y)=1/(2*pi*sigmax*sigmay)*exp(-0.5*(x^2/sigmax^2+y^2/sigmay^2)); > > > > > > In my experiments using Matlab, I am using square grids to deiscretize > the > > > above function. > > > > This is your problem right here. > > Why are you using square grids? > > What's so special about them, especially since what you're > > discretising is naturally a cylindrical polar shape? > > As well as that, you need greater resolution near the centre and less > > at the edges, so your cylindrical polar grid needs to have radial > > spacing that is a function of the radius. > > > > Perhaps you could tell us what your application is? > > What will you do with the Gaussian, once you've discretised it? > > I need to use these sampled grids for 2D convolution... if I use a > non-uniform sampling, how can I handle the subsequent 2D convolution > operation? > > I am intested in learning that if the non-uniform sampling can handle > convolution after the discretization? > > Thank you very much! > > -Lucy
So, once you've discretised, you'll transform to the frequency domain to do the convolution, right? But your function is continuous, so why not do the Fourier transform analytically and discretise in the frequency domain?