#!/usr/bin/awk -f
#
# Copyright (c) 2007 Ariff Abdullah <ariff@FreeBSD.org>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
# SUCH DAMAGE.
#

#
# FIR filter design by windowing method. This might become one of the
# funniest joke I've ever written due to too many tricks being applied to
# ensure maximum precision (well, in fact this is already have the same
# precision granularity compared to its C counterpart). Nevertheless, it's
# working, precise, dynamically tunable based on "presets" (so that
# LIARS^H^H^H^H^H audiophiles can tune it based on their own
# placebo^H^H^H^H^H^H^H needs).
#

#
# Some basic Math functions.
#
function abs(x)
{
	return (((x < 0) ? -x : x) + 0);
}

function fabs(x)
{
	return (((x < 0.0) ? -x : x) + 0.0);
}

function floor(x, r)
{
	r = int(x);
	if (r > x)
		r--;
	return (r + 0);
}

#function xfloor(x, a, r)
#{
#	split(sprintf("%f", x), a, ".");
#	r = abs(a[1]) + 0;
#	if (x < 0.0) {
#		if ((x + r) != 0.0)
#			r++;
#		r = -r;
#	}
#	return (r + 0);
#}

function pow(x, y)
{
	return (exp(1.0 * y * log(1.0 * x)));
}

#
# What the hell...
#
function shl(x, y)
{
	while (y > 0) {
		x *= 2;
		y--;
	}
	return (x);
}

function shr(x, y)
{
	while (y > 0) {
		x /= 2;
		y--;
	}
	return (x);
}

#
# Kaiser linear piecewise functions.
#
function kaiserAttn2Beta(attn, beta)
{
	if (attn < 0.0)
		return (Z_KAISER_BETA_DEFAULT);

	if (attn > 50.0)
		beta = 0.1102 * ((1.0 * attn) - 8.7);
	else if (attn > 21.0)
		beta = (0.5842 * pow((1.0 * attn) - 21.0, 0.4)) +	\
		    (0.07886 * ((1.0 * attn) - 21.0));
	else
		beta = 0.0;

	return (beta);
}

function kaiserBeta2Attn(beta, x, y, i, attn, xbeta)
{
	if (beta < Z_WINDOW_KAISER)
		return (Z_KAISER_ATTN_DEFAULT);
	
	if (beta > kaiserAttn2Beta(50.0))
		attn = ((1.0 * beta) / 0.1102) + 8.7;
	else {
		x = 21.0;
		y = 50.0;
		attn = 0.5 * (x + y);
		for (i = 0; i < 128; i++) {
			xbeta = kaiserAttn2Beta(attn)
			if (beta == xbeta ||				\
			    (i > 63 &&					\
			    fabs(beta - xbeta) < Z_KAISER_EPSILON))
				break;
			if (beta > xbeta)
				x = attn;
			else
				y = attn;
			attn = 0.5 * (x + y);
		}
	}

	return (attn);
}

function kaiserRolloff(len, attn)
{
	return (1.0 - (((1.0 * attn) - 7.95) / (((1.0 * len) - 1.0) * 14.36)));
}

#
# 0th order modified Bessel function of the first kind.
#
function I0(x, s, u, n, h, t)
{
	s = n = u = 1.0;
	h = x * 0.5;

	do {
		t = h / n;
		n += 1.0;
		t *= t;
		u *= t;
		s += u;
	} while (u >= (I0_EPSILON * s));

	return (s);
}

function wname(beta)
{
	if (beta >= Z_WINDOW_KAISER)
		return ("Kaiser");
	else if (beta == Z_WINDOW_BLACKMAN_NUTTALL)
		return ("Blackman - Nuttall");
	else if (beta == Z_WINDOW_NUTTALL)
		return ("Nuttall");
	else if (beta == Z_WINDOW_BLACKMAN_HARRIS)
		return ("Blackman - Harris");
	else if (beta == Z_WINDOW_BLACKMAN)
		return ("Blackman");
	else if (beta == Z_WINDOW_HAMMING)
		return ("Hamming");
	else if (beta == Z_WINDOW_HANN)
		return ("Hann");
	else
		return ("What The Hell !?!?");
}

