C++
Library
since C++98
Intermediate

std::complex

Class template for complex number arithmetic with real and imaginary parts, arithmetic operators, math functions, and C++14 literals.

std::complexsince C++98

A class template in <complex> representing a complex number as a (real, imaginary) pair, providing arithmetic operators, equality comparison, stream I/O, and a complete set of complex-valued mathematical functions.

Overview

std::complex<T> parameterises the scalar type. The standard mandates three explicit specialisations β€” float, double, and long double β€” and since C++11 guarantees binary layout compatibility with C99's _Complex types: a complex<double> occupies exactly two consecutive double values, so reinterpreting to double[2] or passing to a C API expecting double _Complex is well-defined.

The class provides:

  • Full arithmetic: +, -, *, /, compound-assignment forms, and unary negation.
  • Equality comparison: == and !=. No ordering operators exist β€” complex numbers have no natural total order.
  • Stream I/O: operator<< writes (real,imag), operator>> reads the same format.
  • Free-function mathematics: magnitude, argument, squared norm, conjugate, projection, polar construction, and all transcendental functions overloaded for complex arguments.
  • C++14 user-defined literals for concise construction.

Integer and other non-floating-point template arguments are not required to work; implementations may reject them.

Syntax

cpp
#include <complex>

template<typename T>
class std::complex;          // primary template; T should be a floating-point type

std::complex<float>       cf;
std::complex<double>      cd;
std::complex<long double> cl;

Construction methods:

cpp
std::complex<double> a;              // default: (0, 0)
std::complex<double> b(3.0, 4.0);   // (3 + 4i)
std::complex<double> c{3.0, 4.0};   // C++11 brace-init, identical semantics

auto p = std::polar(5.0, 0.9273);   // from magnitude and phase angle (radians)

// C++14 literals β€” requires the using directive
using namespace std::complex_literals;   // C++14
auto d  = 3.0  + 4.0i;    // std::complex<double>
auto df = 3.0f + 4.0if;   // std::complex<float>
auto dl = 3.0l + 4.0il;   // std::complex<long double>

Member access:

cpp
std::complex<double> z(3.0, 4.0);
double r = z.real();   // 3.0
double i = z.imag();   // 4.0
z.real(6.0);           // setter β€” added in C++11
z.imag(8.0);           // setter β€” added in C++11

Examples

Arithmetic and standard math functions

cpp
#include <complex>
#include <iostream>
#include <cmath>

int main() {
    using C = std::complex<double>;

    C z1{3.0, 4.0};
    C z2{1.0, -2.0};

    C sum  = z1 + z2;              // (4, 2)
    C prod = z1 * z2;              // (11, -2)
    C quot = z1 / z2;              // complex division: (z1 * conj(z2)) / norm(z2)

    double mag   = std::abs(z1);   // 5.0       β€” magnitude (L2 norm)
    double phase = std::arg(z1);   // ~0.927 rad β€” phase angle, equivalent to atan2(imag, real)
    double nsq   = std::norm(z1);  // 25.0      β€” SQUARED magnitude: reΒ² + imΒ²
    C      cj    = std::conj(z1);  // (3, -4)   β€” complex conjugate
    C      pr    = std::proj(z1);  // projection onto Riemann sphere β€” C++11

    C root = std::sqrt(C{-1.0, 0.0}); // (0, 1) β€” principal square root
    C ex   = std::exp(C{0.0, M_PI});   // β‰ˆ (-1, Ξ΅) β€” Euler's identity

    std::cout << z1   << '\n';   // (3,4)
    std::cout << prod << '\n';   // (11,-2)
    std::cout << mag  << '\n';   // 5
}

Polynomial evaluation and DFT

cpp
#include <complex>
#include <vector>
#include <cmath>

