brainwagon

Today, I decided to code up a little test of something that I’ve been meaning to do for quite some time: to write an implementation of a tunable bandpass filter. The algorithm is actually pretty darned simple, at least if you have a couple of pieces of code available. Let’s say you want to create a 500hz wide filter. You begin by designing a low pass filter with a cutoff at 250hz or so. If you don’t know how to do that, you can start by using this excellent webpage, which allows you to enter a filter specification and which dumps out code for the filtering operation. Once you have this, the algorithm is pretty simple: you generate two signals (which we will call I and Q) by multiplying your real-valued input signal by a sine and cosine at the frequency you wish to be the center of the bandpass. You then low pass filter each I and Q signals, and then remultiply them again by the sine and cosine as before. You then add the two channels back together, and voila! You’re done. Here’s the spectrogram of an original signal:

It works! Here’s the basic code, compressed down. It uses the libsndfile library to read and write .wav files.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sndfile.h>
SNDFILE*sf_inp,*sf_out;SF_INFO sf_info_inp;SF_INFO sf_info_out;
#define NZEROS 6
#define NPOLES 6
#define GAIN 1.602073263e+06
double cx[NZEROS+1],sx[NZEROS+1];double cy[NPOLES+1],sy[NPOLES+1];static double
filter(double xv[NZEROS+1],double yv[NPOLES+1],double n){xv[0]=xv[1];xv[1]=xv[2
];xv[2]=xv[3];xv[3]=xv[4];xv[4]=xv[5];xv[5]=xv[6];xv[6]=n/GAIN;yv[0]=yv[1];yv[1
]=yv[2];yv[2]=yv[3];yv[3]=yv[4];yv[4]=yv[5];yv[5]=yv[6];yv[6]=(xv[0]+xv[6])+6*(
xv[1]+xv[5])+15*(xv[2]+xv[4])+20*xv[3]+(-0.4675793776*yv[0])+(3.1600883661*yv[1
])+(-8.9237117766*yv[2])+(13.4798197900*yv[3])+(-11.4902571630*yv[4])+(
5.2416002121*yv[5]);return yv[6];}
#define BUFFERSIZE (8192)
main(){int i;double din[BUFFERSIZE];double dout[BUFFERSIZE];sf_count_t total,n;
double ang=0.0,dang;double s,c;sf_inp=sf_open("input.wav",SFM_READ,&sf_info_inp
);sf_info_out.channels=1;sf_info_out.samplerate=sf_info_inp.samplerate;
sf_info_out.format=SF_FORMAT_WAV|SF_FORMAT_PCM_16;sf_out=sf_open("output.wav",
SFM_WRITE,&sf_info_out);dang=2.0*M_PI*1000.0/sf_info_inp.samplerate;total=0;
while(total<sf_info_inp.frames){n=sf_read_double(sf_inp,din,BUFFERSIZE);for(i=0
;i<n;i++){s=sin(ang);c=cos(ang);dout[i]=2.0*(c*filter(cx,cy,c*din[i])+s*filter(
sx,sy,s*din[i]));ang+=dang;ang=fmod(ang,2.0*M_PI);}sf_write_double(sf_out,dout,
n);total+=n;}sf_close(sf_inp);sf_close(sf_out);}