FFT Spectrum Analyzer (cpp, written by Codex)

Codex wrote this cpp project from the task below. It installed and ran on a clean Ubuntu 22.04 machine as written.

Task given to the agent:

TASK: FFT Spectrum Analyzer

Write a program that performs Fast Fourier Transform (FFT) analysis on time-domain signal data, identifying dominant frequencies, computing power spectral density, and supporting windowing functions.

FUNCTIONAL REQUIREMENTS:
- Accept a CSV file path as a command-line argument containing time-domain signal data (columns: time, amplitude)
- Compute the FFT of the signal and extract the frequency spectrum (magnitude and phase)
- Auto-detect the sampling rate from the time column, or accept it via --sample-rate flag
- Identify dominant frequencies: find the top N peaks in the magnitude spectrum (--peaks flag, default: 5) and report their frequencies, magnitudes, and phases
- Compute the Power Spectral Density (PSD) using Welch's method with configurable segment length via --segment flag
- Support windowing functions selectable via --window flag: rectangular (none), Hamming, Hanning, Blackman, and Kaiser (with configurable beta via --beta flag)
- Support inverse FFT via --inverse flag: reconstruct the time-domain signal from frequency-domain data
- Support frequency filtering: apply low-pass, high-pass, or band-pass filters via --filter flag (e.g., --filter low:1000 for 1kHz low-pass) and output the filtered signal
- Export the frequency spectrum data as CSV via --export flag
- Print analysis summary to console: sampling rate, number of samples, frequency resolution, dominant frequencies with magnitudes, and total signal power
- Save the full analysis as JSON with --output flag (default: fft_analysis.json)
- If no input is given, generate a sample signal composed of three sine waves at known frequencies (100Hz, 250Hz, 800Hz) with added white noise, sampled at 4000Hz for 1 second, analyze it, and show that the FFT correctly identifies the three component frequencies
- Handle errors: non-uniform sampling, insufficient data points, and signals with DC offset

Create a complete C++ project for a clean Ubuntu 22.04 machine with only G++ 12+ and CMake 3.22+ installed. Include:
- Source code
- CMakeLists.txt with all dependencies (direct and transitive) pinned to exact versions
- README.md with setup instructions, dependency explanations, build steps, run commands, and expected output

[]

