
// Filter coefficients: https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html

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

#pragma warning(disable : 4996) // unsafe fopen

typedef signed char s8;
typedef unsigned char u8;
typedef signed short s16;
typedef unsigned short u16;
//typedef signed long int s32;		// on 64-bit system use "signed int"
//typedef unsigned long int u32;		// on 64-bit system use "unsigned int"
typedef signed int s32;
typedef unsigned int u32;

typedef unsigned int BOOL;
#define TRUE  1
#define FALSE 0

// setup
#define N_BANDS		16			// number of bands
#define FMIN		16.0		// minimal frequency
#define FMAX		16000.0		// maximal frequency
#define FS			44100.0		// sample frequency
#define Q			10			// filter quality

// middle frequency of the bands (rounded to integer)
int Freq[N_BANDS] = {
	16, 25, 40, 64, 100, 160, 250, 400,
	640, 1000, 1600, 2500, 4000, 6400, 10000, 16000 };

// round double to int
int MyRound(double x)
{
	if (x >= 0)
		return (int)(x + 0.5);
	else
		return (int)(x - 0.5);
}

// convert number to Q26
s32 DblToQ26(double x)
{
	x *= (1<<26);
	s32 n = MyRound(x);
	if (n > (1<<27)-1) n = (1<<27)-1;
	if (n < -(1<<27)+1) n = -(1<<27)+1;
	return n;
}

int main(int argc, char* argv[])
{
	int i;
	double fc;

	// calculate middle frequencies
	printf("Frequencies:\n");
	for (i = 0; i < N_BANDS; i++)
	{
		// fc = fmin * pow(fmax/fmin, i/(N-1))
		// Freq[i] = MyRound(FMIN*pow(FMAX/FMIN, (double)i/(N_BANDS-1)));
		printf("\t%d,", Freq[i]);
	}

	// calculate filter coefficients
	printf("\n\nCoefficients:\n");
	for (i = 0; i < N_BANDS; i++)
	{
		// middle frequency
		fc = Freq[i];

		// audio band-pass biquad: omega = 2 * pi * fc / fs
		double w0 = 2.0 * 3.14159265358979 * fc / FS;
		double cw = cos(w0);
		double sw = sin(w0);

		// alpha = sin(omega) / (2*Q)
		double alpha = sw / (2.0 * Q);

		// non-normalized coefficients (BPF band pass filter)
		double b0 = alpha;
		//double b1 = 0;
		//double b2 = -alpha;
		double a0 = 1 + alpha;
		double a1 = -2 * cw;
		double a2 = 1 - alpha;

		// normalize
		b0 /= a0;
		//b1 /= a0;
		//b2 /= a0;
		a1 /= a0;
		a2 /= a0;

		// convert to Q26
		s32 B0 = DblToQ26(b0);
		//s32 B1 = DblToQ26(b1);
		//s32 B2 = DblToQ26(b2);
		s32 A1 = DblToQ26(a1);
		s32 A2 = DblToQ26(a2);

		// export coefficients (b1 = 0, b2 = -b0)
		//   biguad structure (entries are in Q14 format, -2..+2 = -32767..+32767)
		//   struct BiquadQ14 { s16 b0, a1, a2; ... };
		//    y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]
		printf("\t{ %d, %d, %d, },\n", B0, /* B1, B2,*/ A1, A2 );
	}

	return 0;
}
