The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Signal and image processing

One-dimensional convolution

The conv() function performs a one-dimensional convolution of two vectors. The orientation of the result vector is the same as the orientation of the first input vector.

Usage:

vec conv(x, y, shape);

The shape argument is optional and can be one of the following:

"full": return the full convolution (default setting), with the size equal to x.n_elem + y.n_elem - 1.
"same": return the central part of the convolution, with the same size as vector x.

The convolution operation is also equivalent to finite impulse response (FIR) filtering.

Examples

[[cpp11::register]] list conv1_(const doubles& x, const doubles& y) {
  vec a = as_col(x);
  vec b = as_col(y);

  vec c = conv(a, b);
  vec d = conv(a, b, "same");

  writable::list out(2);
  out[0] = as_doubles(c);
  out[1] = as_doubles(d);

  return out;
}

Two-dimensional convolution

The conv2() function performs a two-dimensional convolution of two matrices. The orientation of the result matrix is the same as the orientation of the first input matrix.

Usage:

mat conv2(A, B, shape);

The shape argument is optional and can be one of the following:

"full": return the full convolution (default setting), with the size equal to size(A) + size(B) - 1.
"same": return the central part of the convolution, with the same size as matrix A.

Caveats

The implementation of 2D convolution in this version is preliminary.

Examples

[[cpp11::register]] list conv2_(const doubles_matrix<>& x,
  const doubles_matrix<>& y) {
  mat a = as_mat(x);
  mat b = as_mat(y);

  mat c = conv2(a, b);
  mat d = conv2(a, b, "same");

  writable::list out(2);
  out[0] = as_doubles_matrix(c);
  out[1] = as_doubles_matrix(d);

  return out;
}

One-dimensional Fast Fourier Transform

The fft() function computes the fast Fourier transform (FFT) of a vector or matrix. The function returns a complex matrix.

Similarly, ifft() computes the inverse fast Fourier transform (IFFT) of a complex matrix.

The transform is done on each column vector of the input matrix.

Usage:

// real or complex
cx_vec Y = fft(X);
cx_vec Y = fft(X, n);

// complex only
cx_mat Z = ifft(cx_mat Y);
cx_mat Z = ifft(cx_mat Y, n);

The optional n argument specifies the transform length:

If n is larger than the length of the input vector, a zero-padded version of the vector is used.
If n is smaller than the length of the input vector, only the first n elements of the vector are used.

Caveats

The transform is fastest when the transform length is a power of 2 (\(2^k,\: k = 1, 2, 3, \ldots\)).
By default, the function uses an internal FFT algorithm based on KISS FFT. With vendoring, it is possible to use the FFTW3 library for faster execution by modifying the cpp11armadillo header to:

...
#include <Rmath.h>

#define ARMA_USE_FFTW3 // add this line
#include <armadillo.hpp>
...

Examples

[[cpp11::register]] list fft1_(const doubles& x) {
  vec a = as_Col(x);

  cx_vec b = fft(a);
  cx_vec c = ifft(b);

  writable::list out(2);
  writable::list out2(2);
  writable::list out3(2);

  out2[0] = as_doubles(real(b));
  out2[1] = as_doubles(imag(b));

  out3[0] = as_doubles(real(c));
  out3[1] = as_doubles(imag(c));

  out[0] = out2;
  out[1] = out3;

  return out;
}

Two-dimensional Fast Fourier Transform

The fft2() function computes the two-dimensional fast Fourier transform (FFT) of a matrix. The function returns a complex matrix.

Similarly, ifft2() computes the inverse fast Fourier transform (IFFT) of a complex matrix.

Usage:

// real or complex
cx_mat Y = fft2(mat X);
cx_mat Y = fft2(mat X, int n_rows, int n_cols);

// complex only
cx_mat Z = ifft2(cx_mat Y);
cx_mat Z = ifft2(cx_mat Y, int n_rows, int n_cols);

The optional n_rows and n_cols arguments specify the transform size:

If n_rows and n_cols are larger than the size of the input matrix, a zero-padded version of the matrix is used.
If n_rows and n_cols are smaller than the size of the input matrix, only the first n_rows and n_cols elements of the matrix are used.

Caveats

The implementation of the 2D transform in this version is preliminary.
The transform is fastest when both n_rows and n_cols are a power of 2 (\(2^k,\: k = 1, 2, 3, \ldots\)).
By default, the function uses an internal FFT algorithm based on KISS FFT. With vendoring, it is possible to use the FFTW3 library for faster execution by modifying the cpp11armadillo header to:

...
#include <Rmath.h>

#define ARMA_USE_FFTW3 // add this line
#include <armadillo.hpp>
...

Examples

[[cpp11::register]] list fft2_(const doubles_matrix<>& x) {
  mat a = as_mat(x);

  cx_mat b = fft2(a);
  cx_mat c = ifft2(b);

  writable::list out(2);
  writable::list out2(2);
  writable::list out3(2);

  out2[0] = as_doubles(real(b));
  out2[1] = as_doubles(imag(b));

  out3[0] = as_doubles(real(c));
  out3[1] = as_doubles(imag(c));

  out[0] = out2;
  out[1] = out3;

  return out;
}

