--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fir/freqresp.c	Mon Aug 02 04:59:46 2021 +0000
@@ -0,0 +1,81 @@
+/*
+ * 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);
+}