function rolloff_round(x)
{
	if (x < 0.67)
		x = 0.67;
	else if (x > 1.0)
		x = 1.0;
	
	return (x);
}

function lpf(imp, n, rolloff, beta, num, i, j, x, nm, ibeta, w)
{
	rolloff = rolloff_round(rolloff + Z_NYQUIST_HOVER * (1.0 - rolloff));
	imp[0] = rolloff;

	#
	# Generate ideal sinc impulses, locate the last zero-crossing and pad
	# the remaining with 0.
	#
	# Note that there are other (faster) ways of calculating this without
	# the misery of traversing the entire sinc given the fact that the
	# distance between each zero crossings is actually the bandwidth of
	# the impulses, but it seems having 0.0001% chances of failure due to
	# limited precision.
	#
	j = n;
	for (i = 1; i < n; i++) {
		x = (M_PI * i) / (1.0 * num);
		imp[i] = sin(x * rolloff) / x;
		if (i != 1 && (imp[i] * imp[i - 1]) <= 0.0)
			j = i;
	}

	for (i = j; i < n; i++)
		imp[i] = 0.0;

	n = j;
	nm = 1.0 * (n - 1);

	if (beta >= Z_WINDOW_KAISER)
		ibeta = I0(beta);

	for (i = 1; i < n; i++) {
		if (beta >= Z_WINDOW_KAISER) {
			# Kaiser window...
			x = (1.0 * i) / nm;
			w = I0(beta * sqrt(1.0 - (x * x))) / ibeta;
		} else {
			# Cosined windows...
			x = (M_PI * i) / nm;
			if (beta == Z_WINDOW_BLACKMAN_NUTTALL) {
				# Blackman - Nuttall
				w = 0.36335819 + (0.4891775 * cos(x)) +	\
				    (0.1365995 * cos(2 * x)) +		\
				    (0.0106411 * cos(3 * x));
			} else if (beta == Z_WINDOW_NUTTALL) {
				# Nuttall
		        	w = 0.355768 + (0.487396 * cos(x)) +	\
				    (0.144232 * cos(2 * x)) +		\
				    (0.012604 * cos(3 * x));
			} else if (beta == Z_WINDOW_BLACKMAN_HARRIS) {
				# Blackman - Harris
				w = 0.422323 + (0.49755 * cos(x)) +	\
				    (0.07922 * cos(2 * x));
			} else if (beta == Z_WINDOW_BLACKMAN) {
				# Blackman
				w = 0.42 + (0.50 * cos(x)) +		\
				    (0.08 * cos(2 * x));
			} else if (beta == Z_WINDOW_HAMMING) {
				# Hamming
				w = 0.54 + (0.46 * cos(x));
			} else if (beta == Z_WINDOW_HANN) {
				# Hann
				w = 0.50 + (0.50 * cos(x));
			} else {
				# What The Hell !?!?
				w = 0.0;
			}
		}
		imp[i] *= w;
	}

	imp["impulse_length"] = n;
	imp["rolloff"] = rolloff;
}

