The Window Method for 
FIR Digital Filter Design

Music 420 / EE 265
	Winter 1992
	Julius O. Smith III
CCRMA, Stanford University

Problem Statement

Find a finite-length impulse response whose DTFT approximates the frequency response of the ideal lowpass filter.

Digital filter theory lies outside the scope of this course, as it is quite a large topic.  However, the convolution theorem is very much within the scope of this course, and every linear, time-invariant filter can be represented as a convolution between the input signal and the filter impulse response.  In the frequency domain, we know by the convolution theorem that the input spectrum is then multiplied pointwise by the spectrum of the impulse response (which we term the frequency response of the filter).    The impulse response is so named because it is what the filter produces as output when its input is the unit sample pulse delta[n].  So, whenever two signals are convolved, one of them can be regarded as the input signal to a filter, and the other can be regarded as the filter impulse response.

A terse reference on digital filter theory is given in handout #4.  A good reference on the topic is the book Digital Filter Design by T.W. Parks and C.S. Burrus, Wiley, 1987.

In this notebook, we examine the problem of designing an FIR lowpass filter by the window method.  We begin with a specification of the desired frequency response which is a rectangular window in the frequency domain, passing low frequencies with unit gain and rejecting high frequencies with a gain of 0.  The inverse DTFT gives us the ideal impulse response, sinc[Pi fc n], which is unfortunately of infinite duration.  We cannot implement a convolution involving an infinite-duration signal.  Therefore, we must window the impulse response.  This notebook illustrates these operations, comparing use of the rectangular window (i.e., the impulse response is simply truncated) to use of the Kaiser window (which tapers the impulse response more gracefully to 0).

Utilities

Some Graphing Tools

pointSize = 0.02;

viewCentered[X_,label_,ptSize_:pointSize]:=Block[{len=Length[X]},
	ListPlot[Transpose[
		{Join[Range[len/2]-1,Range[len/2]-len/2-1],X}],
		PlotStyle->{PointSize[ptSize]}, 
		PlotRange->All, PlotLabel->label]]

view[X_,label_,ptSize_:pointSize]:=ListPlot[Transpose[
	{Range[Length[X]]-1,X}],
		PlotStyle->{PointSize[ptSize]}, 
		PlotRange->All, PlotLabel->label]

Continuous Spectra on the Unit Circle

continuousRect[width_,f_] := If[Abs[f]>(width/2), 0.0, 1.0];

