Fourier Series

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <iomanip>
5
6
// This example computes a0, an, bn for a 2π-periodic function using uniform samples.
7
// We use f(x) = square wave: +1 for sin(x) >= 0, else -1.
8
// Coefficient formulas for 2π-periodic f:
9
//   a0 = (1/π) ∫_{-π}^{π} f(x) dx  ≈ (2/N) Σ f(x_k)
10
//   an = (1/π) ∫ f(x) cos(nx) dx   ≈ (2/N) Σ f(x_k) cos(n x_k)
11
//   bn = (1/π) ∫ f(x) sin(nx) dx   ≈ (2/N) Σ f(x_k) sin(n x_k)
12
13
static inline double PI() { return std::acos(-1.0); }
14
15
double f(double x) {
16
    return (std::sin(x) >= 0.0) ? 1.0 : -1.0; // 2π-periodic square wave
17
}
18
19
int main() {
20
    const int N = 4096;   // number of uniform samples over [0, 2π)
21
    const int M = 15;     // number of harmonics to compute
22
23
    const double twoPI = 2.0 * PI();
24
    const double h = twoPI / N; // step size
25
26
    // Precompute samples
27
    std::vector<double> samples(N);
28
    for (int k = 0; k < N; ++k) {
29
        double xk = (twoPI * k) / N; // x in [0, 2π)
30
        samples[k] = f(xk);
31
    }
32
33
    // Compute coefficients
34
    double a0 = 0.0;
35
    for (int k = 0; k < N; ++k) a0 += samples[k];
36
    a0 = (2.0 / N) * a0; // (1/π) * h * Σ = (1/π) * (2π/N) * Σ = 2/N * Σ
37
38
    std::vector<double> an(M + 1, 0.0), bn(M + 1, 0.0);
39
    for (int n = 1; n <= M; ++n) {
40
        double sumC = 0.0, sumS = 0.0;
41
        for (int k = 0; k < N; ++k) {
42
            double xk = (twoPI * k) / N;
43
            sumC += samples[k] * std::cos(n * xk);
44
            sumS += samples[k] * std::sin(n * xk);
45
        }
46
        an[n] = (2.0 / N) * sumC;
47
        bn[n] = (2.0 / N) * sumS;
48
    }
49
50
    // Print a0/2, then first few an, bn (expect only odd bn nonzero for a square wave)
51
    std::cout << std::fixed << std::setprecision(6);
52
    std::cout << "a0/2 = " << (a0 / 2.0) << "\n";
53
    for (int n = 1; n <= M; ++n) {
54
        std::cout << "n=" << n
55
                  << ": an=" << an[n]
56
                  << ", bn=" << bn[n] << "\n";
57
    }
58
    return 0;
59
}
60

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <iomanip>
5
6
static inline double PI() { return std::acos(-1.0); }
7
8
double f(double x) {
9
    // 2π-periodic square wave
10
    return (std::sin(x) >= 0.0) ? 1.0 : -1.0;
11
}
12
13
// Evaluate S_M(x) = a0/2 + sum_{n=1}^M (an cos(nx) + bn sin(nx))
14
double partial_sum(double x, double a0,
15
                   const std::vector<double>& an,
16
                   const std::vector<double>& bn) {
17
    double s = 0.5 * a0;
18
    int M = (int)an.size() - 1;
19
    for (int n = 1; n <= M; ++n) {
20
        s += an[n] * std::cos(n * x) + bn[n] * std::sin(n * x);
21
    }
22
    return s;
23
}
24
25
int main() {
26
    const int N = 4096; // samples for coefficient estimation
27
    const int M = 25;   // number of harmonics in reconstruction
28
29
    const double twoPI = 2.0 * PI();
30
31
    // 1) Estimate coefficients from uniform samples
32
    std::vector<double> samples(N);
33
    for (int k = 0; k < N; ++k) {
34
        double xk = (twoPI * k) / N;
35
        samples[k] = f(xk);
36
    }
37
38
    double a0 = 0.0;
39
    for (double v : samples) a0 += v;
40
    a0 = (2.0 / N) * a0;
41
42
    std::vector<double> an(M + 1, 0.0), bn(M + 1, 0.0);
43
    for (int n = 1; n <= M; ++n) {
44
        double sumC = 0.0, sumS = 0.0;
45
        for (int k = 0; k < N; ++k) {
46
            double xk = (twoPI * k) / N;
47
            sumC += samples[k] * std::cos(n * xk);
48
            sumS += samples[k] * std::sin(n * xk);
49
        }
50
        an[n] = (2.0 / N) * sumC;
51
        bn[n] = (2.0 / N) * sumS;
52
    }
53
54
    // 2) Evaluate the partial sum at grid points and print CSV: x, S_M(x), f(x)
55
    std::cout << std::fixed << std::setprecision(6);
56
    const int G = 200; // visualization grid
57
    for (int i = 0; i < G; ++i) {
58
        double x = (twoPI * i) / G; // [0, 2π)
59
        double s = partial_sum(x, a0, an, bn);
60
        std::cout << x << "," << s << "," << f(x) << "\n";
61
    }
62
    return 0;
63
}
64

