Module: NumRu::FFTW3

Defined in:
lib/numru/fftw3.rb,
lib/numru/fftw3/version.rb,
ext/numru/fftw3/na_fftw3.c

Overview

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!

Constant Summary collapse

VERSION = 
"1.0.2"
SUPPORT_BIGMEM = 
Qfalse
FORWARD = 
INT2NUM(FFTW_FORWARD)
BACKWARD =

Specifier of backward complex FFT. Integer equal to 1.

INT2NUM(FFTW_BACKWARD)
R2HC =

Specifier of real FFT kind. real to “halfcomplex”

INT2NUM(FFTW_R2HC)
HC2R =

Specifier of real FFT kind. “halfcomplex” to real

INT2NUM(FFTW_HC2R)
DHT =

Specifier of real FFT kind. Discrete Hartley Transform

INT2NUM(FFTW_DHT)
REDFT00 =

Specifier of real FFT kind. cosine (even) transform;

logical data length 2*(n-1); inverse is itself

INT2NUM(FFTW_REDFT00)
REDFT01 =

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

INT2NUM(FFTW_REDFT01)
REDFT10 =

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

INT2NUM(FFTW_REDFT10)
REDFT11 =

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

INT2NUM(FFTW_REDFT11)
RODFT00 =

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

INT2NUM(FFTW_RODFT00)
RODFT01 =

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

INT2NUM(FFTW_RODFT01)
RODFT10 =

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

INT2NUM(FFTW_RODFT10)
RODFT11 =

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

INT2NUM(FFTW_RODFT11)

Class Method Summary collapse

Class Method Details

.fft(narray, dir[, dim, dim, ...]) ⇒ Object

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.



293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
 
# File 'ext/numru/fftw3/na_fftw3.c', line 293

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_bk(na, *dims) ⇒ Object

Backward complex FFT

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



113
114
115
 
# File 'lib/numru/fftw3.rb', line 113

def fft_bk(na, *dims)
  fft(na, BACKWARD, *dims)
end

.fft_fw(na, *dims) ⇒ Object

Forward complex FFT with normalization

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



97
98
99
100
101
102
103
104
105
106
107
 
# File 'lib/numru/fftw3.rb', line 97

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

.fft_r2r(narray, kind[, dim, dim, ...]) ⇒ Object

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 FFTW3.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.



625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
 
# File 'ext/numru/fftw3/na_fftw3.c', line 625

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) );
  }

}