plotContinuousHz[spectrum_,label_] := 
	Plot[spectrum[f], {f,-fs/2,fs/2},
		Axes->{-fs/2,0},
		(* AxesOrigin->{-fs/2, (v2.0) *)
		PlotRange->All, PlotLabel->label];

Indexed Lists of the form {{n0,y0},{n0+1,y1},...}

viewIndexedList[X_,label_]:=ListPlot[X,
		PlotStyle->{PointSize[pointSize]},
		Axes->{X[[1]][[1]],0},
		(* AxesOrigin->{X[[1]][[1]],0}, (v2.0) *)
		PlotRange->All, PlotLabel->label]

Clear[zeroPhaseIndexedListToFFTInput];
zeroPhaseIndexedListToFFTInput[h_] := Block[
(* Assumes FFT size NN is defined.  
   Takes a list of length less than NN. *)
    {len,halfLen,hv,hRight,hLeft,zeroPadding},
	len = Length[h];
	halfLen = (len-1)/2;
	hv = Transpose[h][[2]];
	hRight = Table[hv[[n]],{n, halfLen+1,len}];
	hLeft = Table[hv[[n]],{n,1,halfLen}];
	zeroPadding = 0 Range[NN-len];
	Return[Join[hRight,zeroPadding,hLeft]]];

multiplyIndexedLists[a_,b_] := Block[
	{indices,op1,op2,product},
	indices = Transpose[a][[1]];
	op1 = Transpose[a][[2]];
	op2 = Transpose[b][[2]];
	product = Times[op1,op2];
	Return[Transpose[{indices,product}]]];

FFT

fft[x_] := Chop[N[InverseFourier[x] Sqrt[Length[x]]]];

Kaiser Window

goodBeta[a_] := If[a<21,0,
		  If[a>50, 0.1102 (a-8.7),
		  	0.5842 (a-21)^0.4 + 0.07886 (a-21)]];

M = 41		(* window length *)
nl = -(M-1)/2;	(* Lowest window time index *)
nu = (M-1)/2;	(* Highest window time index *)

I0[x_] := BesselI[0,x];

kaiserUN[n_] := If[ n < nl || n > nu, 0,
		N[ I0[beta Sqrt[1 - (n/nu)^2]]/I0[beta]] ];

kscale := 1/Sum[kaiserUN[n],{n,nl,nu}];

kaiser[n_] := N[ kscale kaiserUN[n] ];

Ideal Lowpass Filter Frequency Response

fs = 1.0;	(* sampling rate in Hz *)
fc = fs/4	(* lowpass cutoff frequency in Hz *)

H[f_] = continuousRect[2 fc/fs,f];

Let's plot an example

plotContinuousHz [H,"Ideal Lowpass Filter Frequency Response"];

Ideal Lowpass Filter Impulse Response

The ideal impulse response is the inverse Fourier transform of the ideal frequency response, or, in the discrete-time case, the inverse DTFT.

h[t_] = Integrate[H[f] Exp[I 2 Pi f t / fs], {f, -fs/2, fs/2}]

= Integrate[Exp[I 2 Pi f t / fs], {f, -fc/2, fc/2}]

It is possible to extend Mathematica to handle simple integrals such as the above, but we'll help it out for now.

Clear[h];
h[t_] = Integrate[Exp[I 2 Pi f t / fs], {f,-fc/2,fc/2}]

Plot[h[t],{t,-40,40},PlotRange->All, Axes->{-40,0},
	PlotLabel->"Ideal Lowpass Filter Impulse Response"];

It is also possible to extend Mathematica with transformation rules such as Euler's identity, to simplify the above expression, but again we'll just enter the simplification manually here.  See the Georgia Tech. DSP package for more sophisticated simplification rules useful in signal processing.

Clear[h];
h[t_] = If[ t==0, fc, Chop[N[ Sin[ Pi fc t] / (Pi t) ]] ];

Plot[h[t] /. h[0]->fc/fs,{t,-40,40},
	PlotRange->All, Axes->{-40,0},
	PlotLabel->"Ideal Lowpass Filter Impulse Response"];

If t above is continuous, then we are looking at the inverse Fourier transform of the continuous rect function in the frequency domain (frequency ranges from minus infinity to infinity).  If we evaluate t only on the integers, i.e., sample it with a sampling rate of fs=1 sample per second, then we obtain the inverse DTFT of the continuous rect function in the frequency domain (frequency ranging from minus Pi to Pi).  We are are free to sample without error because the rect function is fully contained in the interval [-Pi,Pi) so that there is no aliasing.  In other words the sinc function is truly bandlimited which means we can sample it at intervals less than its reciprocal bandwidth without losing any information.

A problem with the above impulse response is that it has infinite duration.  (We arbitrarily truncated it for the plot.)  We must therefore window the impulse response to make it finite.

Kaiser Window Function

Refer to the notebook on the Kasier window.  Here we simply use it.

NN = 1024; 	(* FFT size for display purposes *)
att = 80;	(* Desired stop-band attenuation *)

beta = goodBeta[att]

w = Table[{n,kaiserUN[n]},{n,nl,nu}];

kaiserPlot = viewIndexedList[w,"Kaiser Window"];

Rectangularly Windowed Ideal Impulse Response

hrect = Table[{n,h[n]},{n,nl,nu}];

rectPlot = viewIndexedList[hrect,
	"Rectangular Windowed Ideal Impulse Resp"];

Kaiser-Windowed Ideal Impulse Response

hkaiser = multiplyIndexedLists[hrect,w];

kaiserPlot = viewIndexedList[hkaiser,
     "Kaiser-Windowed Ideal Impulse Response"];

Resulting Frequency Response

Rectangular Window

(* hrectPadded = zeroPhaseWindowToFFTInput[hrect]; bugged *)
	Clear[h];
	h = hrect;
	len = Length[h];
	halfLen = (len-1)/2;
	hRight = Table[h[[n]][[2]],{n, halfLen+1,len}];
	hLeft = Table[h[[n]][[2]],{n,1,halfLen}];
	zeroPadding = 0 Range[NN-len];
	hrectPadded = Join[hRight,zeroPadding,hLeft];

Hrect = fft[hrectPadded];
HrectDB = 20 Log[10,Abs[Hrect]];

HrectPlot = viewCentered [HrectDB,
	"Frequency Response : Rectangular Window",.008];

Kaiser Window

(* hKaiserPadded = zeroPhaseWindowToFFTInput[hkaiser]; bugged *)
	Clear[h];
	h = hkaiser;
	len = Length[h];
	halfLen = (len-1)/2;
	hRight = Table[h[[n]][[2]],{n, halfLen+1,len}];
	hLeft = Table[h[[n]][[2]],{n,1,halfLen}];
	zeroPadding = 0 Range[NN-len];
	hKaiserPadded = Join[hRight,zeroPadding,hLeft];

Hkaiser = fft[hKaiserPadded];
HkaiserDB = 20 Log[10,Abs[Hkaiser]];

viewCentered[HkaiserDB,
	"Frequency Response : Kaiser Window",.008];

We see that Kaiser's formula for beta is quite good here.  However, if the cut-off frequency fc is made less than half the main lobe width, we approach the results found previously in the notebook on the Kaiser window.   Why not use this notebook to try this out for yourself?