This is a digital

resonant filter, similar to that used in

analog synthesizers such as

TB-303.
Disclaimer: I am not an

electrical engineer; I just looked up an algorithm in a textbook on computer music and then implemented it in C.

#include <math.h>
/* The equation for an IIR resonant filter:
S = sampling rate in Hz = OUT_FREQ
c = center frequency of filter
B = bandwidth of filter
x(n) = nth sample of signal
Calculate coefficients:
R = e^(-Pi*B/S) (approximation for pole radius)
G = 1 - R
b1 = 2*R*cos(2*Pi*c/S)
b2 = -R^2
y(n) = G*(x(n) - R*x(n - 2)) + b1*y(n - 1) + b2*y(n - 2)
Algorithm after page 269 of
Moore, F. Richard. _Elements_of_Computer_Music_.
Upper Saddle River: PTR Prentice Hall, 1990.
*/
typedef struct REZFILTER
{
float R, b1, b2;
float xnm1, xnm2, ynm1, ynm2;
} REZFILTER;
/* resonate() **************************
* Compute the next sample of a resonant filter.
*/
float resonate(REZFILTER *r, float xn)
{
float yn;
yn = (1 - r->R)*(xn - r->R*r->xnm2) +
r->b1*r->ynm1 + r->b2*r->ynm2;
r->xnm2 = r->xnm1;
r->xnm1 = xn;
r->ynm2 = r->ynm1;
r->ynm1 = yn;
return yn;
}
/* set_resonate() **********************
* Set a resonant filter's center frequency and bandwidth,
* both in fractions of sample rate. For instance, 1000 Hz
* in a 44100 Hz sample is represented as (1000./44100)
*/
void set_resonate(REZFILTER *r,
float center, float bandwidth)
{
float R = exp(-M_PI * bandwidth);
r->R = R;
r->b1 = 2*R*cos(2 * M_PI * center);
r->b2 = -R*R;
}

To use this filter, first initialize all state elements (xnm1, xnm2, ynm1, ynm2) to 0, as set_resonate() doesn't do this (so as to avoid pops in time-varying filters). Then call set_resonate() with your desired center frequency and bandwidth (as real fractions of the sampling rate).

set_resonate(&my_filter, 1000./OUT_FREQ, 40./OUT_FREQ);

Then call the resonate() function with each sample, in order:

yn += resonate(&my_filter, xn);