This code provides users with commonly used window functions

window_td.cpp

/* window_td.cpp v1.3
 *
 * Copyright (C) 2018, 2021, 2023 BD1ES <bd1es@hotmail.com>
 */

/*
 * Dependencies:
 *      PocketFFT:
 *          https://gitlab.mpcdf.mpg.de/mtr/pocketfft
 *          A customized version of "pocketfft.c" is needed for the ESP32 toolchains.
 *          Download from: https://gist.github.com/bd1es/777a4e2b825f8e0ce44e3bfae0ae643a
 
 *      PGFFT:
 *          https://www.shoup.net/PGFFT
 */

#include <cmath>
#include <bits/stdc++.h>

#include <pocketfft.h>
#include <PGFFT.h>

#include <window_td.h>

// ------------------------------------------------------------------------------------------------

#ifndef M_PI
static const double M_PI = double(3.1415926535897932384626);
#endif

std::vector <double> window_td::gencossum(size_t n, const std::vector <double> &a, form_opt opt)
{
    if(a.size() < 2){
        return boxcar(0);
    }else if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = opt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = a[0];
        for(size_t k = 1; k < a.size(); k++){
            y[i] += a[k] * cos(M_PI*k*(1.0-2.0*i/N));
        }
    }

    return y;
}

double window_td::sinc(double x)
{
    x *= M_PI;
    if(x != 0){
        return sin(x)/x;
    }else{
        return 1;
    }
}

std::vector <double> window_td::boxcar(size_t n)
{
    return std::vector <double> (n, 1);
}

std::vector <double> window_td::rectwin(size_t n)
{
    return boxcar(n);
}

std::vector <double> window_td::bartlett(size_t n, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    size_t N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = i <= N/2 ? 2.0*i/N : 2.0 - 2.0*i/N;
    }

    return y;
}

std::vector <double> window_td::triang(size_t n, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    size_t N = fopt == symmetric ? n : n+1;
    for(size_t i = 0; i < n; ++i){
        y[i] = 1.0 - abs((1.0-N+2.0*i) / (N+N%2));
    }

    return y;
}

std::vector <double> window_td::parzenwin(size_t n, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        double x = abs(i - N/2.0);
        if(x > (N/4.0)){
            y[i] = 2.0*pow(1.0-2.0*x/(1.0+N), 3);
        }else{
            y[i] = 1.0 - 6.0*pow(2.0*x/(1.0+N), 2) * (1.0-2.0*x/(1.0+N));
        }
    }

    return y;
}

std::vector <double> window_td::welchwin(size_t n, form_opt opt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = opt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = 1.0 - pow((2.0*i-N)/N, 2);
    }

    return y;
}

std::vector <double> window_td::bohmanwin(size_t n, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        double x = 2.0 * abs(i-N/2.0) / N;
        y[i] = (1.0-x) * cos(M_PI*x) + sin(M_PI*x)/M_PI;
    }
    y[0] = 0;

    if(fopt == symmetric){
        y[n-1] = 0;
    }

    return y;
}

std::vector <double> window_td::barthannwin(size_t n, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        double x = i/N - 0.5;
        y[i] = 0.62 - 0.48*abs(x) + 0.38*cos(2.0*M_PI*x);
    }

    return y;
}

std::vector <double> window_td::hamming(size_t n, form_opt opt, bool precise, bool optimal)
{
    static const std::vector <double> a[3] = {
        {0.53836, 0.46164},
        {25.0/46, 21.0/46},
        {0.54, 0.46},
    };

    if(optimal){
        return gencossum(n, a[0], opt);
    }else if(precise){
        return gencossum(n, a[1], opt);
    }else{
        return gencossum(n, a[2], opt);
    }
}

std::vector <double> window_td::hann(size_t n, form_opt opt)
{
    static const std::vector <double> a = {0.5, 0.5};

    return gencossum(n, a, opt);
}

std::vector <double> window_td::hanning(size_t n, form_opt opt)
{
    return hann(n, opt);
}

std::vector <double> window_td::blackman(size_t n, form_opt opt, bool exact)
{
    static const std::vector <double> a[2] = {
        {0.42, 0.5, 0.08},
        {7938.0/18608, 9240.0/18608, 1420.0/18608},
    };

    if(!exact){
        return gencossum(n, a[0], opt);
    }else{
        return gencossum(n, a[1], opt);
    }
}

std::vector <double> window_td::nuttallwin(size_t n, form_opt opt, nw_index idx)
{
    // This is from https://pure.mpg.de/rest/items/item_152164_2/component/file_152163/content
    static const std::vector <double> a[7] = {
        {3.0/8, 4.0/8, 1.0/8},                          // Nuttall3
        {0.40897, 0.5, 0.09103},                        // Nuttall3a
        {0.4243801, 0.4973406, 0.0782793},              // Nuttall3b
        {10.0/32, 15.0/32, 6.0/32, 1.0/32},             // Nuttall4
        {0.338946, 0.481973, 0.161054, 0.018027},       // Nuttall4a
        {0.355768, 0.487396, 0.144232, 0.012604},       // Nuttall4b
        {0.3635819, 0.4891775, 0.1365995, 0.0106411},   // Nuttall4c
    };

    if((idx>=NW3) && (idx<=NW4C)){
        return gencossum(n, a[idx], opt);
    }else{
        return boxcar(0);
    }
}

/* This function is equivalent to a single-precision 4-term Hans-Helge cosine-sum window. And it's
 * called "nuttall" in SciPy because it's "a minimum 4-term Blackman-Harris according to Nuttall."
 * In other words, Nuttall created this window in his paper "Some windows with very good side-lobe
 * behaviour, IEEE Trans, Feb 1981."
 */
std::vector <double> window_td::blackmannuttall(size_t n, form_opt opt)
{
    return nuttallwin(n, opt, NW4C);
}

std::vector <double> window_td::blackmanharris(size_t n, form_opt opt)
{
    static const std::vector <double> a = {0.35875, 0.48829, 0.14128, 0.01168};

    return gencossum(n, a, opt);
}

std::vector <double> window_td::flattopwin(size_t n, form_opt opt, bool optimal)
{
    static const std::vector <double> a[2] = {
        {1/4.6402, 1.93/4.6402, 1.29/4.6402, 0.388/4.6402, 0.0322/4.6402},
        // This is from http://www.iowahills.com/Example%20Code/WebFFTCode.txt:
        {1, 1.93293488969227, 1.28349769674027, 0.38130801681619, 0.02929730258511},
    };

    if(!optimal){
        return gencossum(n, a[0], opt);
    }else{
        std::vector <double> y = gencossum(n, a[1], opt);

        double max = y[n/2];
        for(size_t i = 0; i < n; ++i){
            y[i] /= max;
        }

        return y;
    }
}

std::vector <double> window_td::sinwin(size_t n, double alpha, form_opt opt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = opt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = pow(sin(M_PI*double(i)/N), alpha);
    }

    return y;
}

std::vector <double> window_td::lanczoswin(size_t n, double alpha, form_opt opt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = opt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = pow(sinc(2.0*i/N - 1.0), alpha);
    }

    return y;
}

std::vector <double> window_td::tukeywin(size_t n, double r, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    if(r >= 1){
        return hanning(n, fopt);
    }else if(r <= 0){
        return boxcar(n);
    }else{
        size_t N = fopt == symmetric ? n : n+1;
        for(size_t i=0; i < (N/2); ++i){
            y[i] = (1.0 + cos(M_PI*(2.0*i/(N-1)/r-1.0))) / 2.0;
        }
        for(size_t i = floor(r*(N-1)/2.0) + 1; i < ((N+1)/2); ++i){
            y[i] = 1;
        }
        for(size_t i = N/2; i < n; ++i){
            y[i] = y[N-1 -i];
        }
    }

    return y;
}