// Evaluate a real-coefficient polynomial at a complex point β€” Horner's method
std::complex<double> eval_poly(const std::vector<double>& coeffs,
                               std::complex<double> z) {
    std::complex<double> result{};
    for (auto it = coeffs.rbegin(); it != coeffs.rend(); ++it)
        result = result * z + *it;
    return result;
}

// Compute a single DFT bin directly (O(N), for illustration)
std::complex<double> dft_bin(const std::vector<double>& x, std::size_t k) {
    using namespace std::complex_literals;  // C++14
    const std::size_t N = x.size();
    std::complex<double> acc{};
    for (std::size_t n = 0; n < N; ++n) {
        double angle = -2.0 * M_PI * static_cast<double>(k * n) / N;
        acc += x[n] * std::exp(1.0i * angle);
    }
    return acc;
}

Interoperating with C APIs

Since C++11 the standard guarantees that complex<T> has the same binary layout as T[2]:

cpp
#include <complex>
#include <vector>

// Hypothetical C API that takes an array of (re, im) pairs as doubles
extern "C" void fft_inplace(double* data, std::size_t n);

void run_fft(std::vector<std::complex<double>>& buf) {
    // reinterpret_cast here is explicitly permitted by the C++11 standard
    fft_inplace(reinterpret_cast<double*>(buf.data()), buf.size());
}

Best Practices

Choose double by default. std::complex<float> saves memory on large arrays, but precision loss compounds through iterative algorithms. Profile and measure before downgrading.

Use std::polar for phase-domain construction. Writing r * std::cos(theta) and r * std::sin(theta) manually is verbose and slightly less precise; the library version may exploit fused multiply-add or trigonometric co-computation internally.

Compare distances with std::norm, not std::abs. To check whether two complex values are within eps of each other, write std::norm(z1 - z2) < eps * eps. This avoids a square root and is the idiomatic pattern for proximity tests.

Scope complex literals tightly. Put using namespace std::complex_literals; inside a function, not at file scope, to avoid naming conflicts β€” particularly with user-defined or third-party _i literals.

Prefer free-function std::abs over manual hypot. std::abs(z) is numerically safe against intermediate overflow; std::sqrt(z.real()*z.real() + z.imag()*z.imag()) overflows for large components even when the true magnitude is representable.

Common Pitfalls

std::norm is squared magnitude, not the L2 norm. The name is actively misleading. std::norm(z) returns reΒ² + imΒ². To get |z|, use std::abs(z). Confusing these introduces silent numeric errors that only appear at scale.

No std::numeric_limits specialisation. std::numeric_limits<std::complex<double>> is not specialised. It compiles but returns the unspecialised-template defaults β€” typically 0 or false for every member. Query limits through the scalar type: std::numeric_limits<double>::epsilon().

Mixed-precision types do not implicitly convert. std::complex<double> and std::complex<float> are unrelated types. Mixing them in arithmetic is a compile error. Convert explicitly: static_cast<std::complex<double>>(cf).

operator>> format is strict. std::cin >> z expects (3,4) β€” parentheses, a comma, no spaces inside. A bare 3.0 4.0 will not parse correctly; the stream enters a failed state.

i suffix requires the using directive. operator""i, operator""if, and operator""il live in std::literals::complex_literals, not in the global namespace or even in std. A bare using namespace std; does not bring them in. Without using namespace std::complex_literals;, 4.0i is a syntax error.

Transcendental Functions and Branch Cuts

<complex> provides overloads of all standard transcendental functions for complex arguments. Each operates on the principal branch β€” the unique value chosen when the function is multivalued.

Exponential and logarithm

std::exp(z) computes e^z. For z = x + iy this equals e^x Β· (cos y + i sin y) β€” the formula underlying Euler's identity. The exponential function is entire (no branch cuts, defined everywhere in the complex plane).

std::log(z) computes the principal complex logarithm with a branch cut along the negative real axis. For a number in polar form rΒ·e^(iΞΈ) with ΞΈ ∈ (βˆ’Ο€, Ο€], the result is ln r + iΞΈ. Crossing the branch cut (approaching a negative real number from above vs below) gives results that differ in the sign of the imaginary part:

