#include "BlipBuffer.hh"

#include "Math.hh"
#include "cstd.hh"
#include "narrow.hh"
#include "ranges.hh"
#include "xrange.hh"

#include <algorithm>
#include <array>
#include <cassert>
#include <iostream>

namespace openmsx {

// Number of samples in a (pre-calculated) impulse-response wave-form.
static constexpr int BLIP_IMPULSE_WIDTH = 16;

// Derived constants
static constexpr int BLIP_RES = 1 << BlipBuffer::BLIP_PHASE_BITS;


// Precalculated impulse table.
static constexpr auto impulses = [] {
	constexpr int HALF_SIZE = BLIP_RES / 2 * (BLIP_IMPULSE_WIDTH - 1);
	std::array<double, BLIP_RES + HALF_SIZE + BLIP_RES> fImpulse = {};
	std::span<double, HALF_SIZE> out = subspan<HALF_SIZE>(fImpulse, BLIP_RES);
	std::span<double, BLIP_RES>  end = subspan<BLIP_RES >(fImpulse, BLIP_RES + HALF_SIZE);

	// generate sinc, apply hamming window
	double overSample = ((4.5 / (BLIP_IMPULSE_WIDTH - 1)) + 0.85);
	double to_angle = Math::pi / (2.0 * overSample * BLIP_RES);
	double to_fraction = Math::pi / (2 * (HALF_SIZE - 1));
	for (auto i : xrange(HALF_SIZE)) {
		double angle = ((i - HALF_SIZE) * 2 + 1) * to_angle;
		out[i] = cstd::sin<2>(angle) / angle;
		out[i] *= 0.54 - 0.46 * cstd::cos<2>((2 * i + 1) * to_fraction);
	}

	// need mirror slightly past center for calculation
	for (auto i : xrange(BLIP_RES)) {
		end[i] = out[HALF_SIZE - 1 - i];
	}

	// find rescale factor
	double total = 0.0;
	for (auto i : xrange(HALF_SIZE)) {
		total += out[i];
	}
	double rescale = 1.0 / (2.0 * total);

	// integrate, first difference, rescale, convert to float
	constexpr int IMPULSES_SIZE = BLIP_RES * (BLIP_IMPULSE_WIDTH / 2) + 1;
	std::array<float, IMPULSES_SIZE> imp = {};
	double sum = 0.0;
	double next = 0.0;
	for (auto i : xrange(IMPULSES_SIZE)) {
		imp[i] = float((next - sum) * rescale);
		sum += fImpulse[i];
		next += fImpulse[i + BLIP_RES];
	}
	// Original code would now apply a correction on each kernel so that
	// the (integer) coefficients sum up to 'kernelUnit'. I've measured
	// that after switching to float coefficients this correction is
	// roughly equal to std::numeric_limits<float>::epsilon(). So it's no
	// longer meaningful.

	// reshuffle values to a more cache friendly order
	std::array<std::array<float, BLIP_IMPULSE_WIDTH>, BLIP_RES> result = {};
	for (auto phase : xrange(BLIP_RES)) {
		const auto* imp_fwd = &imp[BLIP_RES - phase];
		const auto* imp_rev = &imp[phase];
		auto* p = result[phase].data();
		for (size_t i : xrange(BLIP_IMPULSE_WIDTH / 2)) {
			*p++ = imp_fwd[BLIP_RES * i];
		}
		for (ptrdiff_t i = BLIP_IMPULSE_WIDTH / 2 - 1; i >= 0; --i) {
			*p++ = imp_rev[BLIP_RES * i];
		}
	}
	return result;
}();


BlipBuffer::BlipBuffer()
{
	if (false) {
		for (auto i : xrange(BLIP_RES)) {
			std::cout << "\t{ " << impulses[i][0];
			for (auto j : xrange(1, BLIP_IMPULSE_WIDTH)) {
				std::cout << ", " << impulses[i][j];
			}
			std::cout << " },\n";
		}
	}

	std::ranges::fill(buffer, 0);
}

void BlipBuffer::addDelta(TimeIndex time, float delta)
{
	unsigned tmp = time.toInt() + BLIP_IMPULSE_WIDTH;
	assert(tmp < BUFFER_SIZE);
	availSamp = std::max(availSamp, narrow<ptrdiff_t>(tmp));

	auto phase = time.fractAsInt();
	auto ofst = time.toInt() + offset;
	const float* __restrict impulse = impulses[phase].data();
	if ((ofst + BLIP_IMPULSE_WIDTH) <= BUFFER_SIZE) [[likely]] {
		float* __restrict result = &buffer[ofst];
		for (auto i : xrange(BLIP_IMPULSE_WIDTH)) {
			result[i] += impulse[i] * delta;
		}
	} else {
		for (auto i : xrange(BLIP_IMPULSE_WIDTH)) {
			buffer[(ofst + i) & BUFFER_MASK] += impulse[i] * delta;
		}
	}
}

static constexpr float BASS_FACTOR = 511.0f / 512.0f;

template<size_t PITCH>
void BlipBuffer::readSamplesHelper(float* __restrict out, size_t samples)
{
	assert((offset + samples) <= BUFFER_SIZE);
	auto acc = accum;
	auto ofst = offset;
	for (auto i : xrange(samples)) {
		out[i * PITCH] = acc;
		acc *= BASS_FACTOR;
		acc += buffer[ofst];
		buffer[ofst] = 0.0f;
		++ofst;
	}
	accum = acc;
	offset = ofst & BUFFER_MASK;
}

static bool isSilent(float x)
{
	// 'accum' falls away slowly (via BASS_FACTOR). Because it's a float it
	// takes a long time before it's really zero. But much sooner we can
	// already say it's practically silent. This check is safe when the
	// input is in range [-1..+1] (does not say silent too soon). When the
	// input is in range [-32768..+32767] it takes longer before we switch
	// to silent mode (but still less than a second).
	constexpr float threshold = 1.0f / 32768;
	return std::abs(x) < threshold;
}

// TODO replace 'out + samples + PITCH' with 'stride_view' ???
template<size_t PITCH>
bool BlipBuffer::readSamples(float* __restrict out, size_t samples)
{
	if (availSamp <= 0) {
		#ifdef DEBUG
		// buffer contains all zeros (only check this in debug mode)
		assert(std::ranges::all_of(buffer, [](const auto& b) { return b == 0.0f; }));
		#endif
		if (isSilent(accum)) {
			return false; // muted
		}
		auto acc = accum;
		for (auto i : xrange(samples)) {
			out[i * PITCH] = acc;
			acc *= BASS_FACTOR;
		}
		accum = acc;
	} else {
		availSamp -= narrow<ptrdiff_t>(samples);
		auto t1 = std::min(samples, BUFFER_SIZE - offset);
		readSamplesHelper<PITCH>(out, t1);
		if (t1 < samples) {
			assert(offset == 0);
			auto t2 = samples - t1;
			assert(t2 < BUFFER_SIZE);
			readSamplesHelper<PITCH>(&out[t1 * PITCH], t2);
		}
		assert(offset < BUFFER_SIZE);
	}
	return true;
}

template bool BlipBuffer::readSamples<1>(float*, size_t);
template bool BlipBuffer::readSamples<2>(float*, size_t);

} // namespace openmsx