module NumRu::FFTW3

Ruby wrapper of FFTW3, a fast discrete Fourier transform library. www.fftw.org

Features

How to use

Copy and paste the following code line-by-line using irb. Or you can run it by saving it in a file fftw3test.rb (say) and executing “ruby fftw3test.rb”.

require "numru/fftw3"
include NumRu

na = NArray.float(8,6)   # float -> will be coerced to complex
na[1,1]=1

# <example 1: complex FFT on all dimensions >

fc = FFTW3.fft(na, FFTW3::FORWARD)/na.length  # forward 2D FFT and normalization
nc = FFTW3.fft(fc, FFTW3::BACKWARD)       # backward 2D FFT (complex) --> 
nb = nc.real              # should be equal to na except round errors  
p (nb - na).abs.max       # => 8.970743058303247e-17 (sufficiently small)

# <example 2: complex FFT on all dimensions >
# Same as example 1 but using more user-friendly wrapper of FFTW3.fft

fc = FFTW3.fft_fw(na)     # forward 2D FFT and normalization
nc = FFTW3.fft_bk(fc)     # backward 2D FFT (complex) --> 
nb = nc.real              # should be equal to na except round errors  
p (nb - na).abs.max       # => 8.970743058303247e-17 (sufficiently small)

# <example 3: complex FFT along a dimension >
fc = FFTW3.fft_fw(na, 0)  # forward 1D FFT along the first dim
nc = FFTW3.fft_bk(fc, 0)  # backward 1D FFT along the first dim
p (nc.real - na).abs.max  # => 1.1102230246251565e-16 (sufficiently small)

# <example 4: complex FFT along a dimension >
fc = FFTW3.fft_fw(na, 1)  # forward 1D FFT along the second dim

# <example 5: real FFT along a dimension >

fc = FFTW3.fft_r2r(na, FFTW3::RODFT00, 0)  # not normalized sine transform along the 1st dim
len = 2*(na.shape[0]+1)    # this is the supposed length of this transformation
nc = FFTW3.fft_r2r(fc/len, FFTW3::RODFT00, 0)  # forward==backward transformation
p (nc-na).abs.max          # => 2.220446049250313e-16 (sufficiently small)

# <example 5b: real FFT on all dimensions >

fc = FFTW3.fft_r2r(na, FFTW3::REDFT11)   # unnormalized cosine transform
len = 4*na.length          # from (2*na.shape[0]) * (2*na.shape[1])
nc = FFTW3.fft_r2r(fc/len, FFTW3::REDFT11)   # forward==backward transformation
p (nc-na).abs.max          # => 6.228190483314256e-17 (sufficiently small)

See the FFTW3 manual for the kinds of supported real FFTs. See Ch.2 of www.fftw.org/fftw3_doc/ (www.fftw.org/fftw3.pdf). Virtually all kinds are supported!

Constants

BACKWARD

Specifier of backward complex FFT. Integer equal to 1.

DHT

Specifier of real FFT kind. Discrete Hartley Transform

FORWARD

Specifier of forward complex FFT. Integer equal to -1.

HC2R

Specifier of real FFT kind. “halfcomplex” to real

R2HC

Specifier of real FFT kind. real to “halfcomplex”

REDFT00

Specifier of real FFT kind. cosine (even) transform; logical data length 2*(n-1); inverse is itself

REDFT01

Specifier of real FFT kind. cosine (even) transform; logical data length 2*n; inverse is REDFT10

REDFT10

Specifier of real FFT kind. cosine (even) transform; logical data length 2*n; inverse is REDFT01

REDFT11

Specifier of real FFT kind. cosine (even) transform; logical data length 2*n; inverse is itself

RODFT00

Specifier of real FFT kind. sine (odd) transform; logical data length 2*(n+1); inverse is itself

RODFT01

Specifier of real FFT kind. sine (odd) transform; logical data length 2*n; inverse is RODFT10

RODFT10

Specifier of real FFT kind. sine (odd) transform; logical data length 2*n; inverse is RODFT01

RODFT11

