FPGA Central

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" <[email protected]> wrote in message
news:[email protected]..
> 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" <[email protected]> wrote in message news:<cfp1n4$bp4$[email protected]>...
<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" <[email protected]> wrote in message
> news:[email protected]..
>
>>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"@[email protected]> wrote in message
news:[email protected]..
>
>
> 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" <[email protected]> wrote in message
news:[email protected] om...
> "lucy" <[email protected]> wrote in message
news:<cfp1n4$bp4$[email protected]>...
> <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" <[email protected]> wrote in message news:<cfqusi$mif$[email protected]>...
> "Derek Goring" <[email protected]> wrote in message
> news:[email protected] om...
> > "lucy" <[email protected]> wrote in message
> news:<cfp1n4$bp4$[email protected]>...
> > <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?

"Derek Goring" <[email protected]> wrote in message
news:[email protected] om...
> "lucy" <[email protected]> wrote in message
news:<cfqusi$mif$[email protected]>...
> > "Derek Goring" <[email protected]> wrote in message
> > news:[email protected] om...
> > > "lucy" <[email protected]> wrote in message
> > news:<cfp1n4$bp4$[email protected]>...
> > > <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?

FYI, The conv2 in Matlab is 10 times slower than if I do
ifft2(fft2(A).*fft2(B))),

but the latter is 100 times slower than first do 1D filter on rows then do
1D filter on columns,...


> But your function is continuous, so why
> not do the Fourier transform analytically and discretise in the
> frequency domain?

This is really a good suggestion, suppose I do everything Fourier Transform
analytically and the overall computation result is discreteized in the
frequency domain,

how do I convert the frequency domain samples back to time-domain to get my
time-domain results?

Thanks a lot,

-Lucy

lucy wrote:

> "Derek Goring" <[email protected]> wrote in message
> news:[email protected] om...
>
>>"lucy" <[email protected]> wrote in message
>
> news:<cfqusi$mif$[email protected]>...
>
>>>"Derek Goring" <[email protected]> wrote in message
>>>news:[email protected] e.com...
>>>
>>>>"lucy" <[email protected]> wrote in message
>>>
>>> news:<cfp1n4$bp4$[email protected]>...
>>>
>>>><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?
>
>
> FYI, The conv2 in Matlab is 10 times slower than if I do
> ifft2(fft2(A).*fft2(B))),
>
> but the latter is 100 times slower than first do 1D filter on rows then do
> 1D filter on columns,...
>
>
>
>>But your function is continuous, so why
>>not do the Fourier transform analytically and discretise in the
>>frequency domain?
>
>
> This is really a good suggestion, suppose I do everything Fourier Transform
> analytically and the overall computation result is discreteized in the
> frequency domain,
>
> how do I convert the frequency domain samples back to time-domain to get my
> time-domain results?
>
> Thanks a lot,
>
> -Lucy
>
>
Inverse FFT. You should be able to do an FFT on your data, then do an
IFFT on the result, and get your data back. Or do an FFT on your data,
filter it in the frequency domain (by multiplication) and IFFT back.
This works just as well for the 2D FFT as for 1D.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com