MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# Thread Subject: Confusion over calculating Fourier transforms

Subject: Confusion over calculating Fourier transforms

From: Mat Hunt

### Mat Hunt

Date: 2 Dec, 2009 09:30:21

Message: 1 of 12

Hi,

I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):

>x=-3:0.001:3;
>y=exp(-x.^2);
>f=abs(fft(y));
>plot(f);

I just set something which is completely wrong, what am I doing wrong here?

Mat

Subject: Confusion over calculating Fourier transforms

From: TideMan

### TideMan

Date: 2 Dec, 2009 10:08:04

Message: 2 of 12

On Dec 2, 10:30 pm, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Hi,
>
> I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):
>
> >x=-3:0.001:3;
> >y=exp(-x.^2);
> >f=abs(fft(y));
> >plot(f);
>
> I just set something which is completely wrong, what am I doing wrong here?
>
> Mat

Did you read the response when you posted this exact same question 20
hours ago?

Did you try what Matt J suggested?
And why did you start a new thread, rather than continue with the old
one?

 Subject: Confusion over calculating Fourier transforms From: Rune Allnor Date: 2 Dec, 2009 10:27:20 Message: 3 of 12 On 2 Des, 10:30, "Mat Hunt" wrote: > Hi, > > I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say): > > >x=-3:0.001:3; > >y=exp(-x.^2); > >f=abs(fft(y)); > >plot(f); > > I just set something which is completely wrong, what am I doing wrong here? The computed result is correct. The FFT function might not compute what you think it does, though. Read up on the various flavours of the Fourier transform in general (there are at least four of them), and the DFT in particular. Rune

Subject: Confusion over calculating Fourier transforms

From: Mat Hunt

### Mat Hunt

Date: 2 Dec, 2009 14:02:20

Message: 4 of 12

Hi,

I couldn't find the previous post, my apologies.

To give you an idea where my problem comes from. I am calculating the shape of a free surface of a train of waves, i wrote it as a fourier transform, \eta =\int_{\mathbb{R}}A(k)\exp (ikx)dk. I have calculated the function A(k) which is rather nasty, so i was thinking about using matlab to compute it numerically.

I was using the Gaussian to get the feel of things with matlab.

Mat

Subject: Confusion over calculating Fourier transforms

From: Matt J

### Matt J

Date: 2 Dec, 2009 15:51:03

Message: 5 of 12

"Mat Hunt" <hunt_mat@hotmail.com> wrote in message <hf5c3c$1f5$1@fred.mathworks.com>...
> Hi,
>
> I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):
>
> >x=-3:0.001:3;
> >y=exp(-x.^2);
> >f=abs(fft(y));
> >plot(f);
>
> I just set something which is completely wrong, what am I doing wrong here?

