audio_spectrum.hpp

Go to the documentation of this file.

#pragma once

 
#include "audio_metrics.hpp" // kSilenceFloorDb + buffer types
 
#include <pulp/audio/buffer.hpp>
 
#include <cstdint>
#include <span>
#include <string>
#include <vector>
 
namespace pulp::test::audio {
 
// ── Window functions ──────────────────────────────────────────────────────

enum class Window {
    rectangular,
    hann,
    hamming,
    blackman,
    flat_top,
    kaiser,
};

inline constexpr double kDefaultKaiserBeta = 14.0;

std::string window_name(Window w);

double window_coherent_gain(Window w, int n,
                            double kaiser_beta = kDefaultKaiserBeta);

std::vector<double> window_coefficients(Window w, int n,
                                        double kaiser_beta = kDefaultKaiserBeta);
 
// ── Magnitude / frequency response ────────────────────────────────────────

struct ResponsePoint {
    double hz = 0.0;
    double magnitude_db = 0.0;
};

struct ResponseCurve {
    std::string analyzer = "magnitude_response";
    std::string stimulus;   
    Window window = Window::rectangular;
    int fft_length = 0;
    int analysis_offset = 0;
    double sample_rate = 0.0;
    double bin_hz = 0.0;    
    std::vector<ResponsePoint> full;        
    std::vector<ResponsePoint> checkpoints; 

    double magnitude_db_at(double hz) const;
    double attenuation_db_at(double hz) const { return -magnitude_db_at(hz); }
};

struct ResponseOptions {
    int fft_length = 16384;     
    int analysis_offset = 0;    
    Window window = Window::rectangular;
    int channel = 0;            
    double kaiser_beta = kDefaultKaiserBeta;
};

ResponseCurve response_relative_to_input(
    const pulp::audio::BufferView<const float>& input,
    const pulp::audio::BufferView<const float>& output, double sample_rate,
    std::span<const double> checkpoints_hz, const ResponseOptions& options = {});

ResponseCurve magnitude_spectrum_curve(
    const pulp::audio::BufferView<const float>& signal, double sample_rate,
    std::span<const double> checkpoints_hz, const ResponseOptions& options = {});
 
// ── Phase / group delay ───────────────────────────────────────────────────

struct PhasePoint {
    double hz = 0.0;
    double magnitude_db_rel_peak = 0.0;
    double phase_rad = 0.0;
    double group_delay_samples = 0.0;
    bool defined = false;
    bool phase_defined = false;
};

struct PhaseCurve {
    std::string analyzer = "group_delay";
    std::string stimulus;
    Window window = Window::rectangular;
    int fft_length = 0;
    int analysis_offset = 0;
    double sample_rate = 0.0;
    double bin_hz = 0.0;
    double magnitude_floor_db = 0.0;
    std::vector<PhasePoint> full;        
    std::vector<PhasePoint> checkpoints; 

    bool defined_at(double hz) const;
    bool phase_defined_at(double hz) const;
    double group_delay_samples_at(double hz) const;
    double group_delay_seconds_at(double hz) const;
    double phase_radians_at(double hz) const;
    double magnitude_db_rel_peak_at(double hz) const;
};

struct GroupDelayOptions {
    int fft_length = 16384;
    int analysis_offset = 0;
    Window window = Window::rectangular;
    int channel = 0;
    double magnitude_floor_db = -60.0;
};

PhaseCurve measure_group_delay(
    const pulp::audio::BufferView<const float>& input,
    const pulp::audio::BufferView<const float>& output, double sample_rate,
    std::span<const double> checkpoints_hz,
    const GroupDelayOptions& options = {});
 
// ── THD / THD+N ───────────────────────────────────────────────────────────

struct Harmonic {
    int index = 0;          
    double hz = 0.0;
    double magnitude = 0.0;
    double db_below_fundamental = 0.0; 
};

struct ThdResult {
    std::string analyzer = "thd";
    double fundamental_hz = 0.0;
    Window window = Window::rectangular;
    int fft_length = 0;
    int analysis_offset = 0;
    double sample_rate = 0.0;
    double bin_hz = 0.0;
    bool coherent = false;    
    int num_harmonics = 0;    
    double thd = 0.0;         
    double thd_plus_n = 0.0;  
    std::vector<Harmonic> harmonics; 
 