std::vector< double > window_td::plancktaper(size_t n, double epsilon, form_opt opt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector< double > y(n);

    double tau = opt == symmetric ? n-1 : n;
    double t_cut = tau * (0.5-epsilon);
    for(size_t i = 0; i < n; ++i){
        double t = abs(i - tau/2.0);
        double z = (t_cut-tau/2.0) * (1.0/(t-t_cut) + 1.0/(t-tau/2.0));
        y[i] = t > t_cut ? t < (tau/2.0) ? 1.0/(exp(z)+1.0) : 0 : 1;
    }

    return y;
}

std::vector <double> window_td::gaussian(size_t n, double a, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = exp(-pow(a*(i-N/2.0), 2)/2.0);
    }

    return y;
}

std::vector <double> window_td::gausswin(size_t n, double a, form_opt fopt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt == symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = exp(-pow(a*(2.0*i-N)/N, 2)/2.0);
    }

    return y;
}

std::vector <double> window_td::kaiser(size_t n, double beta, form_opt fopt)
{
    auto i0 = [] (double beta) {
        double sigma = 0, iter = 1;
        beta /= 2.0;
        for(int i = 1; i <= 16; ++i){
            iter *= beta/i;
            sigma += iter*iter;
        }

        return sigma+1;
    };
    
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    double N = fopt==symmetric ? n-1 : n;
    for(size_t i = 0; i < n; ++i){
        y[i] = i0(beta*sqrt(1.0-pow(1.0-2.0*i/N, 2))) / i0(beta);
    }

    return y;
}

/* Implementation of an arbitrary-length complex in-place FFT utilizing PocketFFT or PGFFT.
 * It places DC in bin 0 and adapts the reverse transformation by 1/n for n-FFT. Note that
 * if an internal procedure fails (to request memory,) the function returns false.
 *
 * PocketFFT is distributed under a 3-clause BSD style license. See here:
 *     https://gitlab.mpcdf.mpg.de/mtr/pocketfft/-/blob/master/LICENSE.md
 * Copyright (C) 2008-2018 Max-Planck-Society
 * author Martin Reinecke
 *
 * PGFFT is C++ code for performing a complex FFT, distributed under a FreeBSD license.
 * Copyright (C) 2019, Victor Shoup
 */
enum fft_kernel {k_PocketFFT, k_PGFFT};
static bool _fft(std::vector< std::complex <double>> &v, bool forward, fft_kernel k = k_PGFFT)
{
    auto pocketfft = [forward, &v] (size_t n)
    {
        cfft_plan plan = make_cfft_plan(n);
        if(!plan) return false;

        if(forward){
            cfft_forward(plan, reinterpret_cast <double *> (&v[0]), 1.0);
        }else{
            cfft_backward(plan, reinterpret_cast <double *> (&v[0]), 1.0/n);
        }
        destroy_cfft_plan(plan);

        return true;
    };

    auto pgfft = [forward, &v] (size_t n)
    {
        PGFFT pgfft(n);
        pgfft.apply(&v[0]);
        
        if(!forward){
	    // Reverse the vector.
	    if (n>2) reverse(v.begin()+1, v.end());
        
            for(auto &x : v) x /= n;
        }
        
        return true;
    };

    size_t n = v.size();
    if(n==0) return false;

    if(k == k_PocketFFT){    
        if(!pocketfft(n)) return false;
    }else{
        if(!pgfft(n)) return false;
    }

    return true;
}

/* The Dolph-Chebyshev window W(k) is defined in the frequency domain by the expression:
 *
 *        cheby(n-1, beta * cos(pi * k/n))
 * W(k) = --------------------------------,         0 <= k <= (n-1)
 *                cheby(n-1, beta)
 *
 * with beta = cosh(1/(n-1) * acosh(10^(at/20)),
 *      at = stop band attenuation in dB
 */
std::vector <double> window_td::chebwin(size_t n, double at, form_opt fopt)
{
    // Nth-order Chebyshev polynomial at the point x.
    auto cheby = [] (size_t nth, double x) {
        if(abs(x) > 1){
            return cosh(double(nth) * acosh(std::complex <double> (x)));
        }else{
            return cos(double(nth) * acos(std::complex <double> (x)));
        }
    };

    if(n < 2){
        return boxcar(n);
    }

    std::vector <double> y(n);

    size_t nx = fopt==symmetric ? n : n+1;
    std::vector <std::complex <double>> x(nx);
    double beta = cosh(acosh(pow(10, at/20)) / (nx-1));

    for(size_t k = 0; k < nx; ++k){
        x[k] = cheby(nx-1, beta * cos(M_PI*k/nx));
    }

    // Perform iFFT and normalise the value of y.
    if(nx % 2){
        if(!_fft(x, true)){
            return boxcar(0);
        }

        size_t M = (nx+1)/2;
        size_t j = 0;
        for(size_t i = M-1; i > 0; --i){
            y[j++] = x[i].real() / x[0].real();
        }
        for(size_t i = 0; i < (fopt==symmetric ? M : M-1); ++i){
            y[j++] = x[i].real() / x[0].real();
        }
    }else{
        // Perform e^(j * 0.5 omega) when the length is even.
        for(size_t i = 0; i < nx; ++i){
            x[i] *= exp(std::complex <double> (0, M_PI*i/nx));
        }

        if(!_fft(x, true)){
            return boxcar(0);
        }

        size_t M = nx/2+1;
        size_t j = 0;
        for(size_t i = M-1; i > 1; --i){
            y[j++] = x[i].real() / x[0].real();
        }
        for(size_t i = 0; i < (fopt==symmetric ? M : M-1); ++i){
            y[j++] = x[i].real() / x[0].real();
        }
    }

    return y;
}

// ported from https://github.com/scipy/scipy/blob/master/scipy/signal/windows/windows.py
std::vector <double> window_td::taylorwin(size_t n, size_t nbar, double sll, bool norm,
                                          form_opt opt)
{
    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n, 1);

    double B = pow(10, sll/20);
    double A = acosh(B) / M_PI;
    double s2 = pow(nbar, 2) / (pow(A, 2) + pow((nbar-0.5), 2));

    double N = opt == symmetric ? n : n+1;
    for(size_t mi = 1; mi < nbar; ++mi){
        double numerator = mi % 2 ? 1 : -1;
        double denominator = 2;
        for(size_t i = 1; i < nbar; ++i){
            numerator *= 1 - pow(mi, 2) / (pow(A, 2) + pow(i-0.5, 2)) / s2;
            denominator *= i != mi ? 1 - pow(mi, 2) / pow(i, 2) : 1;
        }

        double Fm = numerator/denominator;
        for(size_t i = 0; i < n; ++i){
            y[i] += 2 * Fm * cos(2*M_PI * mi*((i-N/2.0+0.5)/N));
        }
    }

    if(norm){
        double max = y[n/2];
        for(size_t i = 0; i < n; ++i){
            y[i] /= max;
        }
    }

    return y;
}

std::vector <double> window_td::sprwin(size_t n, double beta,
                                       std::vector <double> (&winfunc)(size_t, double, form_opt),
                                       double alpha, form_opt opt)
{
    auto sum = [] (std::vector <double> &v, size_t n){
        double s = 0;
        for(size_t i = 0; (i<n) && (i<v.size()); ++i) s += v[i];
        return s;
    };

    if(n < 2){
        return boxcar(n);
    }
    std::vector <double> y(n);

    if(opt == periodic) n++;
    size_t N = (n-1)/2 + 1;

    std::vector <double> k = winfunc(N+1, beta, symmetric);

    double denom = sum(k, N+1);

    for(size_t i = 0; i < N; ++i){
        y[i] = pow(sum(k, i+1)/denom, 1/alpha);
    }
    for(size_t i = N; i < y.size(); ++i){
        y[i] = y[2*N-1-i-(n%2)];
    }

    return y;
}

window_td.h