cpp
#include <complex>
#include <numbers>   // C++20
#include <iostream>

int main() {
    // log(-1 + 0i)  β†’ (0, +Ο€)  β€” approached from above
    std::cout << std::log(std::complex<double>(-1.0,  0.0)) << '\n'; // (0,3.14159)

    // log(-1 - 0i)  β†’ (0, -Ο€)  β€” approached from below (other side of the cut)
    std::cout << std::log(std::complex<double>(-1.0, -0.0)) << '\n'; // (0,-3.14159)

    // Recovering polar components via log:
    std::complex<double> z{3.0, 4.0};
    double r     = std::exp(std::log(z).real()); // magnitude: 5.0
    double theta = std::log(z).imag();            // phase angle: ~0.927 rad
    (void)r; (void)theta;
}

std::log10(z) follows the same branch cut convention; it equals std::log(z) / std::log(10.0).

Power and square root

std::pow(z, w) computes z^w via the identity exp(w * log(z)). Because log has a branch cut along the negative real axis, so does pow in its first argument. Integer exponents are not special-cased β€” use repeated multiplication or the overload that accepts a real exponent if precision matters.

std::sqrt(z) returns the principal square root whose real part is non-negative (the right half-plane). This means std::sqrt(-1.0 + 0i) gives (0, 1), not (0, -1).

cpp
#include <complex>
#include <cmath>
#include <iostream>

int main() {
    using C = std::complex<double>;

    // Euler's identity: e^(iΟ€) + 1 = 0
    const double pi = std::acos(-1.0);
    C euler = std::exp(C{0.0, pi});
    std::cout << euler << '\n';   // (-1,1.22465e-16) β€” imaginary part is floating-point noise

    // De Moivre's theorem: (rΒ·e^(iΞΈ))^n = r^n Β· e^(iΒ·nΒ·ΞΈ)
    C z = std::polar(2.0, pi / 6.0);  // magnitude 2, angle 30Β°
    C z3 = std::pow(z, 3.0);           // magnitude 8, angle 90Β° β†’ (0, 8)
    std::cout << z3 << '\n';           // (-2.66454e-15,8)

    // Principal square root β€” real part is always β‰₯ 0
    std::cout << std::sqrt(C{-1.0, 0.0}) << '\n';   // (0,1)
    std::cout << std::sqrt(C{ 0.0, 2.0}) << '\n';   // (1,1)

    // pow with a complex exponent
    C z_pow = std::pow(C{1.0, 1.0}, C{0.5, 0.5});
    std::cout << z_pow << '\n';   // some complex result via exp(wΒ·log(z))
}

Trigonometric and hyperbolic functions

All six trig and hyperbolic functions have complex overloads, plus C++11 added inverse forms (std::asin, std::acos, std::atan, std::asinh, std::acosh, std::atanh). They satisfy the standard identities extended to the complex domain β€” for example, std::cosh(C{0, theta}) equals std::cos(theta) for real theta.

Performance note: when the real part of the argument to std::exp is exactly zero at compile time, prefer std::polar(1.0, theta) or std::complex<double>{std::cos(theta), std::sin(theta)}. The pure-imaginary path through std::exp can be 4–5Γ— slower because most implementations do not special-case it.

Full DFT and Inverse DFT

The Discrete Fourier Transform is the canonical use case for std::complex. The forward transform maps N time-domain samples to N frequency bins; the inverse reconstructs the original signal. Both use twiddle factors of the form e^(βˆ’2Ο€iΒ·kΒ·n / N).

cpp
#include <complex>
#include <vector>
#include <cmath>
#include <numbers>   // C++20 β€” remove and replace pi if targeting C++17