One-dimensional interpolation

The interp1() function performs one-dimensional interpolation of a function specified by vectors X and Y. The function generates a vector YI that contains interpolated values at locations XI.

Usage:

vec interp1(X, Y, XI, YI);
vec interp1(X, Y, XI, YI, method);
vec interp1(X, Y, XI, YI, method, extrapolation_value);

The method argument is optional and can be one of the following:

"nearest": interpolate using single nearest neighbour.
"linear": linear interpolation between two nearest neighbours (default setting).
"*nearest": as per "nearest", but faster by assuming that X and XI are monotonically increasing.
"*linear": as per "linear", but faster by assuming that X and XI are monotonically increasing.

If a location in XI is outside the domain of X, the corresponding value in YI is set to extrapolation_value.

The extrapolation_value argument is optional; by default, it is datum::nan (not-a-number).

Examples

[[cpp11::register]] doubles interp1_(const int& n) {
  vec x = linspace<vec>(0, 3, n);
  vec y = square(x);

  vec xx = linspace<vec>(0, 3, 2 * n);
  vec yy;

  interp1(x, y, xx, yy);             // use linear interpolation by default
  interp1(x, y, xx, yy, "*linear");  // faster than "linear"
  interp1(x, y, xx, yy, "nearest");

  return as_doubles(yy);
}

Two-dimensional interpolation

The interp2() function performs two-dimensional interpolation of a function specified by matrix Z with coordinates given by vectors X and Y. The function generates a matrix ZI that contains interpolated values at the coordinates given by vectors XI and YI.

Usage:

mat interp2(X, Y, Z, XI, YI, ZI);
mat interp2(X, Y, Z, XI, YI, ZI, method);
mat interp2(X, Y, Z, XI, YI, ZI, method, extrapolation_value);

The method argument is optional and can be one of the following:

"nearest": interpolate using nearest neighbours.
"linear": linear interpolation between nearest neighbours (default setting).

If a coordinate in the 2D grid specified by (XI, YI) is outside the domain of the 2D grid specified by (X, Y), the corresponding value in ZI is set to extrapolation_value.

The extrapolation_value argument is optional; by default, it is datum::nan (not-a-number).

Examples

[[cpp11::register]] doubles_matrix<> interp2_(const int& n) {
  mat Z(n, n, fill::randu);

  vec X = regspace(1, Z.n_cols);  // X = horizontal spacing
  vec Y = regspace(1, Z.n_rows);  // Y = vertical spacing

  vec XI = regspace(X.min(), 1.0/2.0, X.max()); // magnify by approx 2
  vec YI = regspace(Y.min(), 1.0/3.0, Y.max()); // magnify by approx 3

  mat ZI;

  interp2(X, Y, Z, XI, YI, ZI); // use linear interpolation by default

  return as_doubles_matrix(ZI);
}

Find polynomial coefficients for data fitting

The polyfit() function finds the polynomial coefficients for data fitting. The function models a 1D function specified by vectors X and Y as a polynomial of order N and stores the polynomial coefficients in a column vector P.

The given function is modelled as:

\[ y = p_0 x^N + p_1 x^{N-1} + p_2 x^{N-2} + \ldots + p_{N-1} x^1 + p_N \]

where \(p_i\) is the \(i\)-th polynomial coefficient. The coefficients are selected to minimise the overall error of the fit (least squares).

The column vector P has \(N+1\) coefficients.

N must be smaller than the number of elements in X.

Usage:

P = polyfit(X, Y, N);
polyfit(P, X, Y, N);

If the polynomial coefficients cannot be found:

P = polyfit(X, Y, N) resets P and returns an error.
polyfit(P, X, Y, N) resets P and returns a bool set to false without an error.

Examples

[[cpp11::register]] doubles polyfit1_(const int& n, const int& m) {
  vec x = linspace<vec>(0, 1, n);
  vec y = 2*pow(x,2) + 2*x + ones<vec>(n);

  vec p = polyfit(x, y, m);

  return as_doubles(p);
}

Evaluate polynomial

The polyval() function evaluates a polynomial. Given a vector P of polynomial coefficients and a vector X containing the independent values of a 1D function, the function generates a vector Y that contains the corresponding dependent values.

For each x value in vector X, the corresponding y value in vector Y is generated using:

\[ y = p_0 x^N + p_1 x^{N-1} + p_2 x^{N-2} + \ldots + p_{N-1} x^1 + p_N \]

where \(p_i\) is the \(i\)-th polynomial coefficient in vector P.

P must contain polynomial coefficients in descending powers (e.g., generated by the polyfit() function).

Usage:

Y = polyval(P, X);

Examples

[[cpp11::register]] doubles polyval1_(const int& n, const int& m) {
  vec x = linspace<vec>(0, 1, n);
  vec y = 2*pow(x,2) + 2*x + ones<vec>(n);

  vec p = polyfit(x, y, m);
  vec q = polyval(p, x);

  return as_doubles(q);
}

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.