# Design of Digital Filters

*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*

## Design of Non-Recursive Filters using the Window Method

The design of non-recursive filters with a finite impulse response (FIR) is a frequent task in practical applications. The designed filter should approximate a desired frequency response as close as possible. First, the design of causal filters is considered. For many applications the resulting filter should have a linear phase characteristic since this results in a constant (and thus frequency independent) group delay. We therefore specialize the design to causal linear-phase filters.

### Causal Filters

Let's assume that the desired frequency characteristics of the filter are given by the frequency response (i.e. the DTFT spectrum) $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$. The corresponding impulse response is computed by its inverse discrete-time Fourier transform (IDTFT)

\begin{equation}
h_\text{d}[k] = \frac{1}{2 \pi} \int\limits_{- \pi}^{\pi} H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega}) \, \mathrm{e}^{\,\mathrm{j}\,\Omega\,k} \; \mathrm{d}\Omega
\end{equation}

In the general case, $h_\text{d}[k]$ will not be a causal FIR. The [Paley-Wiener theorem](https://en.wikipedia.org/wiki/Paley%E2%80%93Wiener_theorem) states, that a causal system $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ can only have zeros at single frequencies. This is not the case for idealized filters, like e.g. the [ideal low-pass filter](https://en.wikipedia.org/wiki/Low-pass_filter#Ideal_and_real_filters). The basic idea of the window method is to truncate the impulse response $h_\text{d}[k]$ in order to derive a causal FIR filter. This can be achieved by applying a window $w[k]$ of finite length $N$ to $h_\text{d}[k]$

\begin{equation}
h[k] = h_\text{d}[k] \cdot w[k]
\end{equation}

where $h[k]$ denotes the impulse response of the designed filter. Its frequency response $H(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ is given by the multiplication theorem of the discrete-time Fourier transform (DTFT)

\begin{equation}
H(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = \frac{1}{2 \pi} \; H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega}) \circledast W(\mathrm{e}^{\,\mathrm{j}\,\Omega})
\end{equation}

where $W(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ denotes the DTFT of the window function $w[k]$. The frequency response $H(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ of the filter is given as the periodic convolution of the desired frequency response $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ and the frequency response of the window function $W(\mathrm{e}^{\,\mathrm{j}\,\Omega})$. The frequency response $H(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ is equal to the desired frequency response $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ only if $W(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = 2 \pi \cdot \delta(\Omega)$. This would require that $w[k] = 1$ for $k = -\infty, \dots, \infty$. Hence for a window $w[k]$ of finite length deviations from the desired frequency response are to be expected.

In order to investigate the effect of truncation on the frequency response $H(\mathrm{e}^{\,\mathrm{j}\,\Omega})$, a particular window is considered. A straightforward choice is the rectangular window $w[k] = \text{rect}_N[k]$ of length $N$. Its DTFT is given as

\begin{equation}
W(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = \mathrm{e}^{-\mathrm{j} \, \Omega \,\frac{N-1}{2}} \cdot \frac{\sin(\frac{N \,\Omega}{2})}{\sin(\frac{\Omega}{2})}
\end{equation}

The frequency-domain properties of the rectangular window have already been discussed for the [leakage effect](../spectral_analysis_deterministic_signals/leakage_effect.ipynb). The rectangular window features a narrow main lobe at the cost of relative high sidelobe level. The main lobe gets narrower with increasing length $N$. The convolution of the desired frequency response with the frequency response of the window function effectively results in smoothing and ringing. While the main lobe will smooth discontinuities of the desired transfer function, the sidelobes result in undesirable ringing effects. The latter can be alleviated by using other window functions. Note that typical [window functions](../spectral_analysis_deterministic_signals/window_functions.ipynb) decay towards their ends and are symmetric with respect to their center. This may cause problems for desired impulse responses with large values towards their ends.

#### Example: Causal approximation of ideal low-pass

The design of an ideal low-pass filter using the window technique is illustrated in the following example. For $|\Omega| < \pi$ the transfer function of the ideal low-pass is given as

\begin{equation}
H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = \begin{cases}
1 & \text{for } |\Omega| \leq \Omega_\text{c} \\
0 & \text{otherwise}
\end{cases}
\end{equation}

where $\Omega_\text{c}$ denotes the corner frequency of the low-pass. An inverse DTFT of the desired transfer function yields

\begin{equation}
h_\text{d}[k] = \frac{\Omega_\text{c}}{\pi} \cdot \text{si}[\Omega_\text{c} \, k]
\end{equation}

The impulse response $h_\text{d}[k]$ is not causal nor FIR. In order to derive a causal and FIR approximation a rectangular window $w[k]$ of length $N$ is applied

\begin{equation}
h[k] = h_\text{d}[k] \cdot \text{rect}_N[k]
\end{equation}

The resulting magnitude and phase response is computed numerically in the following.


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig

N = 32  # length of filter
Omc = np.pi/2

# compute impulse response
k = np.arange(N)
hd = Omc/np.pi * np.sinc(k*Omc/np.pi)
# windowing
w = np.ones(N)
h = hd * w

# frequency response
Om, H = sig.freqz(h)

# plot impulse response
plt.figure(figsize=(10, 3))
plt.stem(h)
plt.title('Impulse response')
plt.xlabel(r'$k$')
plt.ylabel(r'$h[k]$')
# plot magnitude responses
plt.figure(figsize=(10, 3))
plt.plot([0, Omc, Omc], [0, 0, -100], 'r--')
plt.plot(Om, 20 * np.log10(abs(H)))
plt.title('Magnitude response')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|H(e^{j \Omega})|$ in dB')
plt.axis([0, np.pi, -20, 3])
plt.grid()
# plot phase responses
plt.figure(figsize=(10, 3))
plt.plot(Om, np.unwrap(np.angle(H)))
plt.title('Phase')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$\varphi (\Omega)$ in rad')
plt.grid()
```

**Exercises**

* Does the resulting filter have the desired phase?
* Increase the length `N` of the filter. What changes?

### Zero-Phase Filters

Above results show that an ideal-low pass cannot be realized very well with the window method. The reason is that an ideal-low pass has zero-phase, as most of the idealized filters.

Lets assume a general zero-phase filter with transfer function $H_d(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = A(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ with amplitude $A(\mathrm{e}^{\,\mathrm{j}\,\Omega})  \in \mathbb{R}$. Its impulse response $h_d[k] = \mathcal{F}_*^{-1} \{ H_d(e^{j \Omega} \}$ is conjugate complex symmetric

\begin{equation}
h_d[k] = h_d^*[-k]
\end{equation}

due to the symmetry relations of the DTFT. Hence, a transfer function with zero-phase cannot be realized by a causal non-recursive filter. This observation motivates to replace the zero-phase by a linear-phase in such situations. This is illustrated in the following.

### Causal Linear-Phase Filters

The design of non-recursive FIR filters with a linear phase is often desired due to their constant group delay. Let's assume a system with generalized linear phase. For $|\Omega| < \pi$ its transfer function is given as

\begin{equation}
H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) \cdot \mathrm{e}^{-\mathrm{j} \alpha \Omega + \mathrm{j} \beta}
\end{equation}

where $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) \in \mathbb{R}$ denotes the amplitude of the filter, $\alpha$ its linear phase and $\beta$ a constant phase offset. Such a system can be decomposed into two cascaded systems: a zero-phase system with transfer function $A(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ and an all-pass with phase $\varphi(\Omega) = - \alpha \Omega + \beta$. The linear phase term $- \alpha \Omega$ results in the desired constant group delay $t_g(\Omega) = \alpha$. 

The impulse response $h[k]$ of a linear-phase system shows a specific symmetry which can be deduced from the symmetry relations of the DTFT for odd/even symmetry of $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ as

\begin{equation}
h[k] = \pm h[N-1-k]
\end{equation}

for $k=0, 1, \dots, N-1$ where $N \in \mathbb{N}$ denotes the length of the (finite) impulse response. The transfer function of a linear phase filter is given by its DTFT

\begin{equation}
H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = \sum_{k=0}^{N-1} h[k] \, \mathrm{e}^{\,-\mathrm{j}\,\Omega\,k}
\end{equation}

Introducing the symmetry relations of the impulse response $h[k]$ into the DTFT and comparing the result with above definition of a generalized linear phase system reveals four different types of linear-phase systems. These can be discriminated with respect to their phase and magnitude characteristics

| Type | Length $N$ | Impulse Response $h[k]\;$ | Constant Group Delay $\alpha$ in Samples| Constant Phase $\beta$ | Transfer Function $A(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ |
| :---:  | :---: | :---: | :---:| :---: | :---: |
| 1 | odd | $h[k] = h[N-1-k]$ | $\alpha = \frac{N-1}{2} \in \mathbb{N}$ | $\beta = \{0, \pi\}$ | $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = A(\mathrm{e}^{-\,\mathrm{j}\,\Omega})$, all filter characteristics|
| 2 | even| $h[k] = h[N-1-k]$ | $\alpha = \frac{N-1}{2} \notin \mathbb{N}$ | $\beta = \{0, \pi\}$ | $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = A(\mathrm{e}^{-\,\mathrm{j}\,\Omega})$, $A(\mathrm{e}^{\,\mathrm{j}\,\pi}) = 0$, only lowpass or bandpass|
| 3 | odd | $h[k] = -h[N-1-k]$ | $\alpha = \frac{N-1}{2} \in \mathbb{N}$ | $\beta = \{ \frac{\pi}{2}, \frac{3 \pi}{2} \} $ | $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = - A(\mathrm{e}^{-\,\mathrm{j}\,\Omega})$, $A(\mathrm{e}^{\,\mathrm{j}\,0}) = A(\mathrm{e}^{\,\mathrm{j}\,\pi}) =  0$, only bandpass|
| 4 | even | $h[k] = -h[N-1-k]$ | $\alpha = \frac{N-1}{2} \notin \mathbb{N}$ | $\beta = \{ \frac{\pi}{2}, \frac{3 \pi}{2} \} $ | $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = - A(\mathrm{e}^{-\,\mathrm{j}\,\Omega})$, $A(\mathrm{e}^{\,\mathrm{j}\,0}) = 0$, only highpass or bandpass|

These relations have to be considered in the design of a causal linear phase filter. Depending on the desired frequency characteristics $A(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ the suitable type is chosen. The odd/even length $N$ of the filter and the phase (or group delay) is chosen accordingly for the design of the filter. It is also abovious that a filter with zero-phase $\alpha = 0$, e.g. an ideal low-pass, would result in $N=1$.

#### Example: Causal linear-phase approximation of ideal low-pass

We aim at the design of a causal linear-phase low-pass using the window technique. According to the previous example, the desired frequency response has an even symmetry $A(\mathrm{e}^{\,\mathrm{j}\,\Omega}) = A(\mathrm{e}^{\,-\mathrm{j}\,\Omega})$ with $A(\mathrm{e}^{\mathrm{j}\,0}) = 1$. This could be realized by a filter of type 1 or 2. We choose type 1, since the resulting filter then exhibits an integer group delay of $t_g(\Omega) = \frac{N-1}{2}$ samples. Consequently the length of the filter $N$ has to be odd. 

The impulse response $h_\text{d}[k]$ is given by the inverse DTFT of $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$ as

\begin{equation}
h_\text{d}[k] = \frac{\Omega_\text{c}}{\pi} \cdot \text{si} \left[ \Omega_\text{c} \left(k-\frac{N-1}{2} \right) \right]
\end{equation}

The impulse response fulfills the desired symmetry for $k=0,1, \dots, N-1$. A causal and FIR approximation is obtained by applying a window function of length $N$ to the impulse response $h_\text{d}[k]$

\begin{equation}
h[k] = h_\text{d}[k] \cdot w[k]
\end{equation}

Note that the window function $w[k]$ also has to fulfill the desired symmetries.

As already outlined, the chosen window determines the properties of the transfer function $H(\mathrm{e}^{\,\mathrm{j}\,\Omega})$. The [spectral properties of commonly applied windows](../spectral_analysis_deterministic_signals/window_functions.ipynb) have been discussed previously. The width of the main lobe will generally influence the  smoothing of the desired transfer function $H_\text{d}(\mathrm{e}^{\,\mathrm{j}\,\Omega})$, while the sidelobes influence the typical ringing artifacts. This is illustrated in the following.


```python
N = 33  # length of filter
Omc = np.pi/2

# compute impulse response
k = np.arange(N)
hd = Omc/np.pi * np.sinc((k-(N-1)/2)*Omc/np.pi)
# windowing
w1 = np.ones(N)
w2 = np.blackman(N)
h1 = hd * w1
h2 = hd * w2

# frequency responses
Om, H1 = sig.freqz(h1)
Om, H2 = sig.freqz(h2)

# plot impulse response
plt.figure(figsize=(10, 3))
plt.stem(h1)
plt.title('Impulse response (rectangular window)')
plt.xlabel(r'$k$')
plt.ylabel(r'$h[k]$')
# plot magnitude responses
plt.figure(figsize=(10, 3))
plt.plot([0, Omc, Omc], [0, 0, -300], 'r--')
plt.plot(Om, 20 * np.log10(abs(H1)), label='rectangular window')
plt.plot(Om, 20 * np.log10(abs(H2)), label='Blackmann window')
plt.title('Magnitude response')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|H(e^{j \Omega})|$ in dB')
plt.axis([0, np.pi, -120, 3])
plt.legend(loc=3)
plt.grid()
# plot phase responses
plt.figure(figsize=(10, 3))
plt.plot(Om, np.unwrap(np.angle(H1)), label='rectangular window')
plt.plot(Om, np.unwrap(np.angle(H2)), label='Blackmann window')
plt.title('Phase')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$\varphi (\Omega)$ in rad')
plt.legend(loc=3)
plt.grid()
```

**Exercises**

* Does the impulse response fulfill the required symmetries for a type 1 filter?
* Can you explain the differences between the magnitude responses $|H(\mathrm{e}^{\,\mathrm{j}\,\Omega})|$ for the different window functions?
* What happens if you increase the length `N` of the filter?


**Copyright**

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2017*.

24c9f36de7176a73206197798add4d91887be61b
ipynb
Jupyter Notebook
filter_design/window_method.ipynb
swchao/digitalSignalProcessingLecture
89acae62ea710211014912d61a461ca8a3d6d713
MIT
2017-11-14T16:14:37.000Z
2021-05-16T21:01:41.000Z
filter_design/window_method.ipynb
ganlubbq/digital-signal-processing-lecture
f9ac5b2f5500aa612b48d1d920c7cba366c44dba
MIT
filter_design/window_method.ipynb
ganlubbq/digital-signal-processing-lecture
f9ac5b2f5500aa612b48d1d920c7cba366c44dba
MIT
2020-06-26T14:19:29.000Z
2020-12-11T08:31:29.000Z
Qwen/Qwen-72B
1. YES
2. YES