/* window_td.h v1.3
 *
 * This code provides users with commonly used window functions for handling filter designs or
 * developing apps with the ESP-IDF or ESP32 Arduino toolchains. The majority of these windows
 * are designed to be as accurate and easy to use as Octave functions. A few of them also give
 * users the option to change the responses by choosing different window definitions.
 *
 * Copyright (C) 2018, 2021, 2023 BD1ES <bd1es@hotmail.com>
 */

/*
 * This program takes advantage of the following open-source software:
 *      PocketFFT, https://gitlab.mpcdf.mpg.de/mtr/pocketfft
 *      PGFFT, https://www.shoup.net/PGFFT
 *      Octave, https://www.gnu.org/software/octave
 *      Octave Forge, https://octave.sourceforge.io/signal/index.html
 *      SciPy, https://github.com/scipy/scipy/tree/master/scipy/signal/windows
 *
 * References:
 *      Wikipedia, window function:
 *          https://en.wikipedia.org/wiki/Window_function
 
 *      Wikipedia, Kaiser window:
 *          https://en.wikipedia.org/wiki/Kaiser_window
 
 *      G. Heinzel, A. Rüdiger and R. Schilling - Max-Planck-Institut für Gravitationsphysik
 *          Spectrum and spectral density estimation by the Discrete Fourier transform (DFT),
 *          including a comprehensive list of window functions and some new flat-top windows.
 *          https://pure.mpg.de/rest/items/item_152164_2/component/file_152163/content
 */

#ifndef WINDOW_TD_H
#define WINDOW_TD_H

#include <vector>
#include <complex>

namespace window_td
{
    // The symmetry of the window.
    enum form_opt {symmetric, periodic};    // A periodic window of size n equals its symmetrical
                                            // form of size n+1 excluding the last coefficient.

    // -------- Window functions in /m/signal --------

    /* Bartlett
     * /usr/share/octave/4.2.2/m/signal/bartlett.m
     */
    std::vector <double> bartlett(size_t n, form_opt fopt = symmetric);

    /* Blackman
     * /usr/share/octave/4.2.2/m/signal/blackman.m
     *
     * This function generates Blackman coefficients with the 18 dB/oct roll-off. A response of
     * 6 dB/oct for "the exact Blackman" series is also provided as an option.
     */
    std::vector <double> blackman(size_t n, form_opt opt = symmetric, bool exact = false);

    /* Hamming
     * /usr/share/octave/4.2.2/m/signal/hamming.m
     *
     * This function produces Hamming coefficients using the standard expression, with the first
     * term a0 equal to 0.54.
     *
     * There are two options for selecting sidelobe appearances. One example is Hamming's precise
     * series (a0=25/46,) which can be obtained by setting the "precise." Another alternative is
     * the optimal sequence, which appears with an equiripple response (where a0=0.53836) and can
     * be chosen by enabling the "optimal."
     */
    std::vector <double> hamming(size_t n, form_opt opt = symmetric, bool precise = false,
                                 bool optimal = false);

    /* Hann's (raised-cosine) function, also known as the hanning window.
     * /usr/share/octave/4.2.2/m/signal/hanning.m
     */
    std::vector <double> hanning(size_t n, form_opt opt = symmetric);

    // -------- Window functions in the signal package --------

    /* Modified Bartlett-Hann
     * /usr/share/octave/packages/signal-1.3.2/barthannwin.m
     */
    std::vector <double> barthannwin(size_t n, form_opt fopt = symmetric);

    /* Blackman-Harris.
     * /usr/share/octave/packages/signal-1.3.2/blackmanharris.m
     */
    std::vector <double> blackmanharris(size_t n, form_opt opt = symmetric);

    /* Blackman-Nuttall
     * /usr/share/octave/packages/signal-1.3.2/blackmannuttall.m
     */
    std::vector <double> blackmannuttall(size_t n, form_opt opt = symmetric);

    /* Bohman
     * /usr/share/octave/packages/signal-1.3.2/bohmanwin.m
     */
    std::vector <double> bohmanwin(size_t n, form_opt fopt = symmetric);

    /* Boxcar
     * /usr/share/octave/packages/signal-1.3.2/boxcar.m
     */
    std::vector <double> boxcar(size_t n);

    /* Dolph-Chebyshev
     * /usr/share/octave/packages/signal-1.3.2/chebwin.m
     */
    std::vector <double> chebwin(size_t n, double at = 100.0, form_opt fopt = symmetric);

    /* Flat-top
     * /usr/share/octave/packages/signal-1.3.2/flattopwin.m
     *
     * This function produces coefficients with a flat top center using the standard expression.
     * To create a normalized double-precision sequence, please enable the "optimal."
     */
    std::vector <double> flattopwin(size_t n, form_opt opt = symmetric, bool optimal = false);

    /* Gaussian convolution
     * /usr/share/octave/packages/signal-1.3.2/gaussian.m
     */
    std::vector <double> gaussian(size_t n, double a = 1, form_opt fopt = symmetric);

    /* Gaussian
     * /usr/share/octave/packages/signal-1.3.2/gausswin.m
     */
    std::vector <double> gausswin(size_t n, double a = 2.5, form_opt fopt = symmetric);

    /* Hann
     * /usr/share/octave/packages/signal-1.3.2/hann.m
     */
    std::vector <double> hann(size_t n, form_opt opt = symmetric);

    /* Kaiser
     * /usr/share/octave/packages/signal-1.3.2/kaiser.m
     */
    std::vector <double> kaiser(size_t n, double beta = 0.5, form_opt fopt = symmetric);

    /* Nuttall
     * /usr/share/octave/packages/signal-1.3.2/nuttallwin.m
     *
     * This function has been defined to generate coefficients with a window known as "Nuttall4b,"
     * one of seven analogous windows established by Nuttall for "different aims." These windows
     * have indexes defined in "nw_index," allowing users to choose one by changing "idx."
     */
    enum nw_index {NW3=0, NW3A=1, NW3B=2, NW4=3, NW4A=4, NW4B=5, NW4C=6};
    std::vector <double> nuttallwin(size_t n, form_opt opt = symmetric, nw_index idx = NW4B);

    /* Parzen
     * /usr/share/octave/packages/signal-1.3.2/parzenwin.m
     */
    std::vector <double> parzenwin(size_t n, form_opt fopt = symmetric);

    /* Rectangular
     * /usr/share/octave/packages/signal-1.3.2/rectwin.m
     */
    std::vector <double> rectwin(size_t n);

    /* Triangular
     * /usr/share/octave/packages/signal-1.3.2/triang.m
     */
    std::vector <double> triang(size_t n, form_opt fopt = symmetric);

    /* Tukey (cosine-tapered)
     * /usr/share/octave/packages/signal-1.3.2/tukeywin.m
     */
    std::vector <double> tukeywin(size_t n, double r = 0.5, form_opt fopt = symmetric);

    /* Ultraspherical
     * /usr/share/octave/packages/signal-1.3.2/ultrwin.m
     *
     * Notice that:
     * - This function needs that the size n be greater than 7; else, a boxcar will be returned.
     * - If "xmu_" is not NULL, the calculated "XMU" will be placed in its specified location.
     * - The source of this function is in "window_td_ultrwin.cpp" for the license consistency.
     */
    enum ultr_par {beta, att, latt, xmu};
    std::vector <double> ultrwin(size_t n, double mu, double par, ultr_par type, int even_norm = 0,
                                 double *xmu_ = NULL, form_opt fopt = symmetric);

    /* Welch
     * /usr/share/octave/packages/signal-1.3.2/welchwin.m
     */
    std::vector <double> welchwin(size_t n, form_opt opt = symmetric);

    // -------- Some additional window functions --------

    /* Lanczos (sinc)
     */
    std::vector <double> lanczoswin(size_t n, double alpha = 1, form_opt opt = symmetric);

    /* Planck taper
     */
    std::vector <double> plancktaper(size_t n, double epsilon = 0.1, form_opt opt = symmetric);

