Windowed Sinc FIR Generator

Type : LP, HP, BP, BS
References : Posted by Bob Maling
Linked file : wsfir.h
(this linked file is included below)
Notes :
This code generates FIR coefficients for lowpass, highpass, bandpass, and bandstop filters by windowing a sinc function.

The purpose of this code is to show how windowed sinc filter coefficients are generated. Also shown is how highpass, bandpass, and bandstop filters can be made from lowpass filters.

Included windows are Blackman, Hanning, and Hamming. Other windows can be added by following the structure outlined in the opening comments of the header file.
Comments
from : Christian[AT]savioursofsoul[DOT]de
comment : // Object Pascal Port... unit SincFIR; (* Windowed Sinc FIR Generator Bob Maling (BobM.DSP@gmail.com) Contributed to musicdsp.org Source Code Archive Last Updated: April 12, 2005 Translated to Object Pascal by Christian-W. Budde Usage: Lowpass:wsfirLP(H, WindowType, CutOff) Highpass:wsfirHP(H, WindowType, CutOff) Bandpass:wsfirBP(H, WindowType, LowCutOff, HighCutOff) Bandstop:wsfirBS(H, WindowType, LowCutOff, HighCutOff) where: H (TDoubleArray): empty filter coefficient table (SetLength(H,DesiredLength)!) WindowType (TWindowType): wtBlackman, wtHanning, wtHamming CutOff (double): cutoff (0 < CutOff < 0.5, CutOff = f/fs) --> for fs=48kHz and cutoff f=12kHz, CutOff = 12k/48k => 0.25 LowCutOff (double):low cutoff (0 < CutOff < 0.5, CutOff = f/fs) HighCutOff (double):high cutoff (0 < CutOff < 0.5, CutOff = f/fs) Windows included here are Blackman, Blackman-Harris, Gaussian, Hanning and Hamming.*) interface uses Math; type TDoubleArray = array of Double; TWindowType = (wtBlackman, wtHanning, wtHamming, wtBlackmanHarris, wtGaussian); // Window type contstants // Function prototypes procedure wsfirLP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double); procedure wsfirHP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double); procedure wsfirBS(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double); procedure wsfirBP(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double); procedure genSinc(var Sinc : TDoubleArray; const CutOff : Double); procedure wGaussian(var W : TDoubleArray); procedure wBlackmanHarris(var W : TDoubleArray); procedure wBlackman(var W : TDoubleArray); procedure wHanning(var W : TDoubleArray); procedure wHamming(var W : TDoubleArray); implementation // Generate lowpass filter // This is done by generating a sinc function and then windowing it procedure wsfirLP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double); begin genSinc(H, CutOff); // 1. Generate Sinc function case WindowType of // 2. Generate Window function -> lowpass filter! wtBlackman: wBlackman(H); wtHanning: wHanning(H); wtHamming: wHamming(H); wtGaussian: wGaussian(H); wtBlackmanHarris: wBlackmanHarris(H); end; end; // Generate highpass filter // This is done by generating a lowpass filter and then spectrally inverting it procedure wsfirHP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double); var i : Integer; begin wsfirLP(H, WindowType, CutOff); // 1. Generate lowpass filter // 2. Spectrally invert (negate all samples and add 1 to center sample) lowpass filter // = delta[n-((N-1)/2)] - h[n], in the time domain for i:=0 to Length(H)-1 do H[i]:=H[i]*-1.0; H[(Length(H)-1) div 2]:=H[(Length(H)-1) div 2]+1.0; end; // Generate bandstop filter // This is done by generating a lowpass and highpass filter and adding them procedure wsfirBS(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double); var i : Integer; H2 : TDoubleArray; begin SetLength(H2,Length(H)); // 1. Generate lowpass filter at first (low) cutoff frequency wsfirLP(H, WindowType, LowCutOff); // 2. Generate highpass filter at second (high) cutoff frequency wsfirHP(H2, WindowType, HighCutOff); // 3. Add the 2 filters together for i:=0 to Length(H)-1 do H[i]:=H[i]+H2[i]; SetLength(H2,0); end; // Generate bandpass filter // This is done by generating a bandstop filter and spectrally inverting it procedure wsfirBP(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double); var i : Integer; begin wsfirBS(H, WindowType, LowCutOff, HighCutOff); // 1. Generate bandstop filter // 2. Spectrally invert (negate all samples and add 1 to center sample) lowpass filter // = delta[n-((N-1)/2)] - h[n], in the time domain for i:=0 to Length(H)-1 do H[i]:=H[i]*-1.0; H[(Length(H)-1) div 2]:=H[(Length(H)-1) div 2]+1.0; end; // Generate sinc function - used for making lowpass filter procedure genSinc(var Sinc : TDoubleArray; const Cutoff : Double); var i,j : Integer; n,k : Double; begin j:=Length(Sinc)-1; k:=1/j; // Generate sinc delayed by (N-1)/2 for i:=0 to j do if (i=j div 2) then Sinc[i]:=2.0*Cutoff else begin n:=i-j/2.0; Sinc[i]:=sin(2.0*PI*Cutoff*n)/(PI*n); end; end; // Generate window function (Gaussian) procedure wGaussian(var W : TDoubleArray); var i,j : Integer; k : Double; begin j:=Length(W)-1; k:=1/j; for i:=0 to j do W[i]:=W[i]*(exp(-5.0/(sqr(j))*(2*i-j)*(2*i-j))); end; // Generate window function (Blackman-Harris) procedure wBlackmanHarris(var W : TDoubleArray); var i,j : Integer; k : Double; begin j:=Length(W)-1; k:=1/j; for i:=0 to j do W[i]:=W[i]*(0.35875-0.48829*cos(2*PI*(i+0.5)*k)+0.14128*cos(4*PI*(i+0.5)*k)-0.01168*cos(6*PI*(i+0.5)*k)); end; // Generate window function (Blackman) procedure wBlackman(var W : TDoubleArray); var i,j : Integer; k : Double; begin j:=Length(W)-1; k:=1/j; for i:=0 to j do W[i]:=W[i]*(0.42-(0.5*cos(2*PI*i*k))+(0.08*cos(4*PI*i*k))); end; // Generate window function (Hanning) procedure wHanning(var W : TDoubleArray); var i,j : Integer; k : Double; begin j:=Length(W)-1; k:=1/j; for i:=0 to j do W[i]:=W[i]*(0.5*(1.0-cos(2*PI*i*k))); end; // Generate window function (Hamming) procedure wHamming(var W : TDoubleArray); var i,j : Integer; k : Double; begin j:=Length(W)-1; k:=1/j; for i:=0 to j do W[i]:=W[i]*(0.54-(0.46*cos(2*PI*i*k))); end; end.