// Forward DFT: time-domain real samples β†’ complex frequency bins (O(NΒ²))
std::vector<std::complex<double>>
dft(const std::vector<double>& x) {
    const std::size_t N = x.size();
    const double tau = 2.0 * std::numbers::pi;   // 2Ο€
    std::vector<std::complex<double>> X(N);

    for (std::size_t k = 0; k < N; ++k) {
        for (std::size_t n = 0; n < N; ++n) {
            double angle = -tau * static_cast<double>(k * n) / static_cast<double>(N);
            // Use std::polar for the twiddle factor: unit-magnitude complex exponential
            X[k] += x[n] * std::polar(1.0, angle);
        }
    }
    return X;
}

// Inverse DFT: complex frequency bins β†’ real time-domain samples (O(NΒ²))
std::vector<double>
idft(const std::vector<std::complex<double>>& X) {
    const std::size_t N = X.size();
    const double tau = 2.0 * std::numbers::pi;
    std::vector<std::complex<double>> x(N);

    for (std::size_t n = 0; n < N; ++n) {
        for (std::size_t k = 0; k < N; ++k) {
            double angle = tau * static_cast<double>(k * n) / static_cast<double>(N);
            x[n] += X[k] * std::polar(1.0, angle);
        }
        x[n] /= static_cast<double>(N);
    }

    // For a real-valued input the imaginary parts should be negligible after rounding
    std::vector<double> result(N);
    for (std::size_t n = 0; n < N; ++n)
        result[n] = x[n].real();
    return result;
}

// Example: DFT of {2, 3, 5, 7, 11}
// Forward: [28, -3.382+8.784i, -5.618+2.800i, -5.618-2.800i, -3.382-8.784i]
// Inverse round-trips back to {2, 3, 5, 7, 11}

Key design choices in this implementation:

  • std::polar(1.0, angle) constructs the unit-magnitude twiddle factor directly from the phase angle β€” cleaner and at least as fast as writing std::complex<double>{std::cos(angle), std::sin(angle)}.
  • The accumulation happens in std::complex<double> even for a real input, because intermediate results are genuinely complex.
  • The imaginary part of each x[n] in the inverse transform should be near-zero for a real-valued input; taking .real() discards the floating-point rounding residue.

For production code, replace this O(NΒ²) implementation with a radix-2 Cooley-Tukey FFT (O(N log N)) or a library such as FFTW or Intel oneAPI MKL, both of which accept std::complex<double> arrays via the guaranteed memory layout.

Memory Layout, SIMD, and Performance

Array-of-Structs layout and interop

The C++11 standard mandates that complex<T> stores real and imaginary parts in consecutive adjacent memory with no padding. Implementations differ in representation (MSVC uses T[2], GCC/libstdc++ wraps a C _Complex, Clang/libc++ uses two T members) but all satisfy:

cpp
std::complex<double> z{3.0, 4.0};
// Guaranteed: reinterpret_cast<double*>(&z)[0] == 3.0 (real)
//             reinterpret_cast<double*>(&z)[1] == 4.0 (imag)

// For arrays, every pair of consecutive doubles forms one complex number:
std::vector<std::complex<double>> buf(N);
double* raw = reinterpret_cast<double*>(buf.data());
// raw[2*i]   == buf[i].real()
// raw[2*i+1] == buf[i].imag()

This AoS (Array-of-Structs) interleaved layout is what BLAS, LAPACK, and FFTW expect β€” pass reinterpret_cast<double*>(buf.data()) directly without copying:

cpp
// Hypothetical BLAS zscal: scales a complex vector by a scalar
// Signature: void zscal_(int* n, double* alpha, double* x, int* incx);
extern "C" void zscal_(int*, double*, double*, int*);

void scale_vector(std::vector<std::complex<double>>& v, std::complex<double> alpha) {
    int n    = static_cast<int>(v.size());
    int incx = 1;
    // Passes the AoS buffer directly β€” no copy required
    zscal_(&n, reinterpret_cast<double*>(&alpha),
               reinterpret_cast<double*>(v.data()), &incx);
}

