Subsections


FIR bandpass filter

The purpose of this appendix is to make an algorithm that implements a $ M$ -order FIR bandpass filter. The algorithm will be tested in MATLABTM and later implemented in C.

Principle

The FIR filter is designed using windowing. This method is further described in [Oppenheim & Schafer, 1998, page 465 - 477].

The method is to make an ideal filter in the frequency domain, and then translate it into the discrete time domain. However this will give an infinite impulse response. To compensate for this, a window function is multiplied onto the ideal impulse response.

In this appendix a rectangular window is chosen. This is done for simplicity as the window amplitude is either 1 or 0. In FIR filter design the order for the filter is denoted $ M$ and it determines the length of the window. The type of FIR filter used here utilizes linear phase i.e. constant group delay.

An example of a FIR bandpass filter is shown in figure M.1. The filter has a center frequency ($ f_c$ ) of $ 5\textrm{ kHz}$ and a bandwith of 2 kHz. The order ($ M$ ) is 20.

Figure M.1: 20 order FIR bandpass filter.
\includegraphics[width=10cm]{appendix/firFilter/bild/firplotca5kbw2km20.eps}

Mathematical description

Type I FIR filter design by windowing utilizes the ideal form of figure jtkr04:fig:fir02. Examining the ideal filter in the frequency domain, the definition of $ H(e^{j\omega})$ could be written as a modified equation (7.42) from [Oppenheim & Schafer, 1998, p.466]


$\displaystyle H(e^{j\omega} ) = \left\{ \begin{array}{ll} e^{j\omega} &
\left\{...
...c + \omega_b
\end{array} \right. \\
0 & \mathrm{Otherwise} \end{array} \right.$     (M.1)

Figure M.2: Ideal bandpass filter.
\includegraphics{appendix/firFilter/bild/firPerfectForm.eps}

Corresponding to the discrete-time notation of

$\displaystyle h[n] = h_d[n] \cdot w[n]$     (M.2)

with $ h_d[n]$ as the approximated impulse response of the filter, and with $ w[n]$ as the windowing function
$\displaystyle w[n] = \left\{ \begin{array}{ll} 1, & \mathrm{for}  0 \leq n \leq M \\
0, & \mathrm{otherwise} \end{array} \right.$     (M.3)

The product of (M.2) is in the frequency domain equal to the convolution
$\displaystyle H(e^{j\omega}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} H_d(e^{j\omega})\cdot W(e^{j\omega - \phi}) d\phi$     (M.4)

Utilizing equation (7.52) with $ h_d[M-n] = h_d[n]$ and (5.146b) from [Oppenheim & Schafer, 1998, p. 472 and p. 297 respectively] with the amplitude function set to 1, $ H_d(e^{j\omega }$ can be written as

$\displaystyle H_d(e^{j\omega}) = \left\{ \begin{array}{ll} e^{-j\omega \frac{M}...
...c + \omega b
\end{array} \right. \\
0 & \mathrm{Otherwise} \end{array} \right.$     (M.5)

Transforming $ H_d(e^{j\omega})$ to the time domain yields (by use of Euler rewriting)

$\displaystyle h_d[n]$ $\displaystyle =$ $\displaystyle \frac{1}{2\pi} \int_{-\pi}^{\pi} H_d(e^{j\omega}) e^{j\omega n} d\omega$ (M.6)
  $\displaystyle =$ $\displaystyle \frac{1}{2\pi} \int_{-\omega_c - \omega_b}^{-\omega_c + \omega_b}...
...- \omega_b}^{\omega_c + \omega_b} e^{-j\omega\frac{M}{2}} e^{j\omega n} d\omega$ (M.7)
  $\displaystyle =$ $\displaystyle \frac{1}{\pi(n-\frac{M}{2})} \cdot \left[\sin\left((\omega_c + \omega_b)(n-\frac{M}{2})\right) \ldots \right.$  
  $\displaystyle \ldots$ $\displaystyle \left. - \sin\left((\omega_c - \omega_b)(n-\frac{M}{2})\right) \right]$ (M.8)

Multiplying $ h_d[n]$ with the window $ w[n]$ in the time domain, to achieve the definition of $ h[n]$ , and transforming to the Z-domain yields

$\displaystyle H(z)$ $\displaystyle =$ $\displaystyle \sum_{n = -\infty}^{\infty} h_d[n]\cdot w[n]\cdot z^{-n}$ (M.9)
  $\displaystyle =$ $\displaystyle \sum_{n = 0}^{M} h_d[n]\cdot z^{-n}$ (M.10)

due to the magnitude of $ w[n]$ being 1 in $ 0 \leq n \leq M$ and zero otherwise. This is an approximation of the ideal bandpass filter of equation (M.1) and figure jtkr04:fig:fir02.

Using the basic attribute of $ Y(z) = H(z)\cdot X(z)$ and transforming back to the time domain, will yield the filter output

$\displaystyle y[i] = \sum_{n = 0}^{M} h_d[n]\cdot x[i-n]$     (M.11)

The above can be readily implemented in an algorithm, with the value of $ M$ describing the order of the filter. The values of $ h_d[n]$ are reusable for each $ n$ , calculated from the expression in (M.8).

Digitization considerations

The designed FIR filter is specifically intended for digital implementation, and thereby doesn't have any particular implementation considerations, except that according to [Oppenheim & Schafer, 1998, page 386] at least 14 bits must by used for the filter coefficients.

MATLABTM implementation

The MATLABTM implementation relies on the following information to be supplied: The filter algorithm does the following: Note that the first $ M$ samples are using samples before $ t=0$ as zeros. Tested against MATLABTM's internal commands (FIR1), the designed filter produces the same coefficients as MATLABTM.

C implementation considerations

The C-implementation on the selected PC platform imposes the consideration of datatypes and data length.

The system datatype used for handling soundcard samples is set to integers, while the internal workings of the filter relies on float values. The solution is to use float values internally, and export values rounded to integers.

As C contains no built in ``length'' command, the length of the supplied samples must be specified. Furthermore, the data to be processed is supplied in chunks. This results in problems as the filter uses $ M$ samples from the past. When crossing the boundary between chunks this must be addressed.

Conclusion

In this appendix an $ M$ order FIR bandpass filter has been implemented. The order of the filter can be varied according to the demands of a given situation. However the zeroing of the samples before $ t=0$ must be handled for the actual C-implementation.

Jes Toft Kristensen 2005-12-13