from : scoofy[AT]inf[DOT]elte[DOT]hu
comment : The Hanning window is often incorrecty referred to as 'Hanning', since it was named after a guy called Julius von Hann. So it's more appropriate to call it 'Hann' window.

from : rousicoello[AT]gmail[DOT]com
comment : What do I have to do to apply this windowed sinc filter to a signal "Y" for example... what is the code for this? Imagine signal "Y" is a stereo song which means that I have Y1 from channel 1 and Y2 from channel 2.. help me please as soon as posible because at least i want to know how to aplly the filter to any signal...

from : Dave in sinc land
comment : I've seen a BASIC version of the same genSinc code. Shouldn't a 'sin(2.0*Cutoff)' be used when the divide by zero check is 0? ... if (i=j div 2) then Sinc[i]:=2.0*Cutoff else begin n:=i-j/2.0; Sinc[i]:=sin(2.0*PI*Cutoff*n)/(PI*n); end; ...

from : Dave in sinc land
comment : Scrap that, I've just FFT'd the response, and it appears to be correct as it was. Hey it just looked wrong o.k. ;)

from : craig101[DOT]com101[AT]yahoo[DOT]com
comment : Can anyone, provide a sample code on how to apply this to an array of samples?

from : amishman35[AT]cox[DOT]net
comment : Here is a FIR filter that works with up to about 80 coefficients on a Sony PSP running at 280 Megahertz. The multipliers (8191.0 and 32767.0) may need to be incremented by one, or not. You need to call 'fir_filter' with the sample to process in 'sample' and the floating point coefficents in '*coefs' and how long the entire filter is in 'len' and the returned sample is a 16 bit audio (signed short) sample. When it is first run it will normalize the filter. /* code starts here, remove spaces between the includes */ #include < math.h > #include < string.h > #include < stdio.h > #include < stdlib.h > void normalizeCoefs(signed short *firshort, signed short *outfilter, unsigned long len) { // filter gain at uniform frequency intervals double *g=NULL; double *dtemp=NULL; double theta, s, c, sac, sas; double gMax = -100.0; double sc = 10.0/log(10.0); double t = M_PI / len; long i=0; long n=len-1; double normFactor =0; g=(double*)malloc((len+1)*sizeof(double)); dtemp=(double*)malloc((len+1)*sizeof(double)); for (i=0; i<len; i++) { dtemp[i]=firshort[i]/32768.0; } for (i=0; i<len; i++) { theta = i*t; sac = 0.0; sas = 0.0; long k=0; for (k=0; k<len; k++) { c =cos(k*theta); s =sin(k*theta); sac += c*(dtemp[k]); sas += s*(dtemp[k]); } g[i] = sc*log(sac*sac + sas*sas); if(g[i]>gMax) { gMax=g[i]; } } // normalise to 0 dB maximum gain for (i=0; i<len; i++) { g[i] -= gMax; } // normalise coefficients normFactor =0;pow(10.0, -0.05*gMax); for (i=0; i<len; i++) { dtemp[i] *= normFactor; } n=len-1; for (i=0; i<len; i++) { outfilter[i]=dtemp[n--]*32767.0; } free(dtemp); free(g); } signed short generalFIR(short input, short *coefs, unsigned long len,unsigned char resetit) { static signed short *delay_line=NULL; static unsigned int first=1; static signed long reacc=0; unsigned long filtcnt=0; static int consumer=1; static int producer=0; unsigned long tx=0; if(resetit==1) { first=1; } if(first==1) { if(delay_line!=NULL) free(delay_line); producer=0; consumer=1; delay_line=(signed short*)malloc(reallen*sizeof(signed short)); first=0; } reacc=0; delay_line[producer++]=input; if(producer>len-1) { producer-=len; } int filtptr=consumer; for(tx=0;tx<len;tx++) { reacc+=(delay_line[filtptr++]*coefs[tx]); if(filtptr>len-1) { filtptr-=len; } } consumer++; if(consumer>len-1) { consumer-=len; } reacc = reacc >> 15; if(reacc<-32768) reacc= -32768; if(reacc>32767) reacc=32767; return (signed short)reacc; } signed short fir_filter(signed short sample,double *coefs,unsigned long len) { static unsigned char first=1; static signed short *newfir=NULL; if(first==1) { unsigned long i=0; if(newfir!=NULL) { free(newfir); newfir=NULL; } newfir=(signed short*)malloc((filterorder+1) * sizeof(signed short)); if(newfir==NULL) { return 0; } for(k=0;k<len;k++) { newfir[k]=coefs[k] * 8191.0; } normalizeCoefs(newfir,newfir,len); generalFIR(0,newfir,len,1); first=0; } return generalFIR(sample,newfir,len,0); }


