// Quantitative fidelity analysis of captured audio, built to diagnose the exact
// failure modes the audio-mirror path can introduce: a wrong/guessed sample rate
// (pitch shift), periodic under-run/re-prime gaps (a "metallic" / choppy timbre),
// and dropped packets (clicks). It turns "the audio sounds slightly off" into
// numbers: pitch error in cents, SNR/THD, a click count, and a dropout count.
//
// The pitch/SNR/THD metrics assume the input is a single sine tone of a known
// frequency (use coop_tone as the source) -- that's the controlled signal that
// makes pitch error measurable. The discontinuity (click) and dropout metrics are
// content-agnostic, so they still mean something on a recording of real game audio.
//
// Header-only, no Windows / no audio-device dependency, so it is unit-tested with
// synthesized adversarial signals (tests/tone_analysis_test.cpp) and reused by the
// coop_audio_validate tool.
#pragma once

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <vector>

namespace coop
{

// WAVE_FORMAT_* tags we decode (kept local to avoid an mmreg.h dependency, matching audio_mix.hpp).
inline constexpr std::uint32_t kToneFormatPcm = 1;
inline constexpr std::uint32_t kToneFormatFloat = 3;

// One channel's worth of measured fidelity. Fields are NaN/0 when not applicable
// (e.g. pitch metrics need a known expected_hz > 0).
struct ToneReport
{
	bool valid = false;        // enough samples to analyze
	unsigned sample_rate = 0;  // the rate the samples are interpreted at (the *declared* rate)
	std::size_t frames = 0;    // mono frames analyzed
	double duration_sec = 0.0;

	// --- Level ---
	double rms = 0.0;               // 0..1
	double peak = 0.0;              // 0..1
	double clipped_fraction = 0.0;  // fraction of samples at >= 0.999 full-scale

	// --- Pitch (needs a known input tone frequency) ---
	double expected_hz = 0.0;        // the tone frequency that was played
	double dominant_hz = 0.0;        // the fundamental we recovered
	double pitch_error_ratio = 0.0;  // dominant / expected (1.0 = perfect)
	double pitch_error_cents = 0.0;  // 1200*log2(ratio); +/- ~10 cents starts to be audible

	// --- Spectral purity (tone mode) ---
	double snr_db = 0.0;        // fundamental power vs everything else (DC + harmonics excluded from "signal")
	double thd_percent = 0.0;   // harmonics 2..6 vs fundamental

