Design of IIR Allpass Filters: a least squares equation error method

The frequency response of a causal $N^{th}$ order allpass filter is given by $$H(e^{j\omega})=\frac{e^{-jN\omega}A^*(e^{j\omega})}{A(e^{j\omega})}\tag{1}$$ with $^*$ denoting complex conjugation, and with the denominator polynomial given by $$A(e^{j\omega})=1+a_1e^{-j\omega}+\ldots +a_Ne^{-jN\omega}\tag{2}$$ Note that from $(1)$, $|H(e^{j\omega})|=1$, and, consequently $$H(e^{j\omega})=e^{j\phi(\omega)}\tag{3}$$ where $\phi(\omega)$ is the filter's phase response.

The design problem is to find the coefficients $a_n$ in $(2)$ such that the phase $\phi(\omega)$ approximates a given desired phase response $\phi_d(\omega)$. The filter's phase $\phi(\omega)$ is a highly non-linear function of the filter coefficients $a_n$, so it's difficult to directly minimize some norm of the phase error $\phi(\omega)-\phi_d(\omega)$. It is more straightforward to consider the complex error function $H(e^{j\omega})-e^{j\phi_d(\omega)}$. This error function still depends non-linearly on the filter coefficients. However, multiplying this error function by the denominator of $H(e^{j\omega})$ results in an error function which is a linear function of the filter coefficients: $$\epsilon(\omega)=e^{-jN\omega}A^*(e^{j\omega}) - A(e^{j\omega})e^{j\phi_d(\omega)}\tag{4}$$ This error function is usually called "equation error", and it is also used for designing general IIR filters. In this DSP stackexchange post I describe the basics of that method for the general case, and there are also some relevant links.

Minimization of the error function $(4)$ is equivalent to solving the following overdetermined system of equations: $$e^{-jN\omega_i}A^*(e^{j\omega_i})\overset{!}{=}A(e^{j\omega_i})e^{j\phi_d(\omega_i)},\quad i=1,\ldots,M\tag{5}$$ where we've chosen a dense frequency grid $\omega_i$ in the interval $[0,\pi]$ such that there are considerably more equations than unknown filter coefficients. Since the equations in $(5)$ are linear functions of the filter coefficients $a_n$, we can rewrite $(5)$ using matrix/vector notation: $${\mathbf C}{\mathbf a}\overset{!}{=}{\mathbf d}\tag{6}$$ where $\mathbf C$ is a complex valued $M\times N$ matrix with rows $${\mathbf c}_i=\left[e^{-j(N-1)\omega_i}-e^{j(\phi_d(\omega_i)-\omega_i)},e^{-j(N-2)\omega_i}-e^{j(\phi_d(\omega_i)-2\omega_i)},\ldots,1-e^{j(\phi_d(\omega_i)-N\omega_i)}\right]\tag{7}$$ $\mathbf a$ is the $N\times 1$ vector of unknown filter coefficients $a_n$, and $\mathbf d$ is a complex $M\times 1$ vector with elements $e^{j\phi_d(\omega_i)}-e^{-jN\omega_i}$.

We can solve $(6)$ in a least squares sense either by solving the corresponding normal equations or we can directly use the overdetermined system $(6)$ as the input to a QR-solver, as provided by the backslash operator in Matlab or Octave. We could compute the filter coefficients from $(6)$ but we'll add a minor generalization in terms of a weight function which allows to weigh the error differently in different frequency regions. We define a positive weight function $W(\omega)$ on the frequency grid and we define an $M\times M$ diagonal matrix $\mathbf W$ with elements $W(\omega_i)$ on its diagonal. The overdetermined system we are going to solve is then $${\mathbf W}^{1/2}{\mathbf C}{\mathbf a}\overset{!}{=}{\mathbf W}^{1/2}{\mathbf d}\tag{8}$$ We use the square root of the weights in $(8)$ just to be consistent with other weighted least squares methods. Since we're only interested in real-valued solutions, we rewrite $(8)$ as $$\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf C}_R\\{\mathbf W}^{1/2}{\mathbf C}_I\end{array}\right]{\mathbf a}\overset{!}{=}\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf d}_R\\{\mathbf W}^{1/2}{\mathbf d}_I\end{array}\right]\tag{9}$$ where the subscripts $R$ and $I$ denote real and imaginary parts, respectively. The overdetermined system of linear equations $(9)$ is the one that we need to solve in a least squares sense in order to compute the filter coefficients that minimize the weighted equation error.