cmake_minimum_required(VERSION 3.16)
project(fft_spectrum_analyzer LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

add_executable(fft_spectrum_analyzer src/main.cpp)

# FFT Spectrum Analyzer (C++)

## Requirements
- Ubuntu 22.04
- CMake 3.16+
- C++20 compiler

## Build
```bash
cmake -S . -B build
cmake --build build
```

## Run
```bash
./build/fft_spectrum_analyzer signal.csv --peaks 5 --window hamming
./build/fft_spectrum_analyzer signal.csv --sample-rate 4000 --segment 512 --filter low:1000 --export spectrum.csv
./build/fft_spectrum_analyzer signal.csv --inverse --output fft_analysis.json
```

If no input is given, a synthetic signal is generated (100Hz, 250Hz, 800Hz + noise).

#include <algorithm>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <map>
#include <random>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

namespace fs = std::filesystem;

struct Config {
    std::string file;
    double sampleRate = 0;
    bool hasSampleRate = false;
    int peaks = 5;
    int segment = 256;
    std::string window = "rectangular";
    double beta = 8.0;
    bool inverse = false;
    std::string filter;
    std::string exportCsv;
    std::string output = "fft_analysis.json";
};

struct Complex { double re = 0, im = 0; };

static Config parseArgs(int argc, char** argv) {
    Config cfg;
    std::vector<std::string> pos;
    for (int i = 1; i < argc; ++i) {
        std::string a = argv[i];
        if (!a.starts_with("--")) { pos.push_back(a); continue; }
        if (a == "--sample-rate" && i + 1 < argc) { cfg.sampleRate = std::stod(argv[++i]); cfg.hasSampleRate = true; }
        else if (a == "--peaks" && i + 1 < argc) cfg.peaks = std::stoi(argv[++i]);
        else if (a == "--segment" && i + 1 < argc) cfg.segment = std::stoi(argv[++i]);
        else if (a == "--window" && i + 1 < argc) cfg.window = argv[++i];
        else if (a == "--beta" && i + 1 < argc) cfg.beta = std::stod(argv[++i]);
        else if (a == "--inverse") cfg.inverse = true;
        else if (a == "--filter" && i + 1 < argc) cfg.filter = argv[++i];
        else if (a == "--export" && i + 1 < argc) cfg.exportCsv = argv[++i];
        else if (a == "--output" && i + 1 < argc) cfg.output = argv[++i];
        else throw std::runtime_error("Unknown option: " + a);
    }
    if (!pos.empty()) cfg.file = pos[0];
    return cfg;
}

static void writeSample(const fs::path& p) {
    std::ofstream out(p);
    out << "time,amplitude\n";
    int sr = 4000, N = 4000;
    std::mt19937 rng(42);
    std::uniform_real_distribution<double> u(-0.5, 0.5);
    for (int i = 0; i < N; ++i) {
        double t = i / static_cast<double>(sr);
        double x = 1.2 * std::sin(2 * M_PI * 100 * t)
                 + 0.8 * std::sin(2 * M_PI * 250 * t + 0.7)
                 + 0.5 * std::sin(2 * M_PI * 800 * t + 1.2)
                 + u(rng) * 0.35;
        out << t << "," << x << "\n";
    }
}

static void readSignal(const fs::path& p, std::vector<double>& t, std::vector<double>& x) {
    std::ifstream in(p);
    std::string line;
    std::getline(in, line);
    while (std::getline(in, line)) {
        if (line.empty()) continue;
        std::stringstream ss(line);
        std::string a, b;
        std::getline(ss, a, ',');
        std::getline(ss, b, ',');
        t.push_back(std::stod(a));
        x.push_back(std::stod(b));
    }
}

static double autoSampleRate(const std::vector<double>& t) {
    if (t.size() < 2) throw std::runtime_error("Insufficient data points");
    std::vector<double> dt;
    for (std::size_t i = 1; i < t.size(); ++i) dt.push_back(t[i] - t[i-1]);
    double mean = 0;
    for (double v : dt) mean += v;
    mean /= dt.size();
    double var = 0;
    for (double v : dt) var += (v - mean) * (v - mean);
    var /= dt.size();
    if (std::sqrt(var) > 0.05 * mean) throw std::runtime_error("Non-uniform sampling detected");
    return 1.0 / mean;
}

static double besselI0(double x) {
    double sum = 1, term = 1;
    for (int k = 1; k < 25; ++k) {
        term *= (x*x)/(4.0*k*k);
        sum += term;
    }
    return sum;
}

static std::vector<double> window(const std::string& name, int N, double beta) {
    std::vector<double> w(N, 1.0);
    for (int n = 0; n < N; ++n) {
        if (name == "hamming") w[n] = 0.54 - 0.46 * std::cos((2*M_PI*n)/(N-1));
        else if (name == "hanning") w[n] = 0.5 * (1 - std::cos((2*M_PI*n)/(N-1)));
        else if (name == "blackman") w[n] = 0.42 - 0.5*std::cos((2*M_PI*n)/(N-1)) + 0.08*std::cos((4*M_PI*n)/(N-1));
        else if (name == "kaiser") {
            double alpha = (N-1)/2.0;
            double r = (n-alpha)/alpha;
            w[n] = besselI0(beta*std::sqrt(1-r*r))/besselI0(beta);
        }
    }
    return w;
}

static std::vector<Complex> fft(const std::vector<double>& x) {
    int N = static_cast<int>(x.size());
    std::vector<Complex> X(N);
    for (int k = 0; k < N; ++k) {
        double re = 0, im = 0;
        for (int n = 0; n < N; ++n) {
            double ang = -2*M_PI*k*n/N;
            re += x[n] * std::cos(ang);
            im += x[n] * std::sin(ang);
        }
        X[k] = {re, im};
    }
    return X;
}

static std::vector<double> ifftReal(const std::vector<Complex>& X) {
    int N = static_cast<int>(X.size());
    std::vector<double> x(N, 0.0);
    for (int n = 0; n < N; ++n) {
        double re = 0;
        for (int k = 0; k < N; ++k) {
            double ang = 2*M_PI*k*n/N;
            re += X[k].re*std::cos(ang) - X[k].im*std::sin(ang);
        }
        x[n] = re / N;
    }
    return x;
}

static std::string esc(const std::string& s) {
    std::string out;
    for (char c : s) {
        if (c == '\\') out += "\\\\";
        else if (c == '"') out += "\\\"";
        else if (c == '\n') out += "\\n";
        else out += c;
    }
    return out;
}

int main(int argc, char** argv) {
    try {
        Config cfg = parseArgs(argc, argv);
        if (cfg.file.empty()) {
            cfg.file = "sample_signal.csv";
            writeSample(cfg.file);
        }

        std::vector<double> t, xraw;
        readSignal(cfg.file, t, xraw);
        double sr = cfg.hasSampleRate ? cfg.sampleRate : autoSampleRate(t);

        auto w = window(cfg.window, static_cast<int>(xraw.size()), cfg.beta);
        std::vector<double> x(xraw.size());
        for (std::size_t i = 0; i < x.size(); ++i) x[i] = xraw[i] * w[i];

        auto X = fft(x);
        int N = static_cast<int>(X.size());
        int half = N/2;

        std::vector<double> freq(half+1), mag(half+1), phase(half+1);
        for (int k = 0; k <= half; ++k) {
            freq[k] = k * sr / N;
            mag[k] = std::hypot(X[k].re, X[k].im);
            phase[k] = std::atan2(X[k].im, X[k].re);
        }

        std::vector<int> peakIdx;
        for (int i = 1; i < half; ++i) if (mag[i] > mag[i-1] && mag[i] > mag[i+1]) peakIdx.push_back(i);
        std::sort(peakIdx.begin(), peakIdx.end(), [&](int a, int b){ return mag[a] > mag[b]; });
        if ((int)peakIdx.size() > cfg.peaks) peakIdx.resize(cfg.peaks);

        // Welch PSD
        int seg = std::min(cfg.segment, N);
        int step = std::max(1, seg/2);
        int bins = seg/2 + 1;
        std::vector<double> pxx(bins, 0.0), pfreq(bins, 0.0);
        int segCount = 0;
        for (int start = 0; start + seg <= N; start += step) {
            std::vector<double> segx(seg);
            auto ww = window(cfg.window, seg, cfg.beta);
            for (int i = 0; i < seg; ++i) segx[i] = x[start+i] * ww[i];
            auto SX = fft(segx);
            for (int k = 0; k < bins; ++k) {
                pxx[k] += (SX[k].re*SX[k].re + SX[k].im*SX[k].im)/(seg*sr);
            }
            segCount++;
        }
        for (int i = 0; i < bins; ++i) {
            pxx[i] /= std::max(1, segCount);
            pfreq[i] = i * sr / seg;
        }

        std::vector<double> recon, filtered;
        if (cfg.inverse) recon = ifftReal(X);

        if (!cfg.filter.empty()) {
            auto Xf = X;
            for (int k = 0; k <= half; ++k) {
                double f = k * sr / N;
                bool keep = true;
                if (cfg.filter.rfind("low:", 0) == 0) keep = f <= std::stod(cfg.filter.substr(4));
                else if (cfg.filter.rfind("high:", 0) == 0) keep = f >= std::stod(cfg.filter.substr(5));
                else if (cfg.filter.rfind("band:", 0) == 0) {
                    auto s = cfg.filter.substr(5);
                    auto p = s.find('-');
                    double lo = std::stod(s.substr(0,p));
                    double hi = std::stod(s.substr(p+1));
                    keep = f >= lo && f <= hi;
                }
                if (!keep) {
                    Xf[k] = {0,0};
                    if (k != 0 && k != half) Xf[N-k] = {0,0};
                }
            }
            filtered = ifftReal(Xf);
        }

        if (!cfg.exportCsv.empty()) {
            std::ofstream out(cfg.exportCsv);
            out << "frequency,magnitude,phase\n";
            for (int i = 0; i <= half; ++i) out << freq[i] << "," << mag[i] << "," << phase[i] << "\n";
        }

        double power = 0;
        for (double v : xraw) power += v*v;
        power /= xraw.size();

        std::cout << std::fixed << std::setprecision(3)
                  << "Sampling rate: " << sr << " Hz\n"
                  << "Samples: " << N << "\n"
                  << "Resolution: " << (sr/N) << " Hz\n"
                  << "Top peaks:\n";
        for (std::size_t i = 0; i < peakIdx.size(); ++i) {
            int k = peakIdx[i];
            std::cout << (i+1) << ". f=" << freq[k] << "Hz mag=" << mag[k] << " phase=" << phase[k] << "\n";
        }

        std::ofstream json(cfg.output);
        json << "{\n"
             << "  \"samplingRate\": " << sr << ",\n"
             << "  \"samples\": " << N << ",\n"
             << "  \"frequencyResolution\": " << (sr/N) << ",\n"
             << "  \"totalSignalPower\": " << power << ",\n"
             << "  \"dominantFrequencies\": [\n";
        for (std::size_t i = 0; i < peakIdx.size(); ++i) {
            int k = peakIdx[i];
            json << "    {\"frequency\": " << freq[k] << ", \"magnitude\": " << mag[k] << ", \"phase\": " << phase[k] << "}";
            if (i + 1 < peakIdx.size()) json << ",";
            json << "\n";
        }
        json << "  ],\n"
             << "  \"spectrum\": {\"frequency\": [";
        for (int i = 0; i <= half; ++i) { if (i) json << ","; json << freq[i]; }
        json << "], \"magnitude\": [";
        for (int i = 0; i <= half; ++i) { if (i) json << ","; json << mag[i]; }
        json << "], \"phase\": [";
        for (int i = 0; i <= half; ++i) { if (i) json << ","; json << phase[i]; }
        json << "]},\n"
             << "  \"psd\": {\"frequency\": [";
        for (int i = 0; i < bins; ++i) { if (i) json << ","; json << pfreq[i]; }
        json << "], \"power\": [";
        for (int i = 0; i < bins; ++i) { if (i) json << ","; json << pxx[i]; }
        json << "]},\n";

        json << "  \"reconstructed\": ";
        if (recon.empty()) json << "null";
        else {
            json << "[";
            for (std::size_t i = 0; i < recon.size(); ++i) { if (i) json << ","; json << recon[i]; }
            json << "]";
        }

        json << ",\n  \"filteredSignal\": ";
        if (filtered.empty()) json << "null";
        else {
            json << "[";
            for (std::size_t i = 0; i < filtered.size(); ++i) { if (i) json << ","; json << filtered[i]; }
            json << "]";
        }
        json << "\n}\n";

    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << "\n";
        return 1;
    }
    return 0;
}