Specifier of real FFT kind. sine (odd) transform; logical data length 2*n; inverse is itself

VERSION

Public Class Methods

fft(narray, dir [, dim, dim, ...]) click to toggle source

Conducts complex FFT (unnormalized). You can use more user-friendly wrappers #fft_fw and fft_bk.

ARGUMENTS

  • narray (NArray or NArray-compatible Array) : array to be transformed. If real, coerced to complex before transformation. If narray is single-precision and the single-precision version of FFTW3 is installed (before installing this module), this method does a single-precision transform. Otherwise, a double-precision transform is used.

  • dir (FORWARD (which is simply equal to -1; referable as NumRu::FFTW3::FORWARD) or BACKWARD (which is simply 1) ) : the direction of FFT, forward if FORWARD and backward if BACKWARD.

  • optional 3rd, 4th,… arguments (Integer) : Specifies dimensions to apply FFT. For example, if 0, the first dimension is transformed (1D FFT); If -1, the last dimension is used (1D FFT); If 0,2,4, the first, third, and fifth dimensions are transformed (3D FFT); If entirely omitted, ALL DIMENSIONS ARE SUBJECT TO FFT, so 3D FFT is done with a 3D array.

RETURN VALUE

  • a complex NArray

NOTE

  • As in FFTW, return value is NOT normalized. Thus, a consecutive forward and backward transform would multiply the size of data used for transform. You can normalize, for example, the forward transform FFTW.fft(narray, -1, 0, 1) (FFT regarding the first (dim 0) & second (dim 1) dimensions) by dividing with (narray.shape*narray.shape). Likewise, the result of FFTW.fft(narray, -1) (FFT for all dimensions) can be normalized by narray.length.

static VALUE
na_fftw3(int argc, VALUE *argv, VALUE self)
{
  VALUE val;
  volatile VALUE v1;
  struct NARRAY *a1;

  if (argc<2){
    rb_raise(rb_eArgError, "Usage: fft(narray, direction [,dim0,dim1,...])");
  }
  val = argv[0];
  v1 = na_to_narray(val);
  GetNArray(v1,a1);
  if(a1->type <= NA_SFLOAT || a1->type == NA_SCOMPLEX ){
      return( na_fftw3_float(argc, argv, self) );
  } else {
      return( na_fftw3_double(argc, argv, self) );
  }

}
fft_r2r(narray, kind [, dim, dim, ...]) click to toggle source

Conducts real FFT (unnormalized).

ARGUMENTS

RETURN VALUE

  • a real NArray

NOTE

  • As in FFTW, return value is NOT normalized, and the length needed normalize is different for each kind.

static VALUE
na_fftw3_r2r(int argc, VALUE *argv, VALUE self)
{
  VALUE val;
  volatile VALUE v1;
  struct NARRAY *a1;

  if (argc<2){
    rb_raise(rb_eArgError, "Usage: fft_r2r(narray, kinds [,dim0,dim1,...])");
  }
  val = argv[0];
  v1 = na_to_narray(val);
  GetNArray(v1,a1);
  if(a1->type <= NA_SFLOAT || a1->type == NA_SCOMPLEX ){
      return( na_fftw3_r2r_float(argc, argv, self) );
  } else {
      return( na_fftw3_r2r_double(argc, argv, self) );
  }

}

Public Instance Methods

fft_bk(na, *dims) click to toggle source

Backward complex FFT

This method simply calls FFW3.fft(na, FFW3::BACKWARD, *dims)

# File lib/numru/fftw3.rb, line 98
def fft_bk(na, *dims)
  fft(na, BACKWARD, *dims)
end
fft_fw(na, *dims) click to toggle source

Forward complex FFT with normalization

This calls FFW3.fft(na, FFW3::FORWARD, *dims) and normalizes the result by dividing by the length

# File lib/numru/fftw3.rb, line 82
def fft_fw(na, *dims)
  fc = fft(na, FORWARD, *dims)
  if dims.length == 0
    len = na.total
  else
    len = 1
    shape = na.shape
    dims.each{|d| len *= shape[d]}
  end
  fc / len
end