	// --- Time-domain defects (content-agnostic) ---
	unsigned glitch_count = 0;       // discontinuity events (clicks): big isolated sample jumps
	double glitch_rate_per_sec = 0.0;
	unsigned dropout_count = 0;      // gaps: stretches that fall near-silent mid-signal
	double dropout_ms = 0.0;         // total duration of those gaps
};

namespace detail
{

inline constexpr double kPi = 3.14159265358979323846;

// In-place iterative radix-2 Cooley-Tukey FFT. `a.size()` must be a power of two.
inline void fft(std::vector<std::complex<double>>& a)
{
	const std::size_t n = a.size();
	for (std::size_t i = 1, j = 0; i < n; ++i)
	{
		std::size_t bit = n >> 1;
		for (; (j & bit) != 0; bit >>= 1)
		{
			j ^= bit;
		}
		j ^= bit;
		if (i < j)
		{
			std::swap(a[i], a[j]);
		}
	}
	for (std::size_t len = 2; len <= n; len <<= 1)
	{
		const double ang = -2.0 * kPi / static_cast<double>(len);
		const std::complex<double> wlen(std::cos(ang), std::sin(ang));
		for (std::size_t i = 0; i < n; i += len)
		{
			std::complex<double> w(1.0, 0.0);
			for (std::size_t k = 0; k < len / 2; ++k)
			{
				const std::complex<double> u = a[i + k];
				const std::complex<double> v = a[i + k + len / 2] * w;
				a[i + k] = u + v;
				a[i + k + len / 2] = u - v;
				w *= wlen;
			}
		}
	}
}

// Largest power of two <= n (0 for n==0).
inline std::size_t floor_pow2(std::size_t n)
{
	std::size_t p = 1;
	while ((p << 1) != 0 && (p << 1) <= n)
	{
		p <<= 1;
	}
	return n == 0 ? 0 : p;
}

} // namespace detail

// Decode one channel (default: channel 0) of interleaved PCM into normalized [-1,1]
// floats. Supports float32 and int16 (the formats the mirror's mixer handles); returns
// empty for anything else. `bytes` is the byte length of `pcm`.
inline std::vector<float> decode_channel(const std::uint8_t* pcm, std::size_t bytes, std::uint32_t format_tag,
                                         std::uint32_t bits, std::uint32_t channels, std::uint32_t channel = 0)
{
	std::vector<float> out;
	if (pcm == nullptr || channels == 0 || channel >= channels)
	{
		return out;
	}
	if (format_tag == kToneFormatFloat && bits == 32)
	{
		const std::size_t frames = bytes / (channels * 4);
		out.reserve(frames);
		const auto* f = reinterpret_cast<const float*>(pcm);
		for (std::size_t i = 0; i < frames; ++i)
		{
			out.push_back(f[i * channels + channel]);
		}
	}
	else if (format_tag == kToneFormatPcm && bits == 16)
	{
		const std::size_t frames = bytes / (channels * 2);
		out.reserve(frames);
		const auto* s = reinterpret_cast<const std::int16_t*>(pcm);
		for (std::size_t i = 0; i < frames; ++i)
		{
			out.push_back(static_cast<float>(s[i * channels + channel]) / 32768.0f);
		}
	}
	return out;
}

// Analyze a single channel of normalized float samples. `expected_hz` is the known
// input tone frequency (pass 0 to skip the pitch/SNR/THD metrics for non-tone audio;
// the click/dropout/level metrics still apply).
inline ToneReport analyze_tone(const float* samples, std::size_t frames, unsigned sample_rate,
                               double expected_hz)
{
	ToneReport r;
	r.sample_rate = sample_rate;
	r.frames = frames;
	r.expected_hz = expected_hz;
	if (samples == nullptr || frames < 64 || sample_rate == 0)
	{
		return r;
	}
	r.valid = true;
	r.duration_sec = static_cast<double>(frames) / sample_rate;

	// --- Level: RMS, peak, clipping ---
	double sumsq = 0.0;
	double peak = 0.0;
	std::size_t clipped = 0;
	for (std::size_t i = 0; i < frames; ++i)
	{
		const double x = samples[i];
		sumsq += x * x;
		const double a = std::fabs(x);
		peak = std::max(peak, a);
		if (a >= 0.999)
		{
			++clipped;
		}
	}
	r.rms = std::sqrt(sumsq / static_cast<double>(frames));
	r.peak = peak;
	r.clipped_fraction = static_cast<double>(clipped) / static_cast<double>(frames);

	// --- Discontinuity (click) detection: content-agnostic ---
	// A click is an isolated sample-to-sample jump far larger than the signal's typical
	// step. Use the median |first difference| as a robust scale (immune to the tone's own
	// slope and to a few outliers), and flag steps beyond 8x it. Group samples within a
	// short refractory window into one event so a single click isn't counted many times.
	if (frames >= 3)
	{
		std::vector<float> diff(frames - 1);
		for (std::size_t i = 1; i < frames; ++i)
		{
			diff[i - 1] = std::fabs(samples[i] - samples[i - 1]);
		}
		std::vector<float> sorted(diff);
		std::nth_element(sorted.begin(), sorted.begin() + sorted.size() / 2, sorted.end());
		const double median = sorted[sorted.size() / 2];
		const double thresh = std::max(8.0 * median, 0.02 * std::max(peak, 1e-6));
		const std::size_t refractory = std::max<std::size_t>(sample_rate / 1000, 8); // ~1 ms
		std::size_t last_event_end = 0;
		bool have_event = false;
		for (std::size_t i = 0; i < diff.size(); ++i)
		{
			if (diff[i] > thresh)
			{
				if (!have_event || i > last_event_end)
				{
					++r.glitch_count;
				}
				have_event = true;
				last_event_end = i + refractory;
			}
		}
		r.glitch_rate_per_sec = r.glitch_count / std::max(r.duration_sec, 1e-9);
	}

	// --- Dropout detection: stretches that fall near-silent in an otherwise active signal ---
	// Slide a ~5 ms window; flag windows whose RMS drops below 8% of the global RMS. Only
	// meaningful when the signal is actually present (global RMS above a small floor).
	if (r.rms > 1e-4)
	{
		const std::size_t win = std::max<std::size_t>(sample_rate * 5 / 1000, 16); // ~5 ms
		const std::size_t hop = std::max<std::size_t>(win / 2, 1);
		const double silence_thresh = 0.08 * r.rms;
		bool in_gap = false;
		std::size_t gap_first = 0;  // first silent window's start sample
		std::size_t gap_last = 0;   // last silent window's end sample
		std::size_t total_silent_samples = 0;
		auto close_gap = [&]() {
			if (in_gap)
			{
				total_silent_samples += (gap_last - gap_first);
				in_gap = false;
			}
		};
		for (std::size_t start = 0; start + win <= frames; start += hop)
		{
			double ws = 0.0;
			for (std::size_t i = 0; i < win; ++i)
			{
				const double x = samples[start + i];
				ws += x * x;
			}
			const double wr = std::sqrt(ws / static_cast<double>(win));
			if (wr < silence_thresh)
			{
				if (!in_gap)
				{
					++r.dropout_count;
					gap_first = start;
					in_gap = true;
				}
				gap_last = start + win;
			}
			else
			{
				close_gap();
			}
		}
		close_gap();
		// Gap extent (first silent window start .. last silent window end) -- more faithful
		// than counting interior windows, which drops the partially-silent edge windows.
		r.dropout_ms = static_cast<double>(total_silent_samples) * 1000.0 / sample_rate;
	}

	// --- Spectral analysis (pitch / SNR / THD), Hann-windowed FFT ---
	if (expected_hz > 0.0)
	{
		std::size_t n = detail::floor_pow2(frames);
		n = std::min<std::size_t>(n, std::size_t(1) << 18); // cap cost (~5 s @ 48k)
		if (n >= 1024)
		{
			std::vector<std::complex<double>> buf(n);
			for (std::size_t i = 0; i < n; ++i)
			{
				const double w = 0.5 - 0.5 * std::cos(2.0 * detail::kPi * i / (n - 1)); // Hann
				buf[i] = std::complex<double>(samples[i] * w, 0.0);
			}
			detail::fft(buf);
			const std::size_t half = n / 2;
			std::vector<double> mag(half);
			for (std::size_t i = 0; i < half; ++i)
			{
				mag[i] = std::abs(buf[i]);
			}
			const double bin_hz = static_cast<double>(sample_rate) / static_cast<double>(n);

			// Peak bin, ignoring DC/very low bins (skip < 20 Hz).
			std::size_t lo = std::max<std::size_t>(static_cast<std::size_t>(20.0 / bin_hz), 1);
			std::size_t peak_bin = lo;
			for (std::size_t i = lo; i < half; ++i)
			{
				if (mag[i] > mag[peak_bin])
				{
					peak_bin = i;
				}
			}
			// Quadratic (parabolic) interpolation on log-magnitude for a sub-bin estimate
			// (accurate for a Hann-windowed peak).
			double delta = 0.0;
			if (peak_bin > 0 && peak_bin + 1 < half)
			{
				const double a = std::log(mag[peak_bin - 1] + 1e-30);
				const double b = std::log(mag[peak_bin] + 1e-30);
				const double c = std::log(mag[peak_bin + 1] + 1e-30);
				const double denom = (a - 2.0 * b + c);
				if (std::fabs(denom) > 1e-30)
				{
					delta = 0.5 * (a - c) / denom;
					delta = std::max(-0.5, std::min(0.5, delta));
				}
			}
			r.dominant_hz = (static_cast<double>(peak_bin) + delta) * bin_hz;
			if (r.dominant_hz > 0.0)
			{
				r.pitch_error_ratio = r.dominant_hz / expected_hz;
				r.pitch_error_cents = 1200.0 * std::log2(r.pitch_error_ratio);
			}

			// Power in a +/-half_w bin lobe around a target frequency. A bin-centered Hann
			// peak is 3 bins wide, but an off-bin tone leaks into a wider skirt, so the
			// fundamental lobe is generous (+/-8) to contain that skirt -- otherwise the
			// leakage masquerades as noise and floors a clean tone's SNR (~35 dB instead of
			// >60). Harmonic lobes stay tight (+/-3).
			auto lobe_power = [&](double hz, long half_w) {
				const long center = static_cast<long>(std::lround(hz / bin_hz));
				double p = 0.0;
				for (long k = center - half_w; k <= center + half_w; ++k)
				{
					if (k >= 0 && static_cast<std::size_t>(k) < half)
					{
						p += mag[k] * mag[k];
					}
				}
				return p;
			};

			double total_power = 0.0;
			for (std::size_t i = lo; i < half; ++i)
			{
				total_power += mag[i] * mag[i];
			}
			const double fund_power = lobe_power(r.dominant_hz, 8);
			const double residual = std::max(total_power - fund_power, 1e-30);
			r.snr_db = 10.0 * std::log10(std::max(fund_power, 1e-30) / residual);

			double harm_power = 0.0;
			for (int h = 2; h <= 6; ++h)
			{
				const double hz = r.dominant_hz * h;
				if (hz < (sample_rate / 2.0))
				{
					harm_power += lobe_power(hz, 3);
				}
			}
			r.thd_percent = 100.0 * std::sqrt(harm_power / std::max(fund_power, 1e-30));
		}
	}

	return r;
}

} // namespace coop