    /* Power of sine
     */
    std::vector <double> sinwin(size_t n, double alpha = 1, form_opt opt = symmetric);

    /* Taylor
     * This function returns the Taylor window coefficients.
     *
     * Parameters
     * n:    determines the size of the window.
     * nbar: the number of sidelobes that appear almost constant adjacent to the mainlobe.
     *       Note that if the "nbar" is less than three, the function will return a boxcar.
     * sll:  specifies the desired sidelobe suppression in dBs in relation to the DC gain
     *       of the mainlobe.
     * norm: When obtaining canonical Taylor coefficients, "norm" should be turned off. As a
     *       result, the DC gain stays unitary (0 dB) and all sidelobes are "sll" dBs lower.
     * opt:  controls the symmetry of the window.
     *
     * (This function was ported from SciPy.)
     */
    std::vector <double> taylorwin(size_t n, size_t nbar = 4, double sll = 100.0, bool norm = true,
                                   form_opt opt = symmetric);

    /* This is an exploratory function. It will be removed in future releases.
     * This function employs the Kaiser-Bessel-derived (KBD) window and is intended to provide
     * KBD-like window coefficients with tunable center flatness and roll-off characteristics.
     */
    std::vector <double> sprwin(size_t n,                   // the size of the window
                                double beta = 5,            // tune the "winfunc" when possible
                                std::vector <double> (&winfunc)(size_t, double, form_opt) = kaiser,
                                double alpha = 2,           // power of 1/alpha, for "sprwin"
                                form_opt opt = symmetric    // the symmetry of the returned value
                               );

    // -------- Some commonly used functions/templates created for this namespace --------

    /* Creates cosine-sum windows
     */
    std::vector <double> gencossum(size_t n, const std::vector <double> &a, form_opt opt);

    /* An Octave-compatible sinc(x)
     * /usr/share/octave/4.2.2/m/signal/sinc.m
     */
    double sinc(double x);
}

#endif // WINDOW_TD_H

window_td_mockrun.hpp

#ifndef WINDOW_TD_MOCKRUN_HPP
#define WINDOW_TD_MOCKRUN_HPP

/* window_td_morkrun.hpp
 *
 * window_td_morkrun_t:          checking window_td with ESP-IDF or ESP32 Arduino
 * window_td_test_with_octave_t: checking window_td with Octave on a PC
 */

#include <iostream>
#include <fstream>
#include <complex>
#include <random>
#include <thread>
#include <cassert>

#include "pocketfft.h"
#include "window_td.h"

// ------------------------------------------------------------------------------------------------

class window_td_mockrun_t
{
public:

    window_td_mockrun_t(int n) : max_length(n) {epsilon=1e-15;};

