harmonic_entropy_table.hpp

class harmonic_entropy_table
{
public:
	struct config
	{
		int freqnency_resolution = 8 * 1024;
		double frequency_range = 4.0;
		int max_nd = 10000;
		double stddev = 0.02;
		double a = 1.5;
	};

	explicit harmonic_entropy_table(const config& config = config())
		: config_(config)
		, table_(config.freqnency_resolution)
	{
		std::vector<std::complex<double>> weights(config_.freqnency_resolution * 2);

		for (int n = 1; n <= config_.max_nd; ++n)
		{
			for (int d = 1; d <= config_.max_nd; ++d)
			{
				auto nd = n * d;

				if (nd > config_.max_nd)
				{
					continue;
				}

				if (std::gcd(n, d) != 1)
				{
					continue;
				}

				auto f = std::log2(static_cast<double>(n) / static_cast<double>(d)) / config_.frequency_range;
				auto indexf = f * config_.freqnency_resolution;
				auto indexf_floor = std::floor(indexf);
				auto index = static_cast<std::size_t>(indexf_floor);
				auto ratio = indexf - indexf_floor;
				auto w = 1.0 / std::sqrt(nd);

				if (index >= config_.freqnency_resolution)
				{
					continue;
				}

				weights[index] += w * (1.0 - ratio);

				if (index + 1 >= config_.freqnency_resolution)
				{
					continue;
				}

				weights[index + 1] += w * ratio;
			}
		}

		for (std::size_t i = 1; i < config_.freqnency_resolution; ++i)
		{
			weights[config_.freqnency_resolution * 2 - i] = weights[i];
		}

		std::vector<std::complex<double>> distribution(config_.freqnency_resolution * 2);

		for (int i = 0; i < config_.freqnency_resolution - 1; ++i)
		{
			auto f = static_cast<double>(i) * config_.frequency_range / (config_.freqnency_resolution);
			auto s = 1.0 / (config_.stddev * std::sqrt(boost::math::constants::two_pi<double>())) * std::exp(-f * f / (2.0 * config_.stddev * config_.stddev));

			distribution[i] = s;

			if (i != 0)
			{
				distribution[config_.freqnency_resolution * 2 - i] = s;
			}
		}

		auto weights_a = weights;
		auto distribution_a = distribution;

		for (auto& x : weights_a)
		{
			x = std::pow(x, config_.a);
		}

		for (auto& x : distribution_a)
		{
			x = std::pow(x, config_.a);
		}

		auto weights_fft = var::fft(weights);
		auto distribution_fft = var::fft(distribution);

		for (int i = 0; i < config_.freqnency_resolution * 2; ++i)
		{
			weights_fft[i] *= distribution_fft[i];
		}

		weights = var::ifft(weights_fft);

		auto weights_a_fft = var::fft(weights_a);
		auto distribution_a_fft = var::fft(distribution_a);

		for (int i = 0; i < config_.freqnency_resolution * 2; ++i)
		{
			weights_a_fft[i] *= distribution_a_fft[i];
		}

		weights_a = var::ifft(weights_a_fft);

		for (int i = 0; i < config_.freqnency_resolution; ++i)
		{
			table_[i] = 1.0 / (1.0 - config_.a) * std::log2(weights_a[i].real() / std::pow(weights[i].real(), config_.a));
		}

		auto minmax = std::minmax_element(table_.begin(), table_.end());
		auto min = *minmax.first;
		auto max = *minmax.second;

		for (auto& x : table_)
		{
			x = var::invlerp(min, max, x);
		}
	}

	double operator()(double df) const
	{
		auto indexf = std::abs(df) / config_.frequency_range * config_.freqnency_resolution;
		auto indexf_floor = std::floor(indexf);
		auto index = static_cast<std::size_t>(indexf_floor);

		if (index < config_.freqnency_resolution)
		{
			if (index == config_.freqnency_resolution - 1)
			{
				return table_[index];
			}
			else
			{
				auto ratio = indexf - indexf_floor;

				auto e1 = table_[index];
				auto e2 = table_[index + 1];
				return var::lerp(e1, e2, ratio);
			}
		}
		else
		{
			return 1.0;
		}
	};

private:
	config config_;
	std::vector<double> table_;
};