function makeFilter(imp, nmult, rolloff, beta, num,			\
    nwing, mwing, nrolloff, i, dcgain, v, quality)
{
	nwing = ((nmult * num) / 2) + 1;

	lpf(imp, nwing, rolloff, beta, num);

	mwing = imp["impulse_length"];
	nrolloff = imp["rolloff"];
	quality = imp["quality"];

	dcgain = 0.0;
	for (i = num; i < mwing; i += num)
		dcgain += imp[i];
	dcgain *= 2.0;
	dcgain += imp[0];

	for (i = 0; i < nwing; i++) {
		imp[i] /= dcgain;
		imp["fixed", i] = floor((imp[i] * Z_COEFF_ONE) + 0.5);
	}

	if (quality > 2)
		printf("\n");
	printf("/*\n");
	printf(" *   quality = %d\n", quality);
	printf(" *    window = %s\n", wname(beta));
	if (beta >= Z_WINDOW_KAISER) {
		printf(" *             beta: %.2f\n", beta);
		printf(" *             stop: -%.2f dB\n",		\
		    kaiserBeta2Attn(beta));
	}
	printf(" *    length = %d\n", nmult);
	printf(" * bandwidth = %.2f%%", rolloff * 100.0);
	if (rolloff != nrolloff) {
		printf(" + %.2f%% = %.2f%% (nyquist hover: %.2f%%)",	\
		    (nrolloff - rolloff) * 100.0, nrolloff * 100.0,	\
		    Z_NYQUIST_HOVER * 100.0);
	}
	printf("\n");
	printf(" *     drift = %d\n", num);
	printf(" *     width = %d\n", mwing);
	printf(" */\n");
	printf("static int32_t z_coeff_q%d[%d] = {", quality, nwing);
	for (i = 0; i < nwing; i++) {
		if ((i % 7) == 0)
			printf("\n");
		v = imp["fixed", i];
		printf("  %s0x%05x,", (v < 0) ? "-" : " ", abs(v));
	}
	printf("\n};\n\n");
	printf("/*\n");
	printf(" * q%d differences.\n", quality);
	printf(" */\n");
	printf("static int32_t z_dcoeff_q%d[%d] = {", quality, nwing);
	for (i = 1; i <= nwing; i++) {
		if ((i % 7) == 1)
			printf("\n");
		v = -imp["fixed", i - 1];
		if (i != nwing)
			v += imp["fixed", i];
		printf("  %s0x%05x,", (v < 0) ? "-" : " ", abs(v));
	}
	printf("\n};\n");

	return (nwing);
}

function filter_parse(s, a, i, attn, alen)
{
	split(s, a, ":");
	alen = length(a);

	if (alen == 1 || alen == 2) {
		if (a[1] == "nyquist_hover") {
			i = 1.0 * a[2];
			Z_NYQUIST_HOVER = (i > 0.0 && i < 1.0) ? i : 0.0;
			return (-1);
		}
		i = 1;
		if (alen == 1) {
			attn = Z_KAISER_ATTN_DEFAULT;
			Popts["beta"] = Z_KAISER_BETA_DEFAULT;
		} else if (Z_WINDOWS[a[1]] < Z_WINDOW_KAISER) {
			Popts["beta"] = Z_WINDOWS[a[1]];
			i = floor(a[2]) + 1;
			i -= i % 2;
			Popts["nmult"] = i;
			if (i < 28)
				i = 28;
			i = 1.0 - (6.44 / i);
			Popts["rolloff"] = rolloff_round(i);
			return (0);
		} else {
			attn = 1.0 * a[i++];
			Popts["beta"] = kaiserAttn2Beta(attn);
		}
		i = floor(a[i]) + 1;
		i -= i % 2;
		Popts["nmult"] = i;
		if (i > 7 && i < 28)
			i = 27;
		i = kaiserRolloff(i, attn);
		Popts["rolloff"] = rolloff_round(i);

		return (0);
	}

	if (!(alen == 3 || alen == 4))
		return (-1);

	i = 2;

	if (a[1] == "kaiser") {
		if (alen > 2)
			Popts["beta"] = 1.0 * a[i++];
		else
			Popts["beta"] = Z_KAISER_BETA_DEFAULT;
	} else if (Z_WINDOWS[a[1]] < Z_WINDOW_KAISER)
		Popts["beta"] = Z_WINDOWS[a[1]];
	else if (1.0 * a[1] < Z_WINDOW_KAISER)
		return (-1);
	else
		Popts["beta"] = kaiserAttn2Beta(1.0 * a[1]);
	Popts["nmult"] = floor(a[i++]) + 1;
	Popts["nmult"] -= Popts["nmult"] % 2;
	if (a[1] == "kaiser" && alen == 3)
		i = kaiserRolloff(Popts["nmult"],			\
		    kaiserBeta2Attn(Popts["beta"]));
	else
		i = 1.0 * a[i];
	Popts["rolloff"] = rolloff_round(i);

	return (0);
}