Linked files
/*
	Windowed Sinc FIR Generator
	Bob Maling (BobM.DSP@gmail.com)
	Contributed to musicdsp.org Source Code Archive
	Last Updated: April 12, 2005

	Usage:
		Lowpass:	wsfirLP(h, N, WINDOW, fc)
		Highpass:	wsfirHP(h, N, WINDOW, fc)
		Bandpass:	wsfirBP(h, N, WINDOW, fc1, fc2)
		Bandstop:	wsfirBS(h, N, WINDOW, fc1, fc2)

	where:
		h (double[N]):	filter coefficients will be written to this array
		N (int):		number of taps
		WINDOW (int):	Window (W_BLACKMAN, W_HANNING, or W_HAMMING)
		fc (double):	cutoff (0 < fc < 0.5, fc = f/fs)
						--> for fs=48kHz and cutoff f=12kHz, fc = 12k/48k => 0.25
		
		fc1 (double):	low cutoff (0 < fc < 0.5, fc = f/fs)
		fc2 (double):	high cutoff (0 < fc < 0.5, fc = f/fs)


	Windows included here are Blackman, Hanning, and Hamming. Other windows	can be
	added by doing the following:
		1. "Window type constants" section: Define a global constant for the new window.
		2. wsfirLP() function: Add the new window as a case in the switch() statement.
		3. Create the function for the window
		   
		   For windows with a design parameter, such as Kaiser, some modification
		   will be needed for each function in order to pass the parameter.
*/
#ifndef WSFIR_H
#define WSFIR_H