SIMD auto-vectorization and AoS vs SoA

The interleaved AoS layout makes auto-vectorization difficult. To multiply two AoS complex vectors element-wise with SIMD, the CPU must interleave real and imaginary lanes β€” this requires shuffles (e.g. VSHUFPD on AVX) and additional operations for the cross-products in complex multiplication. Compilers typically do not auto-vectorize element-wise complex multiplication on AoS buffers.

For compute-bound loops (e.g. a large element-wise product), a Structure-of-Arrays (SoA) split can be significantly faster:

cpp
#include <complex>
#include <vector>

// AoS β€” natural std::complex layout; poor auto-vectorization for multiply loops
struct ComplexAoS {
    std::vector<std::complex<float>> data;
};

// SoA β€” separate real/imag arrays; compilers can vectorize each loop independently
struct ComplexSoA {
    std::vector<float> re;
    std::vector<float> im;

    // Element-wise multiply: (a+bi)*(c+di) = (ac-bd) + (ad+bc)i
    // Each loop is a simple float[] operation β€” auto-vectorizes with SSE/AVX
    ComplexSoA multiply(const ComplexSoA& other) const {
        const std::size_t n = re.size();
        ComplexSoA result;
        result.re.resize(n);
        result.im.resize(n);
        for (std::size_t i = 0; i < n; ++i) {
            result.re[i] = re[i] * other.re[i] - im[i] * other.im[i];
            result.im[i] = re[i] * other.im[i] + im[i] * other.re[i];
        }
        return result;
    }
};

The trade-off: SoA loses the direct BLAS/FFTW interop guarantee and requires packing/unpacking when calling library functions. Choose AoS for library interop and moderate compute; choose SoA for inner loops where element-wise operations dominate.

User-Defined Literals for Complex Numbers

C++14 introduced suffix literals that make complex constant expressions read like mathematical notation. Three suffixes exist: i for double, if for float, and il for long double. They live in the inline namespace std::literals::complex_literals, reachable via using namespace std::complex_literals; or the broader using namespace std::literals;.

cpp
// C++14
#include <complex>

void demo_literals() {
    using namespace std::complex_literals;  // required in every scope that uses the suffixes

    auto z1 = 3.0 + 4.0i;       // std::complex<double>  {3.0, 4.0}
    auto z2 = 3.0f + 4.0if;     // std::complex<float>   {3.0f, 4.0f}
    auto z3 = 3.0l + 4.0il;     // std::complex<long double> {3.0l, 4.0l}

    // Pure-imaginary constants are expressed with just the suffix
    auto imag_only = 2.5i;       // std::complex<double>  {0.0, 2.5}

    // They compose with all arithmetic operators
    auto rotated = z1 * std::exp(1.0i * 1.5707963);  // z1 rotated by Ο€/2
    (void)z2; (void)z3; (void)imag_only; (void)rotated;
}

The i suffix does not conflict with a variable named i in the same scope. The literal suffix is only triggered when directly attached to a numeric token (e.g. 4.0i), not when i appears as an identifier. To avoid any ambiguity, confine the using directive to the smallest enclosing block β€” a function body rather than a translation-unit scope.

cpp
// C++14
#include <complex>
#include <vector>

// Building a filter kernel using literals β€” reads as natural math
std::vector<std::complex<double>> make_twiddles(std::size_t N) {
    using namespace std::complex_literals;
    const double tau = 6.283185307179586;     // 2Ο€
    std::vector<std::complex<double>> w(N);
    for (std::size_t k = 0; k < N; ++k) {
        double angle = -tau * static_cast<double>(k) / static_cast<double>(N);
        w[k] = std::exp(1.0i * angle);        // unit-magnitude twiddle factor
    }
    return w;
}

Complex Math Function Reference

<complex> overloads every standard transcendental function for complex arguments, extending the real-valued functions in <cmath>. These overloads are found via ADL when the argument type is std::complex<T>; no additional header is required beyond <complex>.