It should be noted that the method described here does not guarantee the stability of the designed filter. If the resulting filter turns out to be unstable, reflecting poles about the unit circle is not an option, because that would ruin the phase approximation. Instead, a new design has to be done whereby the desired phase response is modified by a linear phase term, i.e., by adding (or subtracting) a constant delay. In this way, it is always possible to obtain a stable filter. The value of that delay, however, needs to be determined experimentally.

I'll show a few design examples to illustrate the method described above. The first design is a fractional delay. The goal is to design an allpass filter with a fractional delay of $7.5$ samples. We choose a filter order $N=8$ and no weighting, i.e., $W(\omega)=1$ for all frequencies. The figures below shows the result of this design. The top figure shows the phase error in radians. The phase error is smaller than $0.05$ radians for more than $90$% of the bandwidth. Only close to the Nyquist frequency does the phase error become relatively large. The bottom figure shows the group delay, which approximates the desired $7.5$ samples quite well for about $85$% of the bandwidth.

In the next design example, the allpass filter is used to equalize the phase distortion of a digital Cauer lowpass filter. The desired phase response of the allpass filter is the negative phase of the Cauer filter with some linear phase subtracted, such that the total desired phase response can be approximated well by a causal filter. I used a $5^{th}$ order Cauer lowpass and a $12^{th}$ order allpass equalizer. A good result was achieved with a desired total delay of $10$ samples. The figure below shows the design result. The top figure shows the phase approximation error in the passband, and the bottom plot shows the group delay functions in the passband for the Cauer filter (blue), the allpass equalizer (orange), and the cascade of Cauer filter and allpass filter (yellow). The total group delay is close to the specified $10$ samples over a large range of the passband.

The final example is the design of a Hilbert transformer. A Hilbert transformer shifts the phase of the input signal by $-\pi/2$ while keeping the amplitude unchanged. An ideal Hilbert transformer is a non-causal system, that's why a delay has to be added to the desired phase such that the approximation error of a causal allpass becomes reasonably small. I chose the filter order $N=10$ and a delay of $9$ samples. The top figure below shows the phase error. As is to be expected, the approximation is more difficult close to DC and Nyquist. However, over a large range of the bandwidth, the resulting filter computes an accurate Hilbert transform. The bottom figure shows the group delay, which approximates the desired value of $9$ samples, with larger errors at low and high frequencies.

A simple Matlab/Octave implementation of this design method can be found on GitHub: iir_ap.m

Digital Filter Design: An Overview

I wrote this little article as an answer to a question on dsp.stackexchange.com and I think it gives a good overview for someone who is just starting with the topic. Enjoy reading and leave comments if you wish.

Digital filter design is a very large and mature topic and there is a lot of material available. What I want to try here is to get a beginner started and to make the existing material more accessible. Instead of digital filters I should actually be talking about discrete-time filters because I will not consider coefficient and signal quantization here.