You need to do plot(fftshift(f)) to see the Gaussian as a centered peak (it will be very narrow, so you'll need to zoom in). Right now, you're seeing a circulant shift of this peak.

Subject: Confusion over calculating Fourier transforms

From: Mat Hunt

### Mat Hunt

Date: 2 Dec, 2009 16:57:01

Message: 6 of 12

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hf62d6$k7a$1@fred.mathworks.com>...
> "Mat Hunt" <hunt_mat@hotmail.com> wrote in message <hf5c3c$1f5$1@fred.mathworks.com>...
> > Hi,
> >
> > I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):
> >
> > >x=-3:0.001:3;
> > >y=exp(-x.^2);
> > >f=abs(fft(y));
> > >plot(f);
> >
> > I just set something which is completely wrong, what am I doing wrong here?
>
> You need to do plot(fftshift(f)) to see the Gaussian as a centered peak (it will be very narrow, so you'll need to zoom in). Right now, you're seeing a circulant shift of this peak.

Okay, I have done this BUT it still doesn't account for getting the amplitude incorrect, the height should be less than 10 but the graph shows that it is in the thousands. Something has gone wrong here and I am not too sure what. As I said before, I want to plot the transform \hat{f}(k) as a function of k. So I need to find k also.

Subject: Confusion over calculating Fourier transforms

From: Matt J

### Matt J

Date: 2 Dec, 2009 18:12:03

Message: 7 of 12

"Mat Hunt" <hunt_mat@hotmail.com> wrote in message <hf668t$pak$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <hf62d6$k7a$1@fred.mathworks.com>...
> > "Mat Hunt" <hunt_mat@hotmail.com> wrote in message <hf5c3c$1f5$1@fred.mathworks.com>...
> > > Hi,
> > >
> > > I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):
> > >
> > > >x=-3:0.001:3;
> > > >y=exp(-x.^2);
> > > >f=abs(fft(y));
> > > >plot(f);
> > >
> > > I just set something which is completely wrong, what am I doing wrong here?
> >
> > You need to do plot(fftshift(f)) to see the Gaussian as a centered peak (it will be very narrow, so you'll need to zoom in). Right now, you're seeing a circulant shift of this peak.
>
> Okay, I have done this BUT it still doesn't account for getting the amplitude incorrect, the height should be less than 10 but the graph shows that it is in the thousands. Something has gone wrong here and I am not too sure what. As I said before, I want to plot the transform \hat{f}(k) as a function of k. So I need to find k also.
=====================

As I mention in the duplicate post, the height is not what you expect because the fft is a sum whereas the continuous Fourier Transform is an integral. To approximate integrals, discrete sums must be multiplied by sampling intervals (giving you a Riemman sum).

You also need to be aware of what formula MATLAB uses for the fft (see "help fft"). Additional scale factors of 1/sqrt(pi) might be necessary according to whether you follow the Fourier transform conventions of physicists or of engineers.

Subject: Confusion over calculating Fourier transforms

Date: 3 Dec, 2009 02:02:04

Message: 8 of 12

Right, just like Matt J says.

If you look at the formula, you will realize that your max in the spectrum [which occurs at DC] is given mathematically by:

maxSpectrum = X(k=0) = sum(exp(-x.^2))

where x is just like you have defined at the very beginning [-3:0.001:3].

Look at this sum and compare it to

max(f)

You will see what I mean.

Subject: Confusion over calculating Fourier transforms

From: Mat Hunt

### Mat Hunt

Date: 3 Dec, 2009 11:28:02

Message: 9 of 12

It's been a while since I last did Fourier analysis at university as a pure maths course and coming back now is a bit hard. I define the transform as:

\hat{f}(\xi)=\int_{-\infty}^{\infty}f(x)e^{-i\xi x}dx

Does Matlab define it differently? if so what? I have had a look as all it gives is the discrete Riemann sum and I am not too sure about what they're doing.

I searched online and I came across the following example:

N=128;
t=linspace(0,3,N);
f=2*exp(-3t)
Ts=t(2)-t(1);
Ws=2*pi/Ts;
F=fft(f);
Fc=fftshift(F)*Ts;

W=Ws*(-N/2:(N/2)-1)/N;
Fa=2./(3+W*i);
plot(W,abs(Fa),W,abs(Fc),'.');
ylabel('|F(\omega)|');
title('Fourier Transform Approximation');

This worked quite well, but when I changed f to f=exp(-t.^2) and likewise the Fourier transform Fa to Fa=sqrt(pi)*exp(-0.25*W.^2), the comparison is very very bad.

Again it seem very odd to me.

Mat

Subject: Confusion over calculating Fourier transforms

Date: 3 Dec, 2009 11:35:24

Message: 10 of 12

The formula is here:

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.shtml

You should substitute k = 1 in X(k).

Subject: Confusion over calculating Fourier transforms

From: Greg Heath

### Greg Heath

Date: 3 Dec, 2009 11:51:03

Message: 11 of 12

On Dec 2, 4:30 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Hi,
>
> I need to calculate the Fouier transform of some rather complicated functions. Naively, I would do the following (to compute the Fourier transform of a Gaussian say):
>
> >x=-3:0.001:3;
> >y=exp(-x.^2);
> >f=abs(fft(y));
> >plot(f);
>
> I just set something which is completely wrong, what am I doing wrong here?

1. Matlab defines the FFT over the nonegative domain

x = dx*(0:N-1) = 0:dx:X-dx = linspace(0,X-dx,N)

where X = N*dx.

2. The proper vertical scaling is given by

Y = fft(y)/N;

The best way to verify this is via a sinusoid.

3. An alternate to using ifftshift(x) is to center y via

y = exp(-(x-3).^2) for 0 <= x <=6

Hope this helps.

Greg

Subject: Confusion over calculating Fourier transforms

From: Matt J

### Matt J

Date: 3 Dec, 2009 16:20:07

Message: 12 of 12

"Mat Hunt" <hunt_mat@hotmail.com> wrote in message <hf87c2$3gp$1@fred.mathworks.com>...
> It's been a while since I last did Fourier analysis at university as a pure maths course and coming back now is a bit hard. I define the transform as:
>
> \hat{f}(\xi)=\int_{-\infty}^{\infty}f(x)e^{-i\xi x}dx
>

Then you need be scaling the output of the fft by dx. Also, you need to plot the result on the following axis

Axis= -ceil((N-1)/2):N-1-ceil((N-1)/2)*2*pi/(N*dx)

To better understand wha't going on, it would be a helpful exercise for you to discretize your above integral by making the substitutions

xi=k*(2*pi/(N*dx))
x=n*dx

where k and n are discrete variables ranging over

Axis= -ceil((N-1)/2):N-1-ceil((N-1)/2)

You should find that the integral discretizes to the following MATLAB operations:

\hat{F}=ifftshift( fft(fftshift(F) ) )*dx

where

F=f(Axis*dx)

\hat{F}=\hat{f}(2*pi*Axis/(N*dx))