This vignette is adapted from the official Armadillo documentation.
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.
[[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;
}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.The implementation of 2D convolution in this version is preliminary.
[[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;
}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:
n is larger than the length of the input vector, a zero-padded version of the vector is used.n is smaller than the length of the input vector, only the first n elements of the vector are used....
#include <Rmath.h>
#define ARMA_USE_FFTW3 // add this line
#include <armadillo.hpp>
...[[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;
}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:
n_rows and n_cols are larger than the size of the input matrix, a zero-padded version of the matrix is used.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.n_rows and n_cols are a power of 2 (\(2^k,\: k = 1, 2, 3, \ldots\))....
#include <Rmath.h>
#define ARMA_USE_FFTW3 // add this line
#include <armadillo.hpp>
...[[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;
}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).
[[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);
}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).
[[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);
}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.[[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);
}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);[[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);
}