direct and inverse fast Fourier transform

fft – fast Fourier transform
ifft – inverse fast Fourier transform

Calling sequence

y = fft(x, dim=mode)  
y = fft(x, mode)  
y = ifft(x, dim=mode)  
y = ifft(x, mode)



fft computes the discrete Fourier transform (dft) of a vector x m:

y = Fmx,   Fm  = [e−-m--]0≤k,j≤m− 1
and ifft computes the inverse transform, that is:
      −1      −1   1
y = Fm  x,  Fm  =  m-¯Fm
using fast Fourier transform algorithms. In most cases nsp uses the routines from the fftw3 library (if you compile nsp and fftw library is not found, fftpack is used instead). Note that fftw3 uses fast algorithms (in mlog(m)) for every m even prime (but transforms are faster when m is a power of 2).

The optional named argument dim have the following meaning:


  1. One interpretation of the dft is “approximation of some Fourier coefficients” of a T-periodic (complex valued) function g:
                 ∫ T                         1   i2πkt
ˆgk = (g|ϕk ) =    g(t)ϕ¯k (t)dt with ϕk(t) = √-e -T--
              0                           T
    by the rectangular rule using m points uniformly spaced in [0,T], sj = Δtj = (T∕m)j, j = 0,,m 1:
              T m∑ −1      1    i2πksj-  √T--m∑ −1       2πikj
ˆgk ≃ ck := --    g(sj)√--e−  T   = ----    g(sj)e−  m
          m  j=0        T           m   j=0
    This formula could give only m different complex numbers because ck+m = ck for every k . So it is possible to approximate Fourier coefficients only from k = (m 1)2 to k = (m 1)2 (when m is odd). Moreover negative modes are transposed as “positive ones” using the periodicity ck+m = ck, and ĝk is approximated by cmk, for k = 1,2,. This could be written:
        √ --
    --T-                                 ⊤
c =  m  FmG,  with G = [g(s0),...,g(sm−1)]
    This explains that the dft computes (approximated) Fourier coefficients not the usual order. The function fftshift can be used to recover it.
  2. The dft corresponds also to (exact) trigonometric polynomial interpolation of a periodic function at regularly spaced points. If g is such a T-periodic function, and one tries to interpolate it at the sj previously defined by the trigonometric polynomial:
p(t) :=     ckϕk(t),  n = (m − 1)∕2
    then the same formula to compute c is obtained but always with the special order for the coefficients (first corresponds to c0, second to c1, third is c2, etc, last but one is c2 and last is c1).


simple example

  // build a pure periodic signal from a trigonometric polynomial
  // and recover its coefficients using fft
  function x = poly_trigo(t, T, coef, num)
      x = zeros(size(t));
      for i=1:numel(coef), x = x + coef(i)*exp((2*%pi*%i*num(i)/T)*t)/sqrt(T); end
  // build a signal (real valued)
  T = 2;   // period
  coef = [0.5,   -1,  1, 0.86, 1,  -1, 0.5];  // coef
  num  = [-12  , -4, -1, 0   , 1,   4,  12];  // mode number of the coef
  // use a time discretisation and plot the signal
  m = 101;
  t = linspace(0, T, m+1)';
  x = poly_trigo(t, T, coef, num);
  plot2d(t, real(x), style=2)
  xtitle("my signal")
  // recover the coefficients using fft
  c = (sqrt(T)/m)*real(fft(x(1:m)));
  // reorder
  c = fftshift(c);
  // limit between mode -20 to 20 (modes other than -12, -4, -1, 0, 1, 4, 12
  // should be null or numerically negligeable)
  plot2d3((-20:20)',c(31:71),style=2, rect=[-21,-1.2,21,1.2])
  plot2d3(num+0.2, coef, frameflag=0,style=5)
  xtitle("blue: coefs from fft, red: exact coefs")
  // in fact the coefs got from fft should be exact and this must be so
  // if m >= 2*12+1 = 25. Try with m = 25
  m = 25;
  t = linspace(0, T, m+1)';
  x = poly_trigo(t, T, coef, num);
  c = (sqrt(T)/m)*real(fft(x(1:m)));
  // reorder
  c = fftshift(c);
  plot2d3((-12:12)',c,style=2, rect=[-13,-1.2,13,1.2])
  plot2d3(num+0.2, coef, frameflag=0,style=5)
  xtitle("blue: coefs from fft, red: exact coefs")

See also fftshift , ifftshift

Authors fftw: M. Frigo,S.G. Johnson; nsp interface: Bruno Pincon; fftpack: P.N. Swarztrauber