#include "dsp/minblep.hpp"
#include "dsp/fft.hpp"
#include "dsp/window.hpp"


namespace rack {
namespace dsp {


void minBlepImpulse(int z, int o, float *output) {
	// Symmetric sinc array with `z` zero-crossings on each side
	int n = 2 * z * o;
	float *x = alignedNew<float>(n);
	for (int i = 0; i < n; i++) {
		float p = math::rescale((float) i, 0.f, (float) (n - 1), (float) -z, (float) z);
		x[i] = sinc(p);
	}

	// Apply window
	blackmanHarrisWindow(x, n);

	// Real cepstrum
	float *fx = alignedNew<float>(2*n);
	RealFFT rfft(n);
	rfft.rfft(x, fx);
	// fx = log(abs(fx))
	fx[0] = std::log(std::abs(fx[0]));
	for (int i = 1; i < n; i++) {
		fx[2*i] = std::log(std::hypot(fx[2*i], fx[2*i+1]));
		fx[2*i+1] = 0.f;
	}
	fx[1] = std::log(std::abs(fx[1]));
	// Clamp values in case we have -inf
	for (int i = 0; i < 2*n; i++) {
		fx[i] = std::max(-30.f, fx[i]);
	}
	rfft.irfft(fx, x);
	rfft.scale(x);

	// Minimum-phase reconstruction
	for (int i = 1; i < n / 2; i++) {
		x[i] *= 2.f;
	}
	for (int i = (n + 1) / 2; i < n; i++) {
		x[i] = 0.f;
	}
	rfft.rfft(x, fx);
	// fx = exp(fx)
	fx[0] = std::exp(fx[0]);
	for (int i = 1; i < n; i++) {
		float re = std::exp(fx[2*i]);
		float im = fx[2*i+1];
		fx[2*i] = re * std::cos(im);
		fx[2*i+1] = re * std::sin(im);
	}
	fx[1] = std::exp(fx[1]);
	rfft.irfft(fx, x);
	rfft.scale(x);

	// Integrate
	float total = 0.f;
	for (int i = 0; i < n; i++) {
		total += x[i];
		x[i] = total;
	}

	// Normalize
	float norm = 1.f / x[n - 1];
	for (int i = 0; i < n; i++) {
		x[i] *= norm;
	}

	std::memcpy(output, x, n * sizeof(float));

	// Cleanup
	alignedDelete(x);
	alignedDelete(fx);
}


} // namespace dsp
} // namespace rack