--- title: "Signal and image processing" output: rmarkdown::html_vignette bibliography: "references.bib" vignette: > %\VignetteIndexEntry{Signal and image processing} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` This vignette is adapted from the official Armadillo [documentation](https://arma.sourceforge.net/docs.html). # One-dimensional convolution {#conv} 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: ```cpp 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 ```cpp [[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 {#conv2} 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: ```cpp 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 ```cpp [[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 {#fft-ifft} 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: ```cpp // 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: ```cpp ... #include #define ARMA_USE_FFTW3 // add this line #include ... ``` ## Examples ```cpp [[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 {#fft-ifft-2} 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: ```cpp // 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: ```cpp ... #include #define ARMA_USE_FFTW3 // add this line #include ... ``` ## Examples ```cpp [[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 {#interp1} 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: ```cpp 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 ```cpp [[cpp11::register]] doubles interp1_(const int& n) { vec x = linspace(0, 3, n); vec y = square(x); vec xx = linspace(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 {#interp2} 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: ```cpp 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 ```cpp [[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 {#polyfit} 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: ```cpp 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 ```cpp [[cpp11::register]] doubles polyfit1_(const int& n, const int& m) { vec x = linspace(0, 1, n); vec y = 2*pow(x,2) + 2*x + ones(n); vec p = polyfit(x, y, m); return as_doubles(p); } ``` # Evaluate polynomial {#polyval} 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: ```cpp Y = polyval(P, X); ``` ## Examples ```cpp [[cpp11::register]] doubles polyval1_(const int& n, const int& m) { vec x = linspace(0, 1, n); vec y = 2*pow(x,2) + 2*x + ones(n); vec p = polyfit(x, y, m); vec q = polyval(p, x); return as_doubles(q); } ```