Calculating Characteristic function from PDF

Hello all,

I am trying to calculate the PDF from characteristic function 
and vice-versa for Gaussian pdf. It is well learnt that they 
form Fourier Transform pairs.

My first shot was to calculate Characteristic function from 
Gaussian PDF. I used two methods to calculate but invain.

First method I diretly used the matlab FFT() function, I get peaks
instead of  expected bell shape.

Second method, I tried calculating modified FFT by replacing 
n with x(n) as x-axis is now -4 to +4, but no success.
Following method is used;

I generate Gaussian-pdf in range -4 to +4
x = -4:0.001:4;
gaussian_pdf = exp(-(x.^2)/2)/(sqrt(2*pi));

I use the following function (to calculate Fourier Transform):

for k = 1:N,
   for n = 1:N,
      temp = temp + gaussian_pdf(n)*exp(-j*2*pi*(k-1)*x(n)/(max(x)-min(x)));
   end
   X(k) = temp/N;
   temp = 0;
end


Please comment,if I am doing something wrong here (which I am!!).
Seeking any suggestion.

Regards,
Vimal

Vimal wrote:

> Please comment,if I am doing something wrong here (which I am!!).
> Seeking any suggestion.

How about calculating it with pencil & paper?

It seems to me it's much easier.

bye,

-- 
   Piergiorgio Sartor

Hello Piergiorgio,

Maybe Vimal doesn't know the tricks for finding the FT of a gaussian
function.

Vimal, you will need to complete the square within the argument of the
exponent. Also you will need to know how to find integral(over all space) of
exp(-x^2) dx. You will have other scale factors, but I assume you can do
u-substitution.  This integral is commonly found (over than table look up)
by squaring the integral, writing as an iterated integral, and then
converting to polar coordinates. When you get your final answer, don't
forget to take the square root of the data.

If you do all of this you will be verifying that a gaussian function is an
identity function under fourier transformation.

Clay





"Piergiorgio Sartor" <piergiorgio.sartor@nexgo.REMOVE.THIS.de> wrote in
message news:3f3a40bf$0$258$4d4ebb8e@read.news.de.uu.net...
> Vimal wrote:
>
> > Please comment,if I am doing something wrong here (which I am!!).
> > Seeking any suggestion.
>
> How about calculating it with pencil & paper?
>
> It seems to me it's much easier.
>
> bye,
>
> -- 
>    Piergiorgio Sartor
>

Clay S. Turner wrote:

> Vimal, you will need to complete the square within the argument of the
> exponent. Also you will need to know how to find integral(over all space) of
> exp(-x^2) dx. You will have other scale factors, but I assume you can do
> u-substitution.  This integral is commonly found (over than table look up)
> by squaring the integral, writing as an iterated integral, and then
> converting to polar coordinates. When you get your final answer, don't
> forget to take the square root of the data.

I think it's much more simple.

The point is that the integral of a PDF is always 1.

So, it's not necessary to calculate the integral at
all, it's only necessary to shuffle the terms under it
to get back a PDF multiplied by some constant...

Since the characteristic function is something like
E[e^jux] and since the Gaussian PDF contains already
some exponential, the trick should be quite easy, I guess.

OK, I spoil it, the completition of a square of an
exponential is required...

bye,

-- 

piergiorgio

vimal_bhatia2@yahoo.com (Vimal) wrote in message news:<b6fc6dda.0308130525.979c887@posting.google.com>...
> Hello all,
> 
> I am trying to calculate the PDF from characteristic function 
> and vice-versa for Gaussian pdf. It is well learnt that they 
> form Fourier Transform pairs.
> 
> My first shot was to calculate Characteristic function from 
> Gaussian PDF. I used two methods to calculate but invain.
> 
> First method I diretly used the matlab FFT() function, I get peaks
> instead of  expected bell shape.

Try the following

x = -4:0.001:4;
gaussian_pdf = exp(-(x.^2)/(2*0.000004))/(sqrt(2*pi)*(0.002));
plot(abs(fft(gaussian_pdf)))

You will understand why you were thinking you were doing something wrong.

> 
> Second method, I tried calculating modified FFT by replacing 
> n with x(n) as x-axis is now -4 to +4, but no success.
> Following method is used;
> 
> I generate Gaussian-pdf in range -4 to +4
> x = -4:0.001:4;
> gaussian_pdf = exp(-(x.^2)/2)/(sqrt(2*pi));
>

> Try the following
> 
> x = -4:0.001:4;
> gaussian_pdf = exp(-(x.^2)/(2*0.000004))/(sqrt(2*pi)*(0.002));
> plot(abs(fft(gaussian_pdf)))
> 
> You will understand why you were thinking you were doing something wrong.
> 

Thanks a lot for the replies. I could see it hapening in theory but
when it comes to matlab things doesn't always translate.