The table below summarises the key free functions and their semantics:

FunctionReturnsNotes
std::abs(z)TMagnitude |z| = √(re² + im²); avoids intermediate overflow
std::arg(z)TPhase angle ΞΈ = atan2(im, re); range (βˆ’Ο€, Ο€]
std::norm(z)TSquared magnitude reΒ² + imΒ² β€” not the L2 norm
std::conj(z)complex<T>Conjugate (re, βˆ’im); complex multiply by conj gives norm
std::proj(z)complex<T>Riemann sphere projection; C++11
std::polar(r, ΞΈ)complex<T>Construct from magnitude and phase
std::exp(z)complex<T>e^z = e^re Β· (cos(im) + iΒ·sin(im))
std::log(z)complex<T>Principal log; branch cut on negative real axis
std::log10(z)complex<T>log₁₀(z) = log(z) / log(10)
std::pow(z, w)complex<T>z^w via exp(wΒ·log(z)); inherits log's branch cut
std::sqrt(z)complex<T>Principal root; real part β‰₯ 0
std::sin/cos/tancomplex<T>Extended to complex domain
std::asin/acos/atancomplex<T>Inverse trig; C++11
std::sinh/cosh/tanhcomplex<T>Hyperbolic; cosh(iz) == cos(z) identity holds
std::asinh/acosh/atanhcomplex<T>Inverse hyperbolic; C++11
cpp
// C++11
#include <complex>
#include <cmath>
#include <cassert>

int main() {
    using C = std::complex<double>;

    C z{3.0, 4.0};

    // Magnitude and angle
    double mag   = std::abs(z);          // 5.0
    double phase = std::arg(z);          // atan2(4, 3) β‰ˆ 0.9273 rad
    double nsq   = std::norm(z);         // 25.0  (reΒ² + imΒ²)

    // Conjugate: multiply z * conj(z) to obtain the squared norm as a real number
    C cj = std::conj(z);                 // (3, -4)
    C product = z * cj;                  // (25, 0)  β€” imaginary part cancels
    assert(product.real() == std::norm(z) && product.imag() == 0.0);

    // Exponential and logarithm
    C lz = std::log(z);                  // (ln 5, atan2(4,3))
    C ez = std::exp(lz);                 // round-trips back to (3, 4) within floating-point error

    // Power: z^2 via exp(2*log(z))
    C z2 = std::pow(z, 2.0);            // (-7, 24)  [= (3+4i)Β² = 9-16 + 24i]

    // Square root: principal branch, real part β‰₯ 0
    C sr = std::sqrt(C{-4.0, 0.0});     // (0, 2)  β€” principal square root of -4

    // Trig identity: sinΒ²(z) + cosΒ²(z) == 1 holds in the complex domain
    C sinz = std::sin(z);
    C cosz = std::cos(z);
    C identity = sinz * sinz + cosz * cosz;   // β‰ˆ (1, 0) up to floating-point error

    (void)mag; (void)phase; (void)nsq; (void)lz; (void)ez; (void)z2; (void)sr; (void)identity;
}

std::norm is the correct function to use for fast proximity checks β€” computing std::norm(a - b) < eps * eps avoids the square root that std::abs computes and is therefore cheaper in tight loops. Prefer std::abs only when you need the actual magnitude.

FFT with Complex Arrays and FFTW Interop

std::vector<std::complex<double>> is the natural buffer type for FFT work because its memory layout is already the interleaved (re, im) pair format that every major FFT library expects. No data conversion or extra copy is needed to hand the buffer to FFTW, Intel MKL, or Apple vDSP.

The guaranteed layout β€” two consecutive T values per element β€” means a complex<double>* and an fftw_complex* (which is double[2]) are aliasing-compatible in practice, and the C++ standard explicitly permits the reinterpret_cast.

