module NumRu::FFTW3

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

Features¶ ↑

Uses NArray (github.com/masa16/narray). (Also it supports NumRu::NArray as well, if this library is compiled to use it).
Multi-dimensional complex and real FFT.
Supports both double and single float transforms.

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

narray (NArray or NArray-compatible Array) : array to be transformed. Cannot be complex.
kind : specifies the kind of FFT. One of the following: R2HC, HC2R, DHT, REDFT00, REDFT01, REDFT10, REDFT11, RODFT00, RODFT01, RODFT10, or RODFT11 (referable as NumRu::FFTW3::REDFT10). See Ch.2 of www.fftw.org/fftw3_doc/ (www.fftw.org/fftw3.pdf).
optional 3rd, 4th,… arguments (Integer) : See ::fft.

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