/*
 * This program computes the frequency response of a Calypso DSP FIR filter
 * whose coefficients have been captured somewhere out in the wild.
 * The math has been taken from section 2.2.2 of this article:
 *
 * https://dspguru.com/dsp/faqs/fir/properties/
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define	NCOEFF	31

int coeff_int[NCOEFF];
float coeff_float[NCOEFF];
unsigned nsteps;
float freq_step, omega_step;
float freq, omega;

static void
coeff_to_float()
{
	unsigned n;

	for (n = 0; n < NCOEFF; n++)
		coeff_float[n] = coeff_int[n] / 16384.0f;
}

static void
do_one_freq()
{
	unsigned n;
	float angle_arg, rsum, isum;
	float gain, db;

	angle_arg = 0;
	rsum = 0;
	isum = 0;
	for (n = 0; n < NCOEFF; n++) {
		rsum += coeff_float[n] * cosf(angle_arg);
		isum -= coeff_float[n] * sinf(angle_arg);
		angle_arg += omega;
	}
	gain = hypotf(rsum, isum);
	db = log10f(gain) * 20.0f;
	printf("%.2f\t%f\t%.2f\n", freq, gain, db);
}

main(argc, argv)
	char **argv;
{
	unsigned n;

	if (argc != 3) {
		fprintf(stderr, "usage: %s coeff-file nsteps\n", argv[0]);
		exit(1);
	}
	if (read_fir_coeff_table_int(argv[1], coeff_int) < 0)
		exit(1);	/* error msg already printed */
	coeff_to_float();
	nsteps = strtoul(argv[2], 0, 0);
	if (nsteps < 1) {
		fprintf(stderr, "error: nsteps argument must be >= 1\n");
		exit(1);
	}
	freq_step = 4000.0f / nsteps;
	omega_step = M_PI / nsteps;
	n = 0;
	freq = 0;
	omega = 0;
	for (;;) {
		do_one_freq();
		if (n >= nsteps)
			break;
		n++;
		freq += freq_step;
		omega += omega_step;
	}
	exit(0);
}