apulSoft Blog

Aug 25th, 2025 - cpp math approximation

1D/2D Chebyshev Approximation in C++

Motivation

In audio DSP, smooth functions appear everywhere. Any bounded smooth function can be approximated with Chebyshev polynomials. I've put together some C++ code that creates such approximations: It takes a lambda function, builds the approximation, and generates ready-to-use C++ code to be dropped in as a replacement for the original lambda. This has allowed me to use polynomial approximations for sections of code that used iterative solvers before.

Working on a special type of filter, I needed the same approach in two dimensions. I realized it is possible to Chebyshev-approximate vectors of Chebyshev coefficients, since their magnitudes are usually of similar size. To make it work with coefficient vectors, the code makes use of templates and operator overloading.

1D Example

To approximate $f(x)=cos(x)$ between 0 and $pi/2$:

plap::Chebyshev_Approximation_1D<double, 7> apx;
auto f = [](auto x){ return std::cos(x); };
auto est_err = apx.createForFunctionRange({0.0, M_PI / 2}, f);

When creating the approximation, the code first builds a higher-order approximation, then truncates to the desired order. Usually, this reduces the error slightly. By default, 4 extra coefficients are used, but this can be adjusted in the call.

apx.eval(x) can now be used to approximate cos(x). A direct ready-to-use constructor can be generated with: std::cout << apx << "\n";

auto apx = plap::Chebyshev_Approximation_1D<double, 7>({0, 1.57079632679489656}, 
    {{1.20438940251109283, -0.513625166679106959, -0.103546344262963899, 
      0.0137320342343593502, 0.00135866983808981094, -0.000107263094406512534, 
      -7.04629679416041767e-06}});
auto y = apx.eval(0.7); // = 0.7648419..  (cos(0.7) = 0.7648421..)

The evaluation routine is based on the Clenshaw algorithm. It is fast, accurate, compiler-optimization-friendly, and avoids temporary variables inside the loops. This is important for the 2D case, where extra copies would be expensive.

Using polyCoeffs(), the 1D class can also calculate the coefficients of an equivalent polynomial, which can be evaluated even faster with something like Estrin's scheme. However, if there are many coefficients (≥12) or if precision needs to be as high as possible, the Clenshaw algorithm performs better.

2D Example

To approximate over the xy-area 0|0 ↔ 1|1:

$ f(x,y)=cos(3.5x)*sin(2.8y) + 0.5xy $

auto f = [](auto x, auto y){ 
    return std::cos(3.5 * x) * std::sin(2.8 * y) + 0.5 * x * y; 
};
plap::Chebyshev_Approximation_2D<double, 5, 5> apx2d;
apx2d.createForFunctionRange({0.0, 1.0}, {0.0, 1.0}, f);
cout << apx2d << "\n";

Result:

auto apx2d = plap::Chebyshev_Approximation_2D<double, 5, 5>({0, 1}, {0, 1}, 
    {{{{0.353022410913424, -1.02555847352952, 0.117095016978286, 
        0.201947768190967, -0.00832974686170527}}, 
    {{0.225763665756079, -0.0853372473603274, 0.0193087530379324, 
        0.0333008157237305, -0.0013735599444935}}, 
    {{0.0537644789646134, 0.466599956798963, -0.042833418457298, 
        -0.0738725992331952, 0.00304702575889243}}, 
    {{0.00225829803570464, 0.0195988464166555, -0.00179915488120506, 
        -0.00310291011748709, 0.00012798584527486}}, 
    {{-0.00234987563220483, -0.0203936109785677, 0.00187211348859318, 
        0.00322873808448831, -0.000133175875957688}}}});
auto z = apx2d.eval(0.4, 0.6); // = 0.2825.. (f(0.4, 0.6) = 0.2889..)

Below are 3D plots of the function and the approximation. They look identical.

f(x,y)

apx2d(x,y)

Painting the difference between function and approximation shows the error surface. The error is reasonable, but the two-step process makes x and y axes visible in the error. Increasing the number of coefficients would improve precision significantly - but remember that all coefficients are used in each z evaluation. A 5x5 approximation already involves 25 coefficients.

err ~0.01: f(x,y) - apx2d(x,y)

This 2D variant is not a mathematically ideal solution - better approaches exist for approximating 2D surfaces. However, this code is surprisingly simple and works well in situations where smoothness is more important than precision. I personally used this to approximate correction coefficients for arbitary-slope filters.