BEGIN {
	I0_EPSILON = 1e-21;
	M_PI = atan2(0.0, -1.0);

	Z_FULL_SHIFT    = 30;
	Z_FULL_ONE      = shl(1, Z_FULL_SHIFT);

	#
	# 6, 7, or 8 depending on how much you can trade off between memory
	# consumption (due to large tables) and precision / quality.
	#
	Z_DRIFT_SHIFT   = 7;
	Z_DRIFT_ONE     = shl(1, Z_DRIFT_SHIFT);

	Z_SHIFT         = Z_FULL_SHIFT - Z_DRIFT_SHIFT;
	Z_ONE           = shl(1, Z_SHIFT);
	Z_MASK          = Z_ONE - 1;

	Z_COEFF_SHIFT   = 15 + (Z_DRIFT_SHIFT - 4);
	if (Z_COEFF_SHIFT > 19)
		Z_COEFF_SHIFT = 19;
	else if (Z_COEFF_SHIFT < 15)
		Z_COEFF_SHIFT = 15;

	Z_COEFF_ONE	= shl(1, Z_COEFF_SHIFT);
	Z_COEFF_UNSHIFT = Z_SHIFT - Z_COEFF_SHIFT;

	Z_LINEAR_SHIFT  = 8;
	Z_LINEAR_ONE    = shl(1, Z_LINEAR_SHIFT)

	Z_WINDOW_KAISER           =  0.0;
	Z_WINDOW_BLACKMAN_NUTTALL = -1.0;
	Z_WINDOW_NUTTALL          = -2.0;
	Z_WINDOW_BLACKMAN_HARRIS  = -3.0;
	Z_WINDOW_BLACKMAN         = -4.0;
	Z_WINDOW_HAMMING          = -5.0;
	Z_WINDOW_HANN             = -6.0;

	Z_WINDOWS["blackman_nuttall"] = Z_WINDOW_BLACKMAN_NUTTALL;
	Z_WINDOWS["nuttall"]          = Z_WINDOW_NUTTALL;
	Z_WINDOWS["blackman_harris"]  = Z_WINDOW_BLACKMAN_HARRIS;
	Z_WINDOWS["blackman"]         = Z_WINDOW_BLACKMAN;
	Z_WINDOWS["hamming"]          = Z_WINDOW_HAMMING;
	Z_WINDOWS["hann"]             = Z_WINDOW_HANN;

	Z_KAISER_2_BLACKMAN_BETA  = 8.568611;
	Z_KAISER_2_BLACKMAN_NUTTALL_BETA = 11.98;

	Z_KAISER_ATTN_DEFAULT = 100;
	Z_KAISER_BETA_DEFAULT = kaiserAttn2Beta(Z_KAISER_ATTN_DEFAULT);

	Z_KAISER_EPSILON = 1e-21;

	#
	# This is practically a joke.
	#
	Z_NYQUIST_HOVER = 0.0;

	if (ARGC < 2) {
		ARGC = 1;
		ARGV[ARGC++] = "100:8";
		ARGV[ARGC++] = "100:16";
		ARGV[ARGC++] = "100:32:0.7929";
		ARGV[ARGC++] = "100:64:0.8990";
		ARGV[ARGC++] = "100:128:0.9499";
	}

	printf("#ifndef _Z_COEFF_H_\n");
	printf("#define _Z_COEFF_H_\n\n");
	printf("/*\n");
	printf(" * Generated using awk/mkfilter.awk, heaven, wind and awesome.\n");
	printf(" *\n");
	printf(" * DO NOT EDIT!\n");
	printf(" */\n\n");
	imp["quality"] = 2;
	for (i = 1; i < ARGC; i++) {
		if (filter_parse(ARGV[i]) == 0) {
			beta = Popts["beta"];
			nmult = Popts["nmult"];
			rolloff = Popts["rolloff"];
			ztab[imp["quality"] - 2] =				\
			    makeFilter(imp, nmult, rolloff, beta, Z_DRIFT_ONE);
			imp["quality"]++;
		}
	}

	printf("\n");
	#
	# XXX
	#
	#if (length(ztab) > 0) {
	#	j = 0;
	#	for (i = 0; i < length(ztab); i++) {
	#		if (ztab[i] > j)
	#			j = ztab[i];
	#	}
	#	printf("static int32_t z_coeff_zero[%d] = {", j);
	#	for (i = 0; i < j; i++) {
	#		if ((i % 19) == 0)
	#			printf("\n");
	#		printf("  0,");
	#	}
	#	printf("\n};\n\n");
	#}
	#
	# XXX
	#
	printf("static const struct {\n");
	printf("\tint32_t len;\n");
	printf("\tint32_t *coeff;\n");
	printf("\tint32_t *dcoeff;\n");
	printf("} z_coeff_tab[] = {\n");
	if (length(ztab) > 0) {
		j = 0;
		for (i = 0; i < length(ztab); i++) {
			if (ztab[i] > j)
				j = ztab[i];
		}
		j = length(sprintf("%d", j));
		lfmt = sprintf("%%%dd", j);
		j = length(sprintf("z_coeff_q%d", length(ztab) + 1));
		zcfmt = sprintf("%%-%ds", j);
		zdcfmt = sprintf("%%-%ds", j + 1);

		for (i = 0; i < length(ztab); i++) {
			l = sprintf(lfmt, ztab[i]);
			zc = sprintf("z_coeff_q%d", i + 2);
			zc = sprintf(zcfmt, zc);
			zdc = sprintf("z_dcoeff_q%d", i + 2);
			zdc = sprintf(zdcfmt, zdc);
			printf("\t{ %s, %s, %s },\n", l, zc, zdc);
		}
	} else
		printf("\t{ 0, NULL, NULL }\n");
	printf("};\n\n");
	printf("#define Z_COEFF_TAB_SIZE\t\t\t\t\t\t\\\n");
	printf("\t((int32_t)(sizeof(z_coeff_tab) /");
	printf(" sizeof(z_coeff_tab[0])))\n\n");
	printf("#define Z_FULL_SHIFT\t\t%d\n", Z_FULL_SHIFT);
	printf("#define Z_FULL_ONE\t\t0x%08xU\n", Z_FULL_ONE);
	printf("\n");
	printf("#define Z_DRIFT_SHIFT\t\t%d\n", Z_DRIFT_SHIFT);
	printf("#define Z_DRIFT_ONE\t\t0x%08xU\n", Z_DRIFT_ONE);
	printf("\n");
	printf("#define Z_SHIFT\t\t\t%d\n", Z_SHIFT);
	printf("#define Z_ONE\t\t\t0x%08xU\n", Z_ONE);
	printf("#define Z_MASK\t\t\t0x%08xU\n", Z_MASK);
	printf("\n");
	printf("#define Z_COEFF_SHIFT\t\t%d\n", Z_COEFF_SHIFT);
	printf("#define Z_COEFF_UNSHIFT\t\t%d\n", Z_COEFF_UNSHIFT);
	printf("\n");
	printf("#define Z_LINEAR_SHIFT\t\t%d\n", Z_LINEAR_SHIFT);
	printf("#define Z_LINEAR_ONE\t\t0x%08xU\n", Z_LINEAR_ONE);
	printf("\n");
	printf("#define Z_QUALITY_ZOH\t\t0\n");
	printf("#define Z_QUALITY_LINEAR\t1\n");
	printf("#define Z_QUALITY_SINC\t\t%d\n", ((length(ztab) - 1) / 2) + 2);
	printf("\n");
	printf("#define Z_QUALITY_MIN\t\t0\n");
	printf("#define Z_QUALITY_MAX\t\t%d\n", length(ztab) + 1);
	printf("\n#endif /* !_Z_COEFF_H_ */\n");
}