Module: Convolver

Defined in:
lib/convolver.rb,
lib/convolver/version.rb,
ext/convolver/convolver.c

Overview

Valid cross-correlation operations for Numo::NArray values.

Constant Summary collapse

MAX_RANK =

Maximum number of dimensions supported by the direct native implementation.

16
VERSION =

Current gem version.

'1.0.1'

Class Method Summary collapse

Class Method Details

.convolve(signal, kernel) ⇒ Numo::SFloat

Chooses the likely fastest implementation for a valid cross-correlation.

The inputs must have the same rank, and the kernel must not be larger than the signal in any dimension. The result shape is:

signal.shape.zip(kernel.shape).map { |signal_size, kernel_size| signal_size - kernel_size + 1 }

Parameters:

  • signal (Numo::NArray) —

    input values

  • kernel (Numo::NArray) —

    correlation kernel

Returns:

  • (Numo::SFloat) —

    valid cross-correlation result

Raises:

  • (ArgumentError) —

    if the inputs have incompatible ranks or shapes



25
26
27
28
29
30
31
32
33
34
35
36
 
# File 'lib/convolver.rb', line 25

def convolve(signal, kernel)
  validate_inputs!(signal, kernel)
  return convolve_basic(signal, kernel) if signal.size < 1000 || kernel.size < 100

  basic_time_predicted = predict_convolve_basic_time(signal, kernel)
  return convolve_basic(signal, kernel) if basic_time_predicted < 0.1

  fft_time_predicted = predict_convolve_fft_time(signal, kernel)
  return convolve_fft(signal, kernel) if fft_time_predicted < 2 * basic_time_predicted

  convolve_basic(signal, kernel)
end

.convolve_basic(signal, kernel) ⇒ Object 



25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
 
# File 'ext/convolver/convolver.c', line 25

static VALUE convolver_convolve_basic(VALUE self, VALUE signal, VALUE kernel) {
  volatile VALUE signal_value;
  volatile VALUE kernel_value;
  volatile VALUE result_value;
  narray_t *signal_narray;
  narray_t *kernel_narray;
  int rank;
  int i;
  size_t signal_shape[LARGEST_RANK];
  size_t kernel_shape[LARGEST_RANK];
  size_t result_shape[LARGEST_RANK];
  size_t numo_result_shape[LARGEST_RANK];

  (void)self;

  if (!rb_obj_is_kind_of(signal, numo_cNArray) || !rb_obj_is_kind_of(kernel, numo_cNArray)) {
    rb_raise(rb_eArgError, "signal and kernel must be Numo::NArray values");
  }

  signal_value = rb_funcall(numo_cSFloat, rb_intern("cast"), 1, signal);
  kernel_value = rb_funcall(numo_cSFloat, rb_intern("cast"), 1, kernel);
  if (!RTEST(na_check_contiguous(signal_value))) {
    signal_value = na_copy(signal_value);
  }
  if (!RTEST(na_check_contiguous(kernel_value))) {
    kernel_value = na_copy(kernel_value);
  }
  GetNArray(signal_value, signal_narray);
  GetNArray(kernel_value, kernel_narray);

  if (signal_narray->size == 0 || kernel_narray->size == 0) {
    rb_raise(rb_eArgError, "signal and kernel must not be empty");
  }

  if (signal_narray->ndim != kernel_narray->ndim) {
    rb_raise(rb_eArgError, "signal and kernel must have equal rank");
  }
  if (signal_narray->ndim > LARGEST_RANK) {
    rb_raise(rb_eArgError, "maximum supported rank is %d", LARGEST_RANK);
  }

  rank = signal_narray->ndim;
  copy_shape(rank, signal_narray->shape, signal_shape);
  copy_shape(rank, kernel_narray->shape, kernel_shape);

  for (i = 0; i < rank; i++) {
    if (signal_shape[i] < kernel_shape[i]) {
      rb_raise(rb_eArgError, "kernel must not be larger than signal in any dimension");
    }
    result_shape[i] = signal_shape[i] - kernel_shape[i] + 1;
    numo_result_shape[rank - i - 1] = (size_t)result_shape[i];
  }

  result_value = nary_new(numo_cSFloat, rank, numo_result_shape);
  convolve_raw(
    rank, signal_shape, (float *)na_get_pointer_for_read(signal_value),
    rank, kernel_shape, (float *)na_get_pointer_for_read(kernel_value),
    rank, result_shape, (float *)na_get_pointer_for_write(result_value)
  );

  return result_value;
}

.convolve_fft(signal, kernel) ⇒ Numo::SFloat

Uses PocketFFT to calculate a valid cross-correlation.

Parameters:

  • signal (Numo::NArray) —

    input values

  • kernel (Numo::NArray) —

    correlation kernel

Returns:

  • (Numo::SFloat) —

    valid cross-correlation result

Raises:

  • (ArgumentError) —

    if the inputs have incompatible ranks or shapes



44
45
46
47
48
49
50
 
# File 'lib/convolver.rb', line 44

def convolve_fft(signal, kernel)
  validate_inputs!(signal, kernel)
  ranges = kernel.shape.zip(signal.shape).map { |kernel_size, signal_size| (kernel_size - 1)...signal_size }
  full_convolution = Numo::Pocketfft.fftconvolve(signal, kernel.reverse)

  Numo::SFloat.cast(full_convolution[*ranges])
end

.convolve_fftw3(signal, kernel) ⇒ Numo::SFloat

Deprecated.

Use convolve_fft; Convolver no longer uses FFTW3.

Compatibility alias for the former FFTW3-backed implementation.

Returns:

  • (Numo::SFloat) —

    valid cross-correlation result



56
57
58
59
 
# File 'lib/convolver.rb', line 56

def convolve_fftw3(signal, kernel)
  warn 'Convolver.convolve_fftw3 is deprecated; use .convolve_fft instead', uplevel: 1
  convolve_fft(signal, kernel)
end

.predict_convolve_basic_time(signal, kernel) ⇒ Float

Estimates the relative cost of convolve_basic.

Parameters:

  • signal (Numo::NArray) —

    input values

  • kernel (Numo::NArray) —

    correlation kernel

Returns:

  • (Float) —

    machine-specific relative cost estimate



77
78
79
80
81
 
# File 'lib/convolver.rb', line 77

def predict_convolve_basic_time(signal, kernel)
  validate_inputs!(signal, kernel)
  outputs = result_shape(signal.shape, kernel.shape).inject(:*)
  4.54e-12 * (outputs * signal.size * kernel.size)
end

.predict_convolve_fft_time(signal, kernel) ⇒ Float

Estimates the relative cost of convolve_fft.

Parameters:

  • signal (Numo::NArray) —

    input values

  • kernel (Numo::NArray) —

    correlation kernel

Returns:

  • (Float) —

    machine-specific relative cost estimate



66
67
68
69
70
 
# File 'lib/convolver.rb', line 66

def predict_convolve_fft_time(signal, kernel)
  validate_inputs!(signal, kernel)
  output_size = result_shape(signal.shape, kernel.shape).inject(:*)
  16 * 4.55e-08 * output_size * Math.log(output_size)
end