Working last evening I could could figure out using zero-padding or
extending the pdf far greater than 4, gave me the plot/smooth Gaussian
pdf. I beleive that is done in above MATLAB code as well.

Thanks again folks!!

Dear Clay, 

Can you explain a bit or send some pionters where "complete the
square", and u-substitution explained.

thanks,

Dear Piergiorgio, 

As the integral of PDF if "1", if I normalize any set of data, then
can I use the same trick for virtually any integration. I guess not...
But I guess exponential's can be dealt like this.

cheers

"Vimal" <vimal_bhatia2@yahoo.com> wrote in message
news:b6fc6dda.0308140519.7c470e8b@posting.google.com...
>
> Dear Clay,
>
> Can you explain a bit or send some pionters where "complete the
> square", and u-substitution explained.
>
> thanks,
>

Hello Vimal,

"Completing the square" is a trick that exploits the relation

(x+a)^2 = x^2 + 2ax + a^2

Here you notice that when the 3rd term's coefficient (RHS) equals the square
of one half of the second term (RHS), then the three terms can be collapsed
into a perfect square. (The coef of the x^2 must be unity - if not then just
factor it out)

Example, let's say you are given

5  =  x^2 + 4x

And you want to find x, so just add the square of one half of the middle
term's coef to both sides.  In this case you will be adding 4.

5+4 = x^2 + 4x + 4

or

9 = (x+2)^2

+-3 =  x+2

 1 or -5 = x

Normally one would use the general quadratic formula for these problems. But
if you look up how the formula is derived, you will see that the "completing
of the square" trick is used.

u-substitution is a technique of integration where the integrand is
transformed into something simpler. Namely where the integral fits the
following form:

            (        du  )
integral ( f(u) ---- ) dx    = integral (  f(u) du )
            (        dx  )

Example

Find integral sin^2(x) cos(x) dx

Let u = sin(x) and hence  du/dx is cos(x)

So we are now just integrating  u^2

which results in     u^3
                            --- + c
                             3

and after replacing  u with sin(x), we find:


sin^3(x)
-------- + c
     3


I hope this helps

Clay

Vimal wrote:

> As the integral of PDF if "1", if I normalize any set of data, then
> can I use the same trick for virtually any integration. I guess not...
> But I guess exponential's can be dealt like this.

In this specific Gaussian case, you can revert
the integral into a Gaussian one (and a scale
factor) due to the nature of Gaussian PDF.

So integration does not really need to be
carried out.

You have: Phi = E[e^sx] or similar.

PDF = e^(x^2)*K

So, the integral is:

sum(e^sx * e^(x^2) * K)

and

sum(e^(sx+x^2) * K)

or something like that, check the terms and conditions.

With some magic this can be converted back to a
Gaussian PDF multiplied by something else.

sum(K * e^(something^2) * e^(s^2))

or similar, the result is than obvious, I hope.

bye,

-- 
   Piergiorgio Sartor

"Clay S. Turner" <physicsNOOOOSPPPPAMMMM@bellsouth.net> wrote in message
news:gKM_a.815$sx4.530@fe05.atl2.webusenet.com...
>
> (x+a)^2 = x^2 + 2ax + a^2
>
> Here you notice that when the 3rd term's coefficient (RHS) equals the
square
> of one half of the second term (RHS), then the three terms can be
collapsed
> into a perfect square. (The coef of the x^2 must be unity - if not then
just
> factor it out)

There's nothing like a good transposition.  I meant to say "when the 3rd
term equals the square of one half of the middle term's coefficient".

The 3rd term is "a^2" and the middle term's coef is "2a."

Sorry about that.

Clay

"Clay S. Turner" wrote:
> 
  ...
> 
> Normally one would use the general quadratic formula for these problems. But
> if you look up how the formula is derived, you will see that the "completing
> of the square" trick is used.
> 
  ...

Clay,

That derivation is flawed because it requires the hand waving that sqrt
is a two-valued function (+ and -), when sqrt() is specifically defined
to be positive. A function must be single valued. Moreover, a scale
factor is lost in the usual quadratic formula. (Ask Rick Lyons.) Here's
my derivation. The scale factor is retained because the quadratic isn't
equated to zero. Moreover, I need no special pleading that sqrt() can
sometimes be negative.

Factor the expression ax� + bx + c.
a{x� + (b/a)x + c/a}
a{x� + (b/a)x + (b�/4a�) -(b�/4a�) + c/a}
a{x� + (b/a)x + (b�/4a�) -(b� - 4ac)/4a�)}
a{(x + b/2a)� - [sqrt(b� - 4ac)/2a]�}; a difference of two squares, u,v.
a{u� - v�} factors to a{(u + v)(u - v)}, the three factors -- that's
right, three! -- that, when multiplied, yield the original expression.
Let me repeat (for emphasis, not to gloat). A quadratic expression has
two roots, but three factors. We forget the third at our peril.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������