Re: fft

[Mario Storti]

>>>>> On Thu, 29 Oct 1998 14:25:00 +0100, 
>>>>>      Diego Alberto Arias Prado <address@hidden> said:

>         Hello!
>         There is something in the fft() function that I don't understand. I
> compute the fft of a vector containing the samples of a sinc(), an exp() and
> a rectangular pulse and what I get is a vector of the same size. What I
> don't understand is this: in the positions 1, 2, 4.. and N-3, N-2, N-1 and N
> of the vector (its size is 1 N) I get "big" values and in the central
> positions the values (I mean, far from 1 and N) I get values near to cero.
>         Thus, from transforming a sinc() I get something that hasn't
> anything to see with a rectangular pulse (for example). HELP!
>         Thanks in advance.

Som care has to be taken wth the sample points in order to get the
expexted results with  the fft(). The assumption is that the function
is periodic so that point N, which is at (N-1)*L/N should get the
value of the function at -L/N. 

The simplest way that I found to get thigs work is to take the
coordinates of points as follows:

Dx = L/N

x_j = (j-1)*Dx      for j=1,....,N/2

x_j = (j-1)*Dx-L    for j=N/2+1,...,N

Note that if the function is periodic the shift in length L does not
matter. The way to code this in Octave is:

> n=2^m;
> L=2;
> x=L*(0:n-1)'/n;
> x=[x(1:n/2);x(n/2+1:n)-L];

Attached  below is a demo. 

Hope this helps.

Mario

%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%
Mario Alberto Storti                           | Fax: (54)(42) 55.09.44 |
Centro Internacional de Metodos Computacionales| Tel: (54)(42) 55.91.75 |
  en Ingenieria - CIMEC (INTEC/CONICET-UNL)    |........................|  
INTEC, Guemes 3450 - 3000 Santa Fe, Argentina                           |
Reply: address@hidden, http://venus.unl.edu.ar/gtm-eng.html |

============================================================
n=2^7;
L=2;
x=L*(0:n-1)'/n;
x=[x(1:n/2);x(n/2+1:n)-L];

for k=1:3
  if k==1
    y=sinc(10*x);
    Title="sinc function";
  elseif k==2
    y=abs(x)<.2;
    Title="square pulse";
  else
    y=exp(-(x/.15).^2);
    Title="Gaussian pulse";
  endif
  z=fft(y);
  title(Title);
  plot(x,y);
  disp("press enter")
  pause
  title(["Fourier transform of " Title])
  plot(x,[real(z) imag(z)]);
  pause
endfor
================================================================