#include <math.h>

// Function prototypes
void wsfirLP(double h[], const int &N, const int &WINDOW, const double &fc);
void wsfirHP(double h[], const int &N, const int &WINDOW, const double &fc);
void wsfirBS(double h[], const int &N, const int &WINDOW, const double &fc1, const double &fc2);
void wsfirBP(double h[], const int &N, const int &WINDOW, const double &fc1, const double &fc2);
void genSinc(double sinc[], const int &N, const double &fc);
void wBlackman(double w[], const int &N);
void wHanning(double w[], const int &N);
void wHamming(double w[], const int &N);

// Window type contstants
const int W_BLACKMAN = 1;
const int W_HANNING = 2;
const int W_HAMMING = 3;

// Generate lowpass filter
// 
// This is done by generating a sinc function and then windowing it
void wsfirLP(double h[],		// h[] will be written with the filter coefficients
			 const int &N,		// size of the filter (number of taps)
			 const int &WINDOW,	// window function (W_BLACKMAN, W_HANNING, etc.)
			 const double &fc)	// cutoff frequency
{
	int i;
	double *w = new double[N];		// window function
    double *sinc = new double[N];	// sinc function
    
	// 1. Generate Sinc function
	genSinc(sinc, N, fc);
    
	// 2. Generate Window function
	switch (WINDOW) {
		case W_BLACKMAN:	// W_BLACKMAN
			wBlackman(w, N);
			break;
		case W_HANNING:		// W_HANNING
			wHanning(w, N);
			break;
		case W_HAMMING:		// W_HAMMING
			wHamming(w, N);
			break;
		default:
			break;
	}

	// 3. Make lowpass filter
	for (i = 0; i < N; i++) {
		h[i] = sinc[i] * w[i];
	}

	// Delete dynamic storage
	delete []w;
	delete []sinc;

	return;
}

// Generate highpass filter
//
// This is done by generating a lowpass filter and then spectrally inverting it
void wsfirHP(double h[],		// h[] will be written with the filter coefficients
			 const int &N,		// size of the filter
			 const int &WINDOW,	// window function (W_BLACKMAN, W_HANNING, etc.)
			 const double &fc)	// cutoff frequency
{
	int i;

	// 1. Generate lowpass filter
	wsfirLP(h, N, WINDOW, fc);

	// 2. Spectrally invert (negate all samples and add 1 to center sample) lowpass filter
	// = delta[n-((N-1)/2)] - h[n], in the time domain
	for (i = 0; i < N; i++) {
		h[i] *= -1.0;	// = 0 - h[i]
	}
	h[(N-1)/2] += 1.0;	// = 1 - h[(N-1)/2]

	return;
}