1
#include <iostream>
2
#include <vector>
3
#include <complex>
4
#include <cmath>
5
#include <iomanip>
6
#include <algorithm>
7
8
// Radix-2 iterative FFT computing DFT with negative sign convention:
9
// Forward transform (invert=false): F[k] = sum_{n=0}^{N-1} a[n] * exp(-i 2π k n / N)
10
// Inverse (invert=true): a[n] = (1/N) sum_{k=0}^{N-1} F[k] * exp(+i 2π k n / N)
11
12
void fft(std::vector<std::complex<double>>& a, bool invert) {
13
    int n = (int)a.size();
14
15
    // bit-reversal permutation
16
    for (int i = 1, j = 0; i < n; ++i) {
17
        int bit = n >> 1;
18
        for (; j & bit; bit >>= 1) j ^= bit;
19
        j ^= bit;
20
        if (i < j) std::swap(a[i], a[j]);
21
    }
22
23
    for (int len = 2; len <= n; len <<= 1) {
24
        double ang = 2 * M_PI / len * (invert ? 1 : -1);
25
        std::complex<double> wlen(std::cos(ang), std::sin(ang));
26
        for (int i = 0; i < n; i += len) {
27
            std::complex<double> w(1.0, 0.0);
28
            for (int j = 0; j < len / 2; ++j) {
29
                std::complex<double> u = a[i + j];
30
                std::complex<double> v = a[i + j + len / 2] * w;
31
                a[i + j] = u + v;
32
                a[i + j + len / 2] = u - v;
33
                w *= wlen;
34
            }
35
        }
36
    }
37
38
    if (invert) {
39
        for (int i = 0; i < n; ++i) a[i] /= n;
40
    }
41
}
42
43
static inline double PI() { return std::acos(-1.0); }
44
45
double f(double x) { // example: 2π-periodic square wave
46
    return (std::sin(x) >= 0.0) ? 1.0 : -1.0;
47
}
48
49
int main() {
50
    int N = 4096; // must be a power of two for this FFT
51
    int M = 15;   // number of harmonics to extract
52
53
    // 1) Sample f at x_k = 2π k / N
54
    std::vector<std::complex<double>> a(N);
55
    for (int k = 0; k < N; ++k) {
56
        double xk = (2.0 * PI() * k) / N;
57
        a[k] = std::complex<double>(f(xk), 0.0);
58
    }
59
60
    // 2) Compute DFT via FFT (forward transform with negative sign)
61
    fft(a, false);
62
63
    // Now a[n] holds F[n] = sum f(x_k) e^{-i 2π n k / N}
64
    // For a 2π-periodic function, complex Fourier coefficients are c_n ≈ F[n] / N.
65
    // Convert to real form: a_n = 2 Re(c_n), b_n = -2 Im(c_n) for n ≥ 1; a_0 = 2 Re(c_0).
66
67
    std::vector<double> an(M + 1, 0.0), bn(M + 1, 0.0);
68
69
    std::complex<double> c0 = a[0] / (double)N; // n = 0
70
    double a0 = 2.0 * c0.real();
71
72
    for (int n = 1; n <= M; ++n) {
73
        std::complex<double> cn = a[n] / (double)N; // uses forward DFT definition
74
        an[n] = 2.0 * cn.real();      // cosine coefficient
75
        bn[n] = -2.0 * cn.imag();     // sine coefficient
76
    }
77
78
    // Print coefficients
79
    std::cout << std::fixed << std::setprecision(6);
80
    std::cout << "a0/2 = " << (a0 / 2.0) << "\n";
81
    for (int n = 1; n <= M; ++n) {
82
        std::cout << "n=" << n << ": an=" << an[n] << ", bn=" << bn[n] << "\n";
83
    }
84
85
    return 0;
86
}
87