cpp
// C++11
#include <complex>
#include <vector>
#include <cstddef>

// fftw3.h declares:  typedef double fftw_complex[2];
// Forward-declare just enough to compile without the header installed:
extern "C" {
    typedef double fftw_complex[2];
    typedef struct fftw_plan_s* fftw_plan;
    fftw_plan fftw_plan_dft_1d(int, fftw_complex*, fftw_complex*, int, unsigned);
    void      fftw_execute(fftw_plan);
    void      fftw_destroy_plan(fftw_plan);
}

constexpr int FFTW_FORWARD  = -1;
constexpr int FFTW_ESTIMATE =  64;

// Compute the forward FFT of a complex-valued signal using FFTW
// Input and output are both std::vector<std::complex<double>>
void fftw_forward(const std::vector<std::complex<double>>& in,
                        std::vector<std::complex<double>>& out) {
    const int N = static_cast<int>(in.size());
    out.resize(N);

    // Guaranteed by C++11: &v[0] is layout-compatible with fftw_complex*
    auto* src = reinterpret_cast<fftw_complex*>(
                    const_cast<std::complex<double>*>(in.data()));
    auto* dst = reinterpret_cast<fftw_complex*>(out.data());

    fftw_plan plan = fftw_plan_dft_1d(N, src, dst, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan);
    fftw_destroy_plan(plan);
}

For real-valued input signals, FFTW's r2c transform is twice as fast and produces only N/2+1 unique output bins (the spectrum is Hermitian-symmetric). However, r2c output is still fftw_complex*, so you can map a std::vector<std::complex<double>> of size N/2+1 onto it with the same reinterpret_cast idiom.

cpp
// C++11
#include <complex>
#include <vector>
#include <algorithm>

// Build a complex signal buffer from separate real and imaginary arrays.
// Useful when converting from SoA (Structure-of-Arrays) to the AoS layout
// that FFT libraries expect.
std::vector<std::complex<double>>
pack_complex(const std::vector<double>& re, const std::vector<double>& im) {
    const std::size_t N = re.size();
    std::vector<std::complex<double>> buf(N);
    for (std::size_t i = 0; i < N; ++i)
        buf[i] = {re[i], im[i]};
    return buf;
}

// Unpack back to SoA after an FFT (e.g. for element-wise SIMD processing)
void unpack_complex(const std::vector<std::complex<double>>& buf,
                    std::vector<double>& re, std::vector<double>& im) {
    const std::size_t N = buf.size();
    re.resize(N); im.resize(N);
    for (std::size_t i = 0; i < N; ++i) {
        re[i] = buf[i].real();
        im[i] = buf[i].imag();
    }
}

// Compute the power spectrum (magnitude squared of each FFT bin) in-place
void power_spectrum(const std::vector<std::complex<double>>& X,
                    std::vector<double>& power) {
    power.resize(X.size());
    // std::norm avoids the sqrt β€” correct for power spectrum (|X[k]|Β²)
    std::transform(X.begin(), X.end(), power.begin(),
                   [](const std::complex<double>& x) { return std::norm(x); });
}

Key points for FFTW interop:

  • Allocate the buffer with std::vector<std::complex<double>> (or fftw_alloc_complex for SIMD-aligned plans with FFTW_MEASURE).
  • Create plans before executing; plans are reusable and expensive to build β€” cache them for repeated transforms of the same size.
  • FFTW_ESTIMATE plans avoid the expensive measurement pass but may be slower to execute than FFTW_MEASURE plans on repeated calls.
  • Thread safety: plan creation is not thread-safe; execution of distinct plans on distinct buffers is.

See Also

  • <cmath> β€” scalar overloads of sin, cos, exp, etc.; <complex> provides separate overloads for complex arguments
  • std::valarray β€” expression-template array type that composes naturally with complex<T> elements for vectorised numerics
  • Operator overloading β€” std::complex is a canonical example of arithmetic operator overloading for user-defined numeric types