A non-recursive linear time-invariant (LTI) filter can be described by the following difference equation $$y(n)=h_0x(n)+h_1x(n-1)+\ldots +h_{N-1}x(n-N+1)=\\=\sum_{k=0}^{N-1}h_kx(n-k)\tag{1}$$ where $y(n)$ is the output sequence, $x(n)$ is the input sequence, $n$ is the time index, $h_k$ are the filter coefficients, and $N$ is the filter length (the number of taps). The filter taps $h_k$ are also the impulse response of the filter because if the input signal is an impulse, i.e. $x(n)=\delta(n)$, then $y(n)=h_n$ (if the filter's memory has been initialized with zeros). Equation (1) describes a linear time-invariant finite impulse response (FIR) system. The sum on the right-hand side of (1) is a convolution sum, i.e. the output signal is obtained by convolving the input signal with the impulse response. This is always true, but for IIR filters we cannot explicitly compute the convolution sum because the impulse response is infinitely long, i.e. there are infinitely many coefficients $h_k$. One important advantage of FIR filters is that they are always stable, i.e. for a bounded input sequence, the output sequence is always bounded. Another advantage is that FIR filters can always be realized with an exactly linear phase, i.e. they will not add any phase distortion apart from a pure delay. Furthermore, the design problem is usually easier, as we will see later.

A recursive LTI filter is described by the following difference equation: $$y(n)=b_0x(n)+b_1x(n-1)+\ldots+b_Mx(n-M)-\\ -a_1y(n-1)-\ldots-a_Ny(n-N)\tag{2}$$
Equation (2) shows that the output is not only composed of weighted and delayed input samples, but also of weighted past output samples. In general, the impulse response of such a system is infinitely long, i.e. the corresponding system is an IIR (infinite impulse response) system. However, there are special cases of recursive filters with a finite impulse response. Note that the impulse response is not anymore given by either the coefficients $b_k$ or $a_k$ as in the case of FIR filters. One advantage of IIR filters is that steep filters with high stopband attenuation can be realized with much fewer coefficients (and delays) than in the FIR case, i.e. they are computationally more efficient. However, one needs to be careful with the choice of the coefficients $a_k$ because IIR filter can be unstable, i.e. their output sequence can be unbounded, even with a bounded input sequence.

Filters can be designed according to specifications either in the time (sample) domain or in the frequency domain, or both. I'll restrict myself to specifications in the frequency domain. In this case one needs to consider the frequency responses of FIR and IIR systems. The frequency response of a system is the Fourier transform of its impulse response, assuming that it exists (which is the case for stable systems). The frequency response of an FIR filter is $$H(e^{j\theta})=\sum_{k=0}^{N-1}h_ke^{-jk\theta}\tag{3}$$ where $\theta$ is the discrete-time frequency variable: $$\theta=\frac{2\pi f}{f_s}$$ with the actual frequency $f$ and the sampling frequency $f_s$. From (3) you can see that approximating a desired frequency response by an FIR system is basically a problem of polynomial approximation. For recursive systems we have $$H(e^{j\theta})=\frac{\sum_{k=0}^Mb_ke^{-j\theta}}{1+\sum_{k=1}^Na_ke^{-j\theta}}\tag{4}$$ and you get a rational approximation problem, which is usually much more difficult than the polynomial approximation problem in the case of FIR filters. From (3) and (4) you can see that the frequency response of an FIR filter is of course only a special case of the response of a recursive filter with coefficients $a_k=0$, $k=1,\dots,N$.

Let's now take a quick look at filter design methods. For FIR filters you could take an inverse Fourier transform of the desired frequency response to get the impulse response of the filter, which directly corresponds to the filter coefficients. Since you approximate the desired response by a finite length impulse response you should apply a smooth window to the obtained impulse response to minimize oscillations in the actual frequency response due to Gibbs' phenomenon. This method is called frequency-sampling method. For simple standard filters like ideal lowpass, highpass, bandpass or bandstop filters (and a few others), you could even analytically calculate the exact impulse response by taking the inverse Fourier transform of the ideal desired response: $$h_k=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\theta})e^{jk\theta}d\theta$$ This integral is easy to evaluate for piecewise constant desired responses, as is the case for ideal frequency-selective filters. This will give you an infinitely long, non-causal impulse response, which needs to be windowed and shifted to make it finite and causal. This method is know as window-design. There are of course many other FIR filter design methods. One important numerical method is the famous Parks-McClellan exchange algorithm which designs optimal filters with constant passband and stopband ripples. It is a numerical approximation method and there are many software implementations available, e.g. in Matlab and Octave. The most common IIR design method for frequency selective filters is the bilinear transformation method. This method simply uses analytical formulas for the design of optimal analog filters (such as Butterworth, Chebyshev, Cauer/elliptic, and Bessel filters), and transforms them to the discrete-time domain by applying a bilinear transformation to the complex variable $s$ (analog domain) which maps the (imaginary) frequency axis of the complex $s$-plane to the unit circle in the complex $z$-plane (discrete-time domain). Don't worry if you do not yet know much about complex transfer functions in the analog or discrete-time domain because there are good implementations available of the bilinear transform method, e.g. in Matlab or Octave.

There are of course many more interesting and useful methods, depending on the type of specifications, but I hope that this will get you started and will make any material you come across more understandable. A very good (and free) book covering some basic filter design methods (and a lot more) is Introduction to Signal Processing by Orfanidis. You can find several design examples there. Another great classic book is Digital Filter Design by Parks and Burrus.