MISCELLANEOUS TECHNICAL ARTICLES BY Dr A R COLLINS

FIR Filter Design

Window method of FIR filter design

The basic idea behind the Window method of filter design is that the ideal frequency response of the filter is equal to 1 for all the pass band frequencies, and equal to 0 for all the stop band frequencies. The filter impulse response is obtained by taking the Discrete Fourier Transform (DFT) of the ideal frequency response. Unfortunately, the filter response would be infinitely long since it has to reproduce the infinitely steep discontinuities at the band edges in the ideal frequency response. To create a Finite Impulse Response (FIR) filter, the time domain filter coefficients must be restricted in number by multiplying by a window function of a finite width. The simplest window function is the rectangular window which corresponds to truncating the sequence after a certain number of terms.

Rectangular windowing in the time domain will result in a frequency spectrum with the width of the pass band close to the desired value but with side lobes appearing at the band edges (the effects of time domain windowing on the frequency spectrum are discussed in more detail in the Spectrum Analyser page). To suppress the side lobes and make the filter frequency response approximate more closely to the ideal, the width of the window must be increased and the window function tapered down to zero at the ends. This will increase the width of the transition region between the pass and stop bands but will lower the side lobe levels outside the pass band.

Kaiser-Bessel filter generator

The Kaiser-Bessel window function is simple to calculate and its parameters can be adjusted to produce the desired maximum side lobe level for a near minimal filter length. To demonstrate the power and simplicity of this technique, a Kaiser-Bessel filter generator is provided below.

To use it, set the sample rate (1kHz < Fs < 1MHz) and the type of filter desired; low pass, band pass or high pass, then set the number of points in the filter (N < 500) then set the frequency of ideal filter edges (Fa, Fb) and the minimum attenuation (Att) required in the stop band. Then press the "CALCULATE FILTER" button. The filter coefficients are calculated and plotted along with a graph of the frequency response of the filter. The actually algorithm and the JavaScript code to implement it are presented at the bottom of the page.

Kaiser FIR Design
FIR Listing

Kaiser-Bessel filter design formulae

The methods used in this FIR generator were taken from the paper by J. F. Kaiser, "Nonrecursive digital filter design using I0-sinh window function". In this paper Kaiser presented empirical formulae for calculating the shape parameter of the Kaiser-Bessel window required to achieve a desired stop band side lode attenuation.

The values of window shape parameter \(\alpha\) needed to achieve the attenuation \( Att \) may be calculated from:

$$\alpha = \left\{ \begin{array}{l l} 0.1102(Att - 8.7) & \quad Att > 50\\ 0.5842(Att - 21)^{0.4} + 0.07886(Att - 21) & \quad 21 \le Att \le 50\\ 0 & \quad Att < 21 \end{array}\right.$$

The transition band \(dF\) is defined as the difference in frequency from the edge of the pass band, where ripple is typically <0.2dB, and the edge of stop band, where the side lobes fall below the specified level. The edge frequency of the ideal filter is in the center of this transition region. The relation between the transition bandwidth and the filter length \(M\), is given by:

$$M \ge \frac{Att - 8}{14.36\ dF/{Fs}}$$

The formula to actually generate the Kaiser-Bessel window coefficients is as follows: $$w[j] = \frac {I_0\left(\alpha\sqrt{1-\left(\displaystyle\frac{j-N\!p}{N\!p}\right)^2}\right)}{I_0(\alpha)} \qquad \text{for}\ \ 0 \le j < M$$

where: \(M\) is the number of points in the filter,
\(N\!p = (M - 1)/2\),
\(\alpha\) is the Kaiser-Bessel window shape factor,
\(I_0()\) is the \(0^{th}\) order Bessel function of the first kind.

The two band-edge frequencies, Fa and Fb, determine the shape of the sinc function impulse response of the ideal filter. This sinc function may be calculated analytically. The formulae for this impulse response is: $$A[j] = \frac{\sin \left(2\pi j \displaystyle\frac{F\!b}{F\!s}\right) - \sin \left(2\pi j \displaystyle\frac{F\!a}{F\!s}\right)}{\pi j}$$

Truncating this impulse response to just M points and suppressing the side lobes to the required level is achieved by windowing the ideal impulse response \(A[\ ]\), by the Kaiser-Bessel window \(w[\ ]\). The required filter coefficients are given by: $$H[j] = A[j] w[j] \qquad \text{for} \ \ 0 \le j < M$$

Kaiser-Bessel filter design source code

The JavaScript code fragments to implement these formulae are shown below. The symmetry of the window function about the center point (j=Np) means that the window coefficients need only be calculated for 0<j<=Np.

JavaScript code to calculate the impulse response of the ideal filter.

A[0] = 2*(Fb-Fa)/Fs;
for(j=1; j<=Np; j++)
{
  A[j] = (Math.sin(2*j*pi*Fb/Fs)-Math.sin(2*j*pi*Fa/Fs))/(j*pi);
}

The JavaScript code to calculate the desired shape Kaiser-Bessel window.

if (Att<21)
{
  Alpha = 0;
}
else if (Att>50)
{
  Alpha = 0.1102*(Att-8.7);
}
else
{
  Alpha = 0.5842*Math.pow((Att-21), 0.4)+0.07886*(Att-21);
}
// Window the ideal response with the Kaiser-Bessel window
Inoalpha = Ino(Alpha);
for (j=0; j<=Np; j++)
{
  H[Np+j] = A[j]*Ino(Alpha*Math.sqrt(1-(j*j/(Np*Np))))/Inoalpha;
}
for (j=0; j<Np; j++)
{
  H[j] = H[M-1-j];
}

The code to generate the \(0^{th}\) order Bessel function used in the code above is:

function Ino(x)
{
  /*
   * This function calculates the zeroth order Bessel function
   */
  var d = 0, ds = 1, s = 1;
  do
  {
    d += 2;
    ds *= x*x/(d*d);
    s += ds;
  }
  while (ds > s*1e-6);
  return s;
}

These algorithms are all that is needed to generate the filter coefficients for "low pass", "high pass" and "band pass" filters. The "band stop" filter type requires one additional manipulation. The coefficients start out the same as the equivalent "band pass" filter but then the center coefficient \(H[N\!p]\), is replaced with \(1 - H[N\!p]\), and all the other coefficients are inverted \(H[j] = -H[j]\).

The complete JavaScript source code of the functions is all in the file dspUtils-14.js