    bool run(int maxlength = 20)
    {
        using namespace window_td;

        printf("Checking windows in symmetric and periodic forms\n"
               "with the maximum lengths of %d:\n", maxlength);

        printf("bartlett ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = bartlett(k, symmetric);
            w2 = bartlett(k-1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("blackman ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = blackman(k, symmetric);
            w2 = blackman(k-1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("hamming ...         "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = hamming(k+1, symmetric);
            w2 = hamming(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("hann ...            "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = hann(k+1, symmetric);
            w2 = hann(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        // --------

        printf("barthannwin ...     "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = barthannwin(k+1, symmetric);
            w2 = barthannwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("blackmanharris ...  "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = blackmanharris(k+1, symmetric);
            w2 = blackmanharris(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("blackmannuttall ... "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = blackmannuttall(k+1, symmetric);
            w2 = blackmannuttall(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("bohmanwin ...       "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = bohmanwin(k+1, symmetric);
            w2 = bohmanwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("chebwin ...         "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = chebwin(k+1, 100, symmetric);
            w2 = chebwin(k, 100, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("flattopwin ...      "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = flattopwin(k+1, symmetric);
            w2 = flattopwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("gaussian ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = gaussian(k+1, 1, symmetric);
            w2 = gaussian(k, 1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("gausswin ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = gausswin(k+1, 2.5, symmetric);
            w2 = gausswin(k, 2.5, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("kaiser ...          "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = kaiser(k+1, 0.5, symmetric);
            w2 = kaiser(k, 0.5, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("nuttallwin ...      "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = nuttallwin(k+1, symmetric);
            w2 = nuttallwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("parzenwin ...       "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = parzenwin(k+1, symmetric);
            w2 = parzenwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("triang ...          "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = triang(k+1, symmetric);
            w2 = triang(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("tukeywin ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = tukeywin(k+1, 0.5, symmetric);
            w2 = tukeywin(k, 0.5, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("ultrwin ...         "); fflush(stdout);
        for(int k = 8; k <= maxlength; ++k){
            w1 = ultrwin(k+1, 1, 2.5, beta, 0, NULL, symmetric);
            w2 = ultrwin(k, 1, 2.5, beta, 0, NULL, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("welchwin ...        "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = welchwin(k+1, symmetric);
            w2 = welchwin(k, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        // --------

        printf("lanczoswin ...      "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = lanczoswin(k+1, 1, symmetric);
            w2 = lanczoswin(k, 1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("plancktaper ...     "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = plancktaper(k+1, 0.1, symmetric);
            w2 = plancktaper(k, 0.1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("sinwin ...          "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = sinwin(k+1, 1, symmetric);
            w2 = sinwin(k, 1, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("taylorwin ...       "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = taylorwin(k+1, 4, 100, true, symmetric);
            w2 = taylorwin(k, 4, 100, true, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        printf("sprwin ...          "); fflush(stdout);
        for(int k = 3; k <= maxlength; ++k){
            w1 = sprwin(k+1, 5.0, sinwin, 1.0, symmetric);
            w2 = sprwin(k, 5.0, sinwin, 1.0, periodic);
            if(!check_results(w1, w2)) return false;
        }
        print_sym_rms(w1);

        return true;
    }

protected:

    double epsilon;
    int max_length;
    std::vector <double> w1, w2;

    bool check_results(const std::vector <double> &vs, const std::vector <double> &vp)
    {
        int ns = vs.size(), np = vp.size();

        if((ns-np) != 1){
            printf("Incorrect vector length -- %d and %d!\n", ns, np);
            return false;
        }

        if(abs(vs[0] - vs[ns-1]) > epsilon){
            puts("The given vector is not symmetric!");
            return false;
        }

        for(int i = 0; i < np; ++i){
            if(abs(vs[i] - vp[i]) > epsilon){
                printf("vp.size = %d, which is inconsistent with vs[%d]'s.\n", np, i);
                for(int j = 0; j < np; ++j){
                    printf("%e\t%e\n", vs[j], vp[j]);
                }
                return false;
            }
        }

        return true;
    }

    // This RMS value can be verified using the Octave function rms, such as rms(chebwin(n))
    void print_sym_rms(const std::vector <double> &v)
    {
        double rms = 0;
        for(auto val : v){
            rms += val*val;
        }
        printf("RMS for symmetric %d: %f\n", int(v.size()), sqrt(rms/v.size()));
    }
};

// ------------------------------------------------------------------------------------------------

class window_td_test_with_octave_t
{
public:

    window_td_test_with_octave_t(size_t _n = 1024)
    {
        n = _n;
    }

    void compare_windows()
    {
        using namespace window_td;

        // window functions in Octave /m/signal
        po_data(0, bartlett(n), "w_bartlett_s", "bartlett(%d)");
        po_data(1, bartlett(n, periodic), "w_bartlett_p", "bartlett(%d)");
        po_data(0, blackman(n), "w_blackman_s", "blackman(%d)");
        po_data(1, blackman(n, periodic), "w_blackman_p", "blackman(%d)");
        po_data(0, hamming(n), "w_hamming_s", "hamming(%d)");
        po_data(1, hamming(n, periodic), "w_hamming_p", "hamming(%d)");
        po_data(0, hann(n), "w_hann_s", "hann(%d)");
        po_data(1, hann(n, periodic), "w_hann_p", "hann(%d)");

        // window functions in Octave /packages/signal-1.3.2
        po_data(0, barthannwin(n), "w_barthannwin_s", "barthannwin(%d)");
        po_data(1, barthannwin(n, periodic), "w_barthannwin_p", "barthannwin(%d)");
        po_data(0, blackmanharris(n), "w_blackmanharris_s", "blackmanharris(%d)");
        po_data(1, blackmanharris(n, periodic), "w_blackmanharris_p", "blackmanharris(%d)");
        po_data(0, blackmannuttall(n), "w_blackmannuttall_s", "blackmannuttall(%d)");
        po_data(1, blackmannuttall(n, periodic), "w_blackmannuttall_p", "blackmannuttall(%d)");
        po_data(0, bohmanwin(n), "w_bohmanwin_s", "bohmanwin(%d)");
        po_data(1, bohmanwin(n, periodic), "w_bohmanwin_p", "bohmanwin(%d)");
        po_data(0, boxcar(n), "w_boxcar", "boxcar(%d)");
        po_data(0, chebwin(n), "w_chebwin_s", "chebwin(%d)");
        po_data(1, chebwin(n, 100, periodic), "w_chebwin_p", "chebwin(%d)");
        po_data(0, flattopwin(n), "w_flattopwin_s", "flattopwin(%d)");
        po_data(1, flattopwin(n, periodic), "w_flattopwin_p", "flattopwin(%d)");
        po_data(0, gaussian(n), "w_gaussian_s", "gaussian(%d)");
        po_data(1, gaussian(n, 1, periodic), "w_gaussian_p", "gaussian(%d)");
        po_data(0, gausswin(n), "w_gausswin_s", "gausswin(%d)");
        po_data(1, gausswin(n, 2.5, periodic), "w_gausswin_p", "gausswin(%d)");
        po_data(0, hann(n), "w_hann_s", "hann(%d)");
        po_data(1, hann(n, periodic), "w_hann_p", "hann(%d)");
        po_data(0, kaiser(n), "w_kaiser_s", "kaiser(%d)");
        po_data(1, kaiser(n, 0.5, periodic), "w_kaiser_p", "kaiser(%d)");
        po_data(0, nuttallwin(n), "w_nuttallwin_s", "nuttallwin(%d)");
        po_data(1, nuttallwin(n, periodic), "w_nuttallwin_p", "nuttallwin(%d)");
        po_data(0, parzenwin(n), "w_parzenwin_s", "parzenwin(%d)");
        po_data(1, parzenwin(n, periodic), "w_parzenwin_p", "parzenwin(%d)");
        po_data(0, rectwin(n), "w_rectwin", "rectwin(%d)");
        po_data(0, triang(n), "w_triang_s", "triang(%d)");
        po_data(1, triang(n, periodic), "w_triang_p", "triang(%d)");
        po_data(0, tukeywin(n), "w_tukeywin_s", "tukeywin(%d)");
        po_data(1, tukeywin(n, 0.5, periodic), "w_tukeywin_p", "tukeywin(%d)");
        po_data(0, ultrwin(n, -1, 100, latt, 0, NULL), "w_ultrwin_s", "ultrwin(%d, -1, 100, 'l')");
        po_data(1, ultrwin(n, -1, 100, latt, 0, NULL, periodic), "w_ultrwin_p",
                                                                 "ultrwin(%d, -1, 100, 'l')");
        po_data(0, welchwin(n), "w_welchwin_s", "welchwin(%d)");
        po_data(1, welchwin(n, periodic), "w_welchwin_p", "welchwin(%d)");

        po_mfile(__func__, -200);
    }

    void check_blackman()
    {
        using namespace window_td;

        static const std::vector <bool> optimal = {false, true};
        for(size_t i = 0; i < optimal.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(blackman(n, symmetric, optimal[i]), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_hamming()
    {
        using namespace window_td;

        struct opts_t{bool precise; bool optimal;};
        static const std::vector <opts_t> opts = {
            {false, false},     // precise = false, optimal = false
            {true, false},      // precise = true, optimal = false
            {false, true},      // precise = false, optimal = true
        };

        for(size_t i = 0; i < opts.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(hamming(n, symmetric, opts[i].precise, opts[i].optimal), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_nuttallwin_3_terms()
    {
        using namespace window_td;

        for(size_t i = NW3; i <= NW3B; ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(nuttallwin(n, symmetric, nw_index(i)), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_nuttallwin_4_terms()
    {
        using namespace window_td;

        for(size_t i = NW4; i <= NW4C; ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(nuttallwin(n, symmetric, nw_index(i)), data_fn);
        }
        pf_mfile(__func__, -200, n);
    }

    void check_flattopwin()
    {
        using namespace window_td;

        static const std::vector <bool> optimal = {false, true};
        for(size_t i = 0; i < optimal.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(flattopwin(n, symmetric, optimal[i]), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_gaussian()
    {
        using namespace window_td;

        static const std::vector <double> a = {0.02, 0.04, 0.06, 0.08, 0.1};
        for(size_t i = 0; i < a.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(gaussian(n, a[i]), data_fn);
        }
        pf_mfile(__func__, -170, n);
    }

    void check_gausswin()
    {
        using namespace window_td;

        static const std::vector <double> a = {2, 2.7, 3.2, 3.7, 4.2};
        for(size_t i = 0; i < a.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(gausswin(n, a[i]), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_tukeywin()
    {
        using namespace window_td;

        static const std::vector <double> r = {0.1, 0.3, 0.5, 0.7, 0.9};
        for(size_t i = 0; i < r.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(tukeywin(n, r[i]), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_kaiser()
    {
        using namespace window_td;

        static const std::vector <double> beta = {3.4, 5.65, 7.85, 10, 12.3};
        for(size_t i = 0; i < beta.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, i);
            pf_data(kaiser(n, beta[i]), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_chebwin()
    {
        using namespace window_td;

        for(size_t at = 40; at <= 120; at += 20){
            std::string data_fn = snformat("%s_%d_dB", __func__, at);
            pf_data(chebwin(n, at), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_ultrwin()
    {
        using namespace window_td;

        double xmu;

        // latt
        static const std::vector <double> att_db = {40, 60, 80, 100, 120};
        for(size_t i = 0; i < att_db.size(); ++i){
            std::string data_fn = snformat("%s_latt_%d", __func__, int(i));
            pf_data(ultrwin(n, -1, att_db[i], latt, 0, &xmu), data_fn);
        }
        pf_mfile(snformat("%s_latt", __func__).c_str(), -140, n);

        // att
        for(size_t i = 0; i < att_db.size(); ++i){
            std::string data_fn = snformat("%s_att_%d", __func__, int(i));
            pf_data(ultrwin(n, -1, att_db[i], att, 0, &xmu), data_fn);
        }
        pf_mfile(snformat("%s_att", __func__).c_str(), -140, n);

        // beta, generating windows equavalent to Dolph-Chebyshev sequences
        static const std::vector <double> betaval = {1.77, 2.47, 3.2, 3.93, 4.66};
        for(size_t i = 0; i < betaval.size(); ++i){
            std::string data_fn = snformat("%s_cheb_%d", __func__, int(i));
            pf_data(ultrwin(n, 0, betaval[i], beta, 0, &xmu), data_fn);
        }
        pf_mfile(snformat("%s_cheb", __func__).c_str(), -140, n);

        // beta, generating windows equavalent to Tapio Saramäki sequences
        mlines.clear();
        for(size_t i = 0; i < betaval.size(); ++i){
            std::string data_fn = snformat("%s_sara_%d", __func__, int(i));
            pf_data(ultrwin(n, 1, betaval[i], beta, 0, &xmu), data_fn);
        }
        pf_mfile(snformat("%s_sara", __func__).c_str(), -140, n);
    }

    void check_lanczoswin()
    {
        using namespace window_td;

        for(size_t i = 1; i <= 5; i ++){
            std::string data_fn = snformat("%s_%d", __func__, int(i));
            pf_data(lanczoswin(n, i), data_fn);
        }
        pf_mfile(__func__, -300, n);
    }

    void check_plancktaper()
    {
        using namespace window_td;

        for(size_t i = 1; i <= 5; i ++){
            std::string data_fn = snformat("%s_%d", __func__, int(i));
            pf_data(plancktaper(n, i/10.0), data_fn);
        }
        pf_mfile(__func__, -200, n);
    }

    void check_sinwin()
    {
        using namespace window_td;

        static const std::vector <double> alpha = {0.5, 1, 2, 4, 8};
        for(size_t i = 0; i < alpha.size(); i ++){
            std::string data_fn = snformat("%s_%d", __func__, int(i));
            pf_data(sinwin(n, alpha[i]), data_fn);
        }
        pf_mfile(__func__, -300, n);
    }

    void check_taylorwin()
    {
        using namespace window_td;

        for(size_t nbar = 2; nbar <= 6; ++nbar){
            std::string data_fn = snformat("%s_%dbars", __func__, nbar);
            pf_data(taylorwin(n, nbar), data_fn);
        }
        pf_mfile(__func__, -140, n);
    }

    void check_sprwin()
    {
        using namespace window_td;

        static const std::vector <double> b = {2, 8, 24, 100};
        for(size_t i = 0; i < b.size(); ++i){
            std::string data_fn = snformat("%s_%d", __func__, int(i));
            pf_data(sprwin(n, b[i], kaiser, 2.0, symmetric), data_fn);
        }
        pf_mfile(__func__, -200, n);
    }

    void check_HFTxx_D()
    {
        using namespace window_td;

        // Heinzel described a set of "new flat-top windows" in his paper on the net. These cosine-
        // sum windows offer better amplitude flatness and sidelobe attenuations than classic flat-
        // top windows, and are therefore more suitable for precise signal analysis.
        auto HFTxx_D = [] (size_t n, size_t w_type, form_opt opt = symmetric) {
            static const std::vector <double> a[5] = {  // HFT70, HFT95, HFT90D, HFT116D, HFT144D
                {1, 1.90796, 1.07349, 0.18199},
                {1, 1.9383379, 1.3045202, 0.4028270, 0.0350665},
                {1, 1.942604, 1.340318, 0.440811, 0.043097},
                {1, 1.9575375, 1.4780705, 0.6367431, 0.1228389, 0.0066288},
                {1, 1.96760033, 1.57983607, 0.81123644, 0.22583558, 0.02773848, 0.00090360},
            };

            return gencossum(n, a[w_type], opt);
        };

        for(size_t i = 0; i < 5; ++i){
            std::string data_fn = snformat("%s_%d", __func__, int(i));
            pf_data(HFTxx_D(n, i), data_fn);
        }
        pf_mfile(__func__, -200, n);
    }

protected:

    size_t n;
    std::vector <std::string> mlines;

    template <typename... Args>
    std::string snformat(const char* format, Args... args)
    {
        int length = std::snprintf(nullptr, 0, format, args...);
        assert(length >= 0);

        std::vector <char> buf(length+1);
        std::snprintf(&buf[0], length+1, format, args...);

        return std::string(&buf[0], length);
    }

    void po_data(int p, const std::vector <double> &y,  const std::string &fn, const char *of)
    {
        std::ofstream fo;
        std::string txt_file_name = fn + ".txt";

        fo.open(txt_file_name);
        for(auto val : y){
            fo << snformat("%.16e\n", val);
        }
        fo.close();

        mlines.push_back("");
        mlines.push_back("load " + txt_file_name);

        std::string ref, statement, line;
        if(p){
            ref = snformat("%s(1:%d)", of, n);
            ref = snformat(ref.c_str(), n+1);
        }else{
            ref = snformat("%s", of);
            ref = snformat(ref.c_str(), n);
        }
        statement = snformat("20 * log10(rms(%s-ref) / rms(ref))", fn.c_str());
        line = snformat("20*log10(rms(%s-%s)/rms(%s))", fn.c_str(), ref.c_str(), ref.c_str());

        mlines.push_back("ref = " + ref + ";");
        mlines.push_back("ans = " + statement + ";");
        mlines.push_back("if ans > th");
        mlines.push_back("    printf(\"" + line + " = \\n%f dB\\n\", ans)");
        mlines.push_back("    printf(\"    Warning: the accumulated error exceeds %f dB.\\n\", th)");
        mlines.push_back("endif");
    }

    void po_mfile(const std::string &fn, double err_threshold)
    {
        std::string mfile_name = fn + ".m";

        std::ofstream fo;
        fo.open(mfile_name);

        fo << "clear\n";
        fo << snformat("th = %f;\n", err_threshold);
        fo << "pkg load signal\n";
        fo << "printf(\"Checking window_td functions with Octave...\\n\")\n";
        for(auto line : mlines){
            fo << line << std::endl;
        }
        mlines.clear();
        fo << "printf(\"Checked.\\n\")\n";

        fo.close();
    }

    void pf_data(const std::vector< double >&y,  const std::string &fn)
    {
        std::ofstream fo;
        std::string txt_file_name = fn + ".txt";

        fo.open(txt_file_name);
        for(auto val : y){
            fo << snformat("%.16e\n", val);
        }
        fo.close();

        mlines.push_back(fn);
    }

    void pf_mfile(const std::string &fn, int slobe, size_t npoints)
    {
        std::ofstream fo;
        std::string mfile_name = fn + ".m";
        fo.open(mfile_name);

        fo << "clear\n";
        fo << "c='krgbm'; clf; figure(1)\n";
        fo << "subplot(2, 1, 1)\n";
        for(size_t i = 0; i < mlines.size(); ++i){
            fo << "load " + mlines[i] + ".txt;\n";
            fo << "[W,f]=freqz(" + mlines[i] + ");\n";
            fo << "plot(f/pi, real(20*log10(W/max(W) + 1e-50)), c(1+mod(";
            fo << i << ", length(c)))); hold on\n";
        }
        fo << "grid on; axis([0 1 " << slobe << " 0]); hold off\n";

        fo << "subplot(2, 1, 2)\n";
        fo << "ymin = 0;\n";
        for(size_t i = 0; i < mlines.size(); ++i){
            fo << mlines[i] + " = " + mlines[i] + "/max(" + mlines[i] + ");\n";
            fo << "cmin = min(" + mlines[i] + ");\n";
            fo << "ymin = min(ymin, cmin);\n";
            fo << "plot(" + mlines[i] + ", c(1+mod(";
            fo << i << ", length(c)))); hold on\n";
        }
        mlines.clear();
        fo << "grid on; axis([0 " << npoints << " ymin 1]); hold off\n";

        fo.close();
    }
};

#endif // WINDOW_TD_MOCKRUN_HPP

window_td_ultrwin.cpp

/* window_td_ultrwin.cpp v1.3
 *
 * Copyright (C) 2018, 2021, 2023 BD1ES <bd1es@hotmail.com>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; see the file COPYING.  If not, see
 * <https://www.gnu.org/licenses/>.
 */

//-------------------------------------------------------------------------------------------------
// The program below is part of __ultrwin__.cc ----------------------------------------------------
// See: https://octave.sourceforge.io/signal/index.html
// It has been modified a bit to compile with window_td. See lines 48-51.

/*

Copyright (C) 2013 Rob Sykes <robs@users.sourceforge.net>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, see
<https://www.gnu.org/licenses/>.

*/

#include <cmath>
#include <cstdlib>
#include <cstdio>
/*
#include <octave/oct.h>
#include <octave/lo-specfun.h>
*/

#if !defined M_PI
#define M_PI 3.14159265358979323846
#endif

#if defined (__cplusplus) && __cplusplus > 201402L
#  define ATTR_FALLTHROUGH [[fallthrough]]
#elif defined (__GNUC__) && __GNUC__ >= 7
#  define ATTR_FALLTHROUGH __attribute__ ((__fallthrough__))
#else
#  define ATTR_FALLTHROUGH ((void) 0)
#endif

#if DEBUG_ULTRWIN
#define SHOW1(x) fprintf(stderr, "%c%+.3e", " 1"[(x) > .5], (x) - ((x) > .5))
#endif

#define EPSILON     (1./0x100000) /* For caller rounding error. */
#define BETA_MAX    (12*(1+EPSILON))
#define XMU_MIN     (.99*(1-EPSILON))


static double *
ultraspherical_win(int n, double mu, double xmu)
{
  double * w = NULL;
  int bad = n < 1 || xmu < XMU_MIN || (!mu && xmu == 1) ||
    (n > BETA_MAX * 2 && xmu * cos(M_PI * BETA_MAX / n) > 1);
  if (!bad && (w = (double *)malloc(sizeof(*w)* n))) {
    int i, j, k, l = 0, m = (n + 1) / 2, met;
    double * divs = w + m - 1, c = 1 - 1 / (xmu * xmu), t, u, v[64], vp, s;
    for (i = 0; i < (int)(sizeof(v) / sizeof(v[0])); v[i++] = 0);
    if (n > 1) for (i = 0; i < m; l = j - (j <= i++)) {
      vp = *v, s = *v = i? (*v + v[1]) * mu * (divs[i] = 1./i) : 1;
      for (met = 0, j = 1, u = 1; ; ++l, v[l] = vp * (i - l) / (mu + l - 1)) {
        #define _ t = v[j], v[j] += vp, vp = t, t = s, s += \
            v[j] * (u *= c * (n - i - j) * divs[j]), met = s && s == t, ++j,
        for (k = ((l-j+1) & ~7) + j; j < k && !met; _ _ _ _ _ _ _ _ (void)0);
        for (; j <= l && !met; _ (void)0);
        #undef _
        if (met || !(j <= i)) break;
      }
      w[i] = s / (n - i - 1);
    }
    else w[0] = 1;
    u = 1 / w[i = m - 1], w[i] = 1;
    for (--i ; i >= 0; u *= (n - 2 - i + mu) / (n - 2 - i), w[i] *= u, --i);
    for (i = 0; i < m; w[n - 1 - i] = w[i], ++i);
  }
  return w;
}


typedef struct {double f, fp;} uspv_t;


static uspv_t
ultraspherical_polyval(int n, double mu, double x, double const *divs)
{
  double fp = n > 0? 2 * x * mu : 1, fpp = 1, f;
  uspv_t result;
  int i, k;
  #define _ f = (2*x*(i+mu)*fp - (i+2*mu-1)*fpp) * divs[i+1], fpp=fp, fp=f, ++i,
  for (i = 1, k = i + ((n - i) & ~7); i < k; _ _ _ _ _ _ _ _ (void)0);
  for (; i < n; _ (void)0);
  #undef _
  result.f = fp, result.fp = fpp;
  return result;
}


#define MU_EPSILON  (1./0x4000)
#define EQ(mu,x)    (fabs((mu)-(x)) < MU_EPSILON)


static uspv_t
ultraspherical_polyval2(     /* With non-+ve integer protection */
    int n, double mu, double x, double const * divs)
{
  int sign = (~(int)floor(mu) & ~(~2u/2))? 1:-1; /* -ve if floor(mu) <0 & odd */
  uspv_t r;
  if (mu < MU_EPSILON && EQ(mu,(int)mu))
    mu = floor(mu + .5) + MU_EPSILON * ((int)mu > mu? -1:1);
  r = ultraspherical_polyval(n, mu, x, divs);
  r.f *= sign, r.fp *= sign;
  return r;
}


static double
find_zero(int n, double mu, int l, double extremum_mag, double ripple_ratio,
          double lower_bound, double const *divs)
{
  double dx, x0, t, x, epsilon = 1e-10, one_over_deriv, target = 0;
  int i, met = 0;
  if (!divs)
    return 0;
  if (!l) {
    double r = ripple_ratio;   /* FIXME: factor in weighted extremum_mag here */
    x = r > 1 ? cosh(acosh(r) / n) : cos(acos(r) / n); /* invert chebpoly-1st */
    x0 = x *= lower_bound / cos(M_PI * .5 / n) + epsilon;
    target = log(extremum_mag * ripple_ratio);
  }
  else {
    double cheb1 = cos(M_PI * (l - .5) / n), cheb2 = cos(M_PI * l / (n + 1));
    if (mu < 1 - l && EQ((int)(mu+.5),mu+.5)) x = met = 1;
    else if (EQ(mu,0)) x = cheb1, met = 1;               /* chebpoly-1st-kind */
    else if (EQ(mu,1)) x = cheb2, met = 1;               /* chebpoly-2nd-kind */
    else x = (cheb1 * cheb2) / (mu * cheb1 + (1 - mu) * cheb2);
    x0 = x;
  }
  for (i = 0; i < 24 && !met; ++i, met = fabs(dx) < epsilon) {/*Newton-Raphson*/
    uspv_t r = ultraspherical_polyval2(n, mu, x, divs);
    if (!(t = ((2*mu + n-1) * r.fp - n*x * r.f)))         /* Fail if div by 0 */
      break;
    one_over_deriv = (1 - x*x) / t;    /* N-R slow for deriv~=1, so take log: */
    if (!l) {                               /* d/dx(f(g(x))) = f'(g(x)).g'(x) */
      one_over_deriv *= r.f;                             /* d/dx(log x) = 1/x */
      if (r.f <= 0)                                 /* Fail if log of non-+ve */
        break;
      if (x + (dx = (target - log(r.f)) * one_over_deriv) <= lower_bound)
        dx = (lower_bound - x) * .875;
      x += dx;
    }
    else x += dx = -r.f * one_over_deriv;
#if DEBUG_ULTRWIN
    fprintf(stderr, "1/deriv=%9.2e dx=%9.2e x=", one_over_deriv, dx);
    SHOW1(x); fprintf(stderr, "\n");
#endif
  }
#if DEBUG_ULTRWIN
  fprintf(stderr, "find_zero(n=%i mu=%g l=%i target=%g r=%g x0=",
      n, mu, l, target, ripple_ratio);
  SHOW1(x0); fprintf(stderr, ") %s ", met? "converged to" : "FAILED at");
  SHOW1(x); fprintf(stderr, " in %i iterations\n", i);
#else
  static_cast<void>(x0);
#endif
  return met? x : 0;
}


static double *
make_divs(int n, double **divs)
{
  int i;
  if (!*divs) {
    *divs = (double *)malloc(n * sizeof(**divs));
    if (*divs)
      for (i = 0; i < n; (*divs)[i] = 1./(i+1), ++i);
  }
  return *divs? *divs - 1 : 0;
}


#define DIVS make_divs(n, &divs)


typedef enum {uswpt_Xmu, uswpt_Beta, uswpt_AttFirst, uswpt_AttLast} uswpt_t;


double *
ultraspherical_window(int n, double mu, double par, uswpt_t type, int even_norm,
                      double *xmu_)
{
  double * w = 0, xmu = 0, * divs = 0, last_extremum_pos = 0;

  if (n > 0 && fabs(mu) <= (8*(1+EPSILON))) switch (type) {
    case uswpt_Beta:
      xmu = mu == 1 && par == 1? 1 : par < .5 || par > BETA_MAX? 0 :
        find_zero(n-1, mu, 1, 0, 0, 0, DIVS) / cos(M_PI * par / n);
      break;

    case uswpt_AttFirst:
      if (par < 0) break;
      ATTR_FALLTHROUGH;

    case uswpt_AttLast:
      if (type == uswpt_AttLast && mu >= 0 && par < 0);
      else if (!EQ(mu,0)) {
        int extremum_num =
          type == uswpt_AttLast? (int)((n-2)/2 +.5) : 1 + EQ(mu,-1.5);
        double extremum_pos =
          find_zero(n-2, mu+1, extremum_num, 0, 0, 0, DIVS);
        double extremum_mag = !extremum_pos? 0 :
          fabs(ultraspherical_polyval2(n-1, mu, extremum_pos, DIVS).f);
        double xmu_lower_bound = !extremum_mag? 0 :
          find_zero(n-1, mu, 1, 0, 0, 0, DIVS); /* 1st null */
        xmu = !xmu_lower_bound? 0 : find_zero(
            n-1, mu, 0, extremum_mag, pow(10, par/20), xmu_lower_bound, DIVS);
        last_extremum_pos =
          type == uswpt_AttLast? extremum_pos : last_extremum_pos;
      }
      else xmu = cosh(acosh(pow(10, par/20))/(n-1)); /* Cheby 1st kind */
      break;

    default: case uswpt_Xmu: xmu = par; break;
  }
#if DEBUG_ULTRWIN
  fprintf(stderr, "n=%i mu=%.3f xmu=%.16g\n", n, mu, xmu);
#endif

  if (xmu > 0)
    w = ultraspherical_win(n, mu, xmu);

  if (w && (~n & !!even_norm) && n > 2 && !(mu == 1 && xmu == 1)) {
    int i = n / 2 - 1, j = 1;
    double * d = DIVS, t = 0, s = -1, x = even_norm == 1? 0 : last_extremum_pos?
      last_extremum_pos : find_zero(n-2, mu+1, i, 0, 0, 0, d);
    x = x? M_PI/2 - acos(x/xmu) : 0;
    for (; i >= 0; t += w[i] * d[j] * (s=-s) * (x?cos(j*x):1), --i, j += 2);
    for (t = M_PI/4 / t, i = 0; t < 1 && i < n; w[i] *= t, ++i);
#if DEBUG_ULTRWIN
    fprintf(stderr, "%snorm DFT(w.sinc πx) @ %g %.16g\n", t<1? "":"NO ", 2*x,t);
#endif
  }
  free(divs);
  if (xmu_)
    *xmu_ = xmu;
  return w;
}

// The program above is part of __ultrwin__.cc ----------------------------------------------------
//-------------------------------------------------------------------------------------------------


#include <window_td.h>

/* window_td::ultrwin
 *
 * As mentioned in "ultrawin.m"
 *
 * -- Function File: [W, XMU] = ultrwin (M, MU, BETA)
 * -- Function File: [W, XMU] = ultrwin (M, MU, ATT, "att")
 * -- Function File: [W, XMU] = ultrwin (M, MU, LATT, "latt")
 * -- Function File: W = ultrwin (M, MU, XMU, "xmu")
 * Return the coefficients of an Ultraspherical window of length M.
 * The parameter MU controls the window's Fourier transform's
 * side-lobe to side-lobe ratio, and the third given parameter
 * controls the transform's main-lobe width/side-lobe-ratio; normalize
 * W such that the central coefficient(s) value is unitary.
 *
 * By default, the third parameter is BETA, which sets the main lobe
 * width to BETA times that of a rectangular window.  Alternatively,
 * giving ATT or LATT sets the ripple ratio at the first or last
 * side-lobe respectively, or giving XMU sets the (un-normalized)
 * window's Fourier transform according to its canonical definition:
 *
 *              (MU)
 *      W(k) = C   [ XMU cos(pi k/M) ],  k = 0, 1, ..., M-1,
 *              M-1
 *
 * where C is the Ultraspherical (a.k.a.  Gegenbauer) polynomial,
 * which can be defined using the recurrence relationship:
 *
 *       (l)    1                  (l)                    (l)
 *      C (x) = - [ 2x(m + l - 1) C   (x) - (m + 2l - 2) C   (x) ]
 *       m      m                  m-1                    m-2
 *
 *                                 (l)        (l)
 *      for m an integer > 1, and C (x) = 1, C (x) = 2lx.
 *                                 0          1
 *
 * For given BETA, ATT, or LATT, the corresponding (determined) value
 * of XMU is also returned.
 *
 * The Dolph-Chebyshev and Saramaki windows are special cases of the
 * Ultraspherical window, with MU set to 0 and 1 respectively.  Note
 * that when not giving XMU, stability issues may occur with MU <=
 * -1.5.  For further information about the window, see
 *
 *    * Kabal, P., 2009: Time Windows for Linear Prediction of Speech.
 *      Technical Report, Dept.  Elec.  & Comp.  Eng., McGill
 *      University.
 *    * Bergen, S., Antoniou, A., 2004: Design of Ultraspherical
 *      Window Functions with Prescribed Spectral Characteristics.
 *      Proc.  JASP, 13/13, pp.  2053-2065.
 *    * Streit, R., 1984: A two-parameter family of weights for
 *      nonrecursive digital filters and antennas.  Trans.  ASSP, 32,
 *      pp.  108-118.
 */

std::vector <double> window_td::ultrwin(size_t n, double mu, double par, ultr_par type,
                                        int even_norm, double *xmu_, form_opt fopt)
{
    if(n < 8){
        return boxcar(n);
    }

    // For enum {uswpt_Xmu, uswpt_Beta, uswpt_AttFirst, uswpt_AttLast} uswpt_t;
    uswpt_t par_type;
    switch(type){
    case xmu:
        par_type = uswpt_Xmu;
        break;
    case beta:
        par_type = uswpt_Beta;
        break;
    case att:
        par_type = uswpt_AttFirst;
        break;
    case latt:
        par_type = uswpt_AttLast;
        break;
    default:
        return boxcar(0);
    }

    std::vector <double> y(n);
    int N = fopt == symmetric ? n : n+1;

    double xmu;
    double *w = ultraspherical_window (N, mu, par, par_type, even_norm, &xmu);
    if(!w){
#if DEBUG_ULTRWIN
        fprintf(stderr, "ultrwin: parameter(s) out of range.\n");
#endif
        y = boxcar(0);
    }else{
        y.assign(w, w+n);
        free(w);
    }
    if(xmu_) *xmu_ = xmu;

    return y;
}

wnd_main.cpp

#include <iostream>
#include <cstdio>

#include <window_td_mockrun.hpp>

int main()
{
    puts("Testing window_td...\n");
    
    window_td_mockrun_t mock(256);
    if(!mock.run(128)){
        puts("Error occurred while testing!");
        return -1;
    }
    
    window_td_test_with_octave_t oct;
    oct.compare_windows();
    oct.check_blackman();
    oct.check_hamming();
    oct.check_nuttallwin_3_terms();
    oct.check_nuttallwin_4_terms();
    oct.check_flattopwin();
    oct.check_gaussian();
    oct.check_gausswin();
    oct.check_tukeywin();
    oct.check_kaiser();
    oct.check_chebwin();
    oct.check_ultrwin();
    oct.check_lanczoswin();
    oct.check_plancktaper();
    oct.check_sinwin();
    oct.check_taylorwin();
    oct.check_sprwin();
    oct.check_HFTxx_D();
    
    puts("Completed.");
    
    return 0;
}

wnd_sketch_sep22a.ino

#include "window_td_mockrun.hpp"

void setup() {
  // put your setup code here, to run once:
  
  Serial.begin(115200);
  while(!Serial);
  Serial.println("Ready");
}

void loop() {
  // put your main code here, to run repeatedly:
  
  puts("Testing window_td ...\n");
  window_td_mockrun_t mock(256);
  if(!mock.run(8)){
      puts("Error occurred while testing!");
      return;
  }
  puts("Completed.");

  // Sleep and wait for a system reset.
  while(1) sleep(100);
}