I don't have any "notes", but I do have this test code.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
double fc = 820; // max is 1/6th of sample rate
double q = 5; // allow range 0.7 to 5
double hz = 870;
double amplitude = 0.25;
double samplerate = 48e6 / 1088.0;
int length=190;
double pi = 3.141592654;
//
http://www.musicdsp.org/showArchiveComment.php?ArchiveID=92
static int32_t imult(int32_t a, int32_t b);
int rs = 30;
int main()
{
int i;
double fmult, damp;
double sample, prevsample=0;
double low, band, high;
double lowout, bandout, highout;
int64_t ifmult, idamp;
int32_t isample, iprevsample=0;
int32_t ilow, iband, ihigh;
int32_t ilowout, ibandout, ihighout;
double elow, eband, ehigh;
double scale = 32768.0 * 4096;
fmult = 2 * sin(pi * fc / (samplerate * 2)); // fmult is 0 to 1.0, for 0 to samplerate/6
damp = 1 / q;
ifmult = fmult * 65536.0 * 16384.0;
idamp = damp * 65536.0 * 16384.0;
// filter is stable when: fmult*fmult + 2*damp < 4.0
low = 0.0;
band = 0.0;
ilow = 0;
iband = 0;
for (i=0; i<length; i++) {
sample = sin(pi*2.0 * hz * i / samplerate) * amplitude;
printf("%6.3f", sample);
// compute the filter at floating point
low = low + fmult * band;
high = (sample+prevsample)/2.0 - low - damp*band;
prevsample = sample;
band = band + fmult * high;
lowout = low;
highout = high;
bandout = band;
low = low + fmult * band;
high = sample - low - damp*band;
band = band + fmult * high;
#if 1
lowout = (lowout + low) / 2.0;
bandout = (bandout + band) / 2.0;
highout = (highout + high) / 2.0;
#else
lowout = lowout;
bandout = bandout;
highout = highout;
#endif
printf(" %6.3f, %6.3f, %6.3f", lowout, bandout, highout);
isample = (int32_t)(sample * scale);
printf("%9d", isample >> 12);
// compute the filter at integers
ilow = ilow + imult(ifmult, iband);
ihigh = ((isample+iprevsample)>>1) - ilow - imult(idamp, iband);
iprevsample = isample;
iband = iband + imult(ifmult, ihigh);
ilowout = ilow;
ihighout = ihigh;
ibandout = iband;
ilow = ilow + imult(ifmult, iband);
ihigh = isample - ilow - imult(idamp, iband);
iband = iband + imult(ifmult, ihigh);
#if 1
ilowout = (ilowout + ilow) >> 1;
ibandout = (ibandout + iband) >> 1;
ihighout = (ihighout + ihigh) >> 1;
#else
ilowout = ilow;
ihighout = ihigh;
ibandout = iband;
#endif
// how much error is there between the intgers and floats?
elow = 100.0 * (lowout - (double)ilowout / scale) / lowout;
eband = 100.0 * (bandout - (double)ibandout / scale) / bandout;
ehigh = 100.0 * (highout - (double)ihighout / scale) / highout;
printf(" %7d, %7d, %7d (%.2f%%,%.2f%%,%.2f%%)",
ilowout >> 12, ibandout >> 12, ihighout >> 12, elow, eband, ehigh);
printf("\n");
}
return 0;
}
static int32_t imult(int32_t a, int32_t b)
{
int64_t result;
int32_t out;
result = (int64_t)a * (int64_t)b;
#if 1
if (result & 0x80000000) result += (int64_t)0x0000000080000000;
out = result >> 32;
return out << 2;
#else
return result >> rs;
#endif
}