    double thd_db() const;        
    double thd_plus_n_db() const; 
    double thd_percent() const { return thd * 100.0; }
};

struct ThdOptions {
    int fft_length = 16384;
    int analysis_offset = 0;
    int num_harmonics = 8;  
    int channel = 0;
    float amplitude = 0.5f;
    Window non_coherent_window = Window::hann;
    double kaiser_beta = kDefaultKaiserBeta;
};

ThdResult measure_thd(const pulp::audio::BufferView<const float>& signal,
                      double fundamental_hz, double sample_rate,
                      const ThdOptions& options = {});
 
// ── Tone projection ───────────────────────────────────────────────────────
//
// Everything below this line is FFT-free. A windowed FFT answers "what is in
// this signal?" and pays for the answer with leakage: a 0 dB tone smears a
// skirt across its neighbors that no window removes, only reshapes (see the
// `Window` table above — even Kaiser at β = 14 only pushes it to −124 dB).
// Least-squares tone projection answers the narrower question "how much of
// THIS exact frequency is present?" and has no skirt at all, because it never
// truncates anything into a window: it fits the basis over the whole segment
// and subtracts. That is what makes a −100 dBc claim measurable with margin.
//
// Use the FFT analyzers to DISCOVER content at unknown frequencies. Use these
// when the frequencies of interest are known in advance — which, for an
// oscillator's harmonics and their fold-back images, they are.

struct ToneFit {
    double cycles_per_sample = 0.0; 
    double sin_gain = 0.0;          
    double cos_gain = 0.0;          
    double amplitude = 0.0;         
    double fitted_energy = 0.0;     
    double residual_energy = 0.0;   
};

ToneFit fit_tone(std::span<const double> samples, double cycles_per_sample);

double tone_residual_db(std::span<const double> samples,
                        double cycles_per_sample);

double tone_gain_db(std::span<const double> samples, double cycles_per_sample,
                    double input_amplitude);
 
// ── Alias / image analysis ────────────────────────────────────────────────

double fold_frequency(double frequency_hz, double sample_rate);

enum class ComponentClass {
    harmonic,
    alias,
};

struct AliasComponent {
    int index = 0;              
    ComponentClass component_class = ComponentClass::harmonic;
    double source_hz = 0.0;     
    double hz = 0.0;            
    double amplitude = 0.0;
    double db_below_fundamental = 0.0; 
    bool in_band = true;
    bool resolved = true;
};

struct AliasOptions {
    double max_alias_frequency_hz = 20000.0;

    int num_harmonics = 64;
 
    int analysis_offset = 0;  
    int analysis_length = 0;
    int channel = 0;

    double min_separation_bins = 2.0;
};

struct AliasReport {
    std::string analyzer = "alias";
    double fundamental_hz = 0.0;
    double sample_rate = 0.0;
    double nyquist_hz = 0.0;
    double max_alias_frequency_hz = 0.0; 
    int num_harmonics = 0;
    int analysis_offset = 0;
    int analysis_length = 0;   
    double bin_hz = 0.0;       

    double fundamental_amplitude = 0.0;
    std::vector<AliasComponent> components;

    double worst_alias_db = kSilenceFloorDb;
    double worst_alias_hz = 0.0;      
    int worst_alias_index = 0;        

    double full_band_worst_alias_db = kSilenceFloorDb;
    double full_band_worst_alias_hz = 0.0;

    double noise_db = kSilenceFloorDb;

    double detection_floor_db = kSilenceFloorDb;

    bool has_unresolved_in_band_alias = false;
};

AliasReport measure_aliasing(const pulp::audio::BufferView<const float>& signal,
                             double fundamental_hz, double sample_rate,
                             const AliasOptions& options = {});

struct ModulationSidebandOptions {
    int max_order = 8;
    double tolerance_bins = 8.0;
    double band_low_hz = 20.0;
    double band_high_hz = 20000.0;
    Window window = Window::kaiser;
    double kaiser_beta = kDefaultKaiserBeta;
    int fft_length = 65536;
    int analysis_offset = 0; 
    int channel = 0;
};

struct ModulationSidebandReport {
    double carrier_hz = 0.0;
    double modulation_hz = 0.0;
    double carrier_power = 0.0;
    double sideband_db = kSilenceFloorDb;
    double residual_db = kSilenceFloorDb;
    int fft_length = 0;
    double bin_hz = 0.0;
};

ModulationSidebandReport
measure_modulation_sidebands(const pulp::audio::BufferView<const float>& signal, double carrier_hz,
                             double modulation_hz, double sample_rate,
                             const ModulationSidebandOptions& options = {});
 
} // namespace pulp::test::audio