// Generate bandstop filter
//
// This is done by generating a lowpass and highpass filter and adding them
void wsfirBS(double h[],		// h[] will be written with the filter taps
			 const int &N,		// size of the filter
			 const int &WINDOW,	// window function (W_BLACKMAN, W_HANNING, etc.)
			 const double &fc1,	// low cutoff frequency
			 const double &fc2)	// high cutoff frequency
{
	int i;
	double *h1 = new double[N];
	double *h2 = new double[N];

	// 1. Generate lowpass filter at first (low) cutoff frequency
	wsfirLP(h1, N, WINDOW, fc1);

	// 2. Generate highpass filter at second (high) cutoff frequency
	wsfirHP(h2, N, WINDOW, fc2);

	// 3. Add the 2 filters together
	for (i = 0; i < N; i++) {
		h[i] = h1[i] + h2[i];
	}

	// Delete dynamic memory
	delete []h1;
	delete []h2;

	return;
}

// Generate bandpass filter
//
// This is done by generating a bandstop filter and spectrally inverting it
void wsfirBP(double h[],		// h[] will be written with the filter taps
			 const int &N,		// size of the filter
			 const int &WINDOW,	// window function (W_BLACKMAN, W_HANNING, etc.)
			 const double &fc1,	// low cutoff frequency
			 const double &fc2)	// high cutoff frequency
{
	int i;

	// 1. Generate bandstop filter
	wsfirBS(h, N, WINDOW, fc1, fc2);

	// 2. Spectrally invert bandstop (negate all, and add 1 to center sample)
	for (i = 0; i < N; i++) {
		h[i] *= -1.0;	// = 0 - h[i]
	}
	h[(N-1)/2] += 1.0;	// = 1 - h[(N-1)/2]

	return;
}

// Generate sinc function - used for making lowpass filter
void genSinc(double sinc[],		// sinc[] will be written with the sinc function
			 const int &N,		// size (number of taps)
			 const double &fc)	// cutoff frequency
{
	int i;
	const double M = N-1;
	double n;

	// Constants
	const double PI = 3.14159265358979323846;

	// Generate sinc delayed by (N-1)/2
	for (i = 0; i < N; i++) {
		if (i == M/2.0) {
			sinc[i] = 2.0 * fc;
		}
		else {
			n = (double)i - M/2.0;
			sinc[i] = sin(2.0*PI*fc*n) / (PI*n);
		}
	}        

	return;
}

// Generate window function (Blackman)
void wBlackman(double w[],		// w[] will be written with the Blackman window
			   const int &N)	// size of the window
{
	int i;
	const double M = N-1;
	const double PI = 3.14159265358979323846;

	for (i = 0; i < N; i++) {
		w[i] = 0.42 - (0.5 * cos(2.0*PI*(double)i/M)) + (0.08*cos(4.0*PI*(double)i/M));
	}

	return;
}

// Generate window function (Hanning)
void wHanning(double w[],		// w[] will be written with the Hanning window
			  const int &N)		// size of the window
{
	int i;
	const double M = N-1;
	const double PI = 3.14159265358979323846;

	for (i = 0; i < N; i++) {
		w[i] = 0.5 * (1.0 - cos(2.0*PI*(double)i/M));
	}

	return;
}

// Generate window function (Hamming)
void wHamming(double w[],		// w[] will be written with the Hamming window
			  const int &N)		// size of the window
{
	int i;
	const double M = N-1;
	const double PI = 3.14159265358979323846;

	for (i = 0; i < N; i++) {
		w[i] = 0.54 - (0.46*cos(2.0*PI*(double)i/M));
	}

	return;
}

#endif