Fourier Coefficients (C++)

Series

This page demonstrates C++ implementations for computing Fourier series coefficients using numerical integration. The trapezoidal rule is used, which converges exponentially for periodic functions.

Fourier Series

For a function $f (x)$ periodic on $[0, P]$ , the Fourier series is:

f (x) \sim \frac{a _{0}}{2} + n = 1 \sum \infty [a_{n} cos (\frac{2 nπ x}{P}) + b_{n} sin (\frac{2 nπ x}{P})]

The coefficients are computed as:

a_{0} = \frac{2}{P} \int_{0}^{P} f (x) d x

a_{n} = \frac{2}{P} \int_{0}^{P} f (x) cos (\frac{2 nπ x}{P}) d x

b_{n} = \frac{2}{P} \int_{0}^{P} f (x) sin (\frac{2 nπ x}{P}) d x

C++ Implementation

The following code implements Fourier coefficient calculation with a functor-based design:

#include <functional>
#include <vector>
#include <numbers>
#include <cmath>
#include <iostream>

// Helper function: linspace
std::vector<double> linspace(double start, double end, size_t num) {
    std::vector<double> result(num);
    double dx = (end - start) / (num - 1);
    for (size_t i = 0; i < num; ++i) {
        result[i] = start + i * dx;
    }
    return result;
}

// Fourier basis functions
namespace sp {
    struct fourierSine {
        int n;
        double P;
        fourierSine(int n_, double P_) : n(n_), P(P_) {}
        double operator()(double x) {
            double pi = std::numbers::pi_v<double>; 
            return std::sin(2.0*n * pi * x / P);
        }
    };

    struct fourierCosine {
        int n;
        double P;
        fourierCosine(int n_, double P_) : n(n_), P(P_) {}
        double operator()(double x) {
            double pi = std::numbers::pi_v<double>; 
            return std::cos(2.0*n * pi * x / P);
        }
    };
}

// Trapezoidal rule template
template<typename F>
double trap(F f, double a, double b, int N) {
    double h = (b - a) / N;
    double s = 0.5 * (f(a) + f(b));
    for (int i = 1; i < N; ++i)
        s += f(a + i*h);
    return s * h;
}

// Inner product of two functions
double innerProduct(std::function<double(double)> f, 
                    std::function<double(double)> g, 
                    double a, double b, int N) {
    auto integrand = [&f, &g](double x) {
        return f(x) * g(x);
    };
    return trap(integrand, a, b, N);
}

// Fourier coefficient calculations
double a0Fourier(std::function<double(double)> f, double P, int n, int N) {
    return 2.0/P * trap(f, 0.0, P, N);
}

double aFourier(std::function<double(double)> f, double P, int n, int N) {
    auto cos_wrapper = sp::fourierCosine(n, P);
    return 2.0/P * innerProduct(f, cos_wrapper, 0.0, P, N);
}

double bFourier(std::function<double(double)> f, double P, int n, int N) {
    auto sin_wrapper = sp::fourierSine(n, P);
    return 2.0/P * innerProduct(f, sin_wrapper, 0.0, P, N);
}

// Fourier interpolant class
struct fourierInterp {
    std::vector<double> an;
    std::vector<double> bn;
    double a0;
    double L;
    
    fourierInterp(std::vector<double> an_, 
                  std::vector<double> bn_, 
                  double a0_, 
                  double L_)
        : an(an_), bn(bn_), a0(a0_), L(L_) {}

    double operator()(double x) {
        double result = a0 / 2.0;
        double pi = std::numbers::pi_v<double>; 
        
        for (int i = 0; i < an.size(); i++) {
            int n = i + 1;
            result += an[i]*std::cos(2.0*n * pi * x / L) 
                    + bn[i]*std::sin(2.0*n * pi * x / L);
        }
        return result;
    }
};

// Example usage
int main() {
    // Define function to approximate: f(x) = x^2
    auto f = [](double x){ return x*x; };
    
    double L = 6.28;  // Period
    int num_coeffs = 20;
    int integration_points = 20000;
    
    // Compute coefficients
    double a0 = a0Fourier(f, L, 1, integration_points);
    std::vector<double> a(num_coeffs);
    std::vector<double> b(num_coeffs);
    
    for(int n = 0; n < num_coeffs; n++) {
        a[n] = aFourier(f, L, n+1, integration_points);
        b[n] = bFourier(f, L, n+1, integration_points);
    }
    
    // Create interpolant
    fourierInterp finterp(a, b, a0, L);
    
    // Evaluate on grid
    std::vector<double> xs = linspace(0, L, 1000);
    std::vector<double> ys(xs.size());
    for(int i = 0; i < xs.size(); i++) {
        ys[i] = finterp(xs[i]);
    }
    
    // Print some results
    std::cout << "a0 = " << a0 << std::endl;
    std::cout << "First few coefficients:" << std::endl;
    for(int i = 0; i < 5; i++) {
        std::cout << "a[" << i+1 << "] = " << a[i] 
                  << ", b[" << i+1 << "] = " << b[i] << std::endl;
    }
    
    return 0;
}

Key Features

Functor-based design for basis functions
Trapezoidal rule with template for flexibility
Kahan summation in interpolant for numerical stability
Generic function interface using std::function

Convergence

The trapezoidal rule converges exponentially for periodic functions, making it ideal for Fourier coefficient calculation. The error decreases as $O (e^{- c N})$ for some constant $c$ , where $N$ is the number of integration points.