FIR Filters with prescribed magnitude and phase responses: weighted least squares in the frequency domain

In this post I'm going to explain the least squares design of finite impulse response (FIR) filters in the frequency domain. The least squares criterion is a practically useful criterion (minimizing error energy) and, in the case of FIR filters, it also leads to a simple problem statement the solution of which can be computed by solving a system of linear equations.

The least squares design of linear phase filters is a special case, and its solution is implemented in the Matlab/Octave function $\texttt{firls.m}$. Here I want to describe the general case of approximating a complex desired frequency response, i.e., a response that prescribes not only the magnitude but also the phase response of the filter. If the prescribed phase response is linear and if the corresponding delay matches the chosen filter length, then the resulting filter will have a linear phase. Otherwise, the designed filter will have a non-linear phase response approximating the prescribed phase response.

The described method cannot only design filters with real-valued coefficients but it can also design FIR filters with complex coefficients. This is simply achieved by prescribing a complex desired frequency response that is not conjugate symmetric. Let $D(e^{j\omega})$ denote the complex desired frequency response. A complex FIR filter is obtained if $$D(e^{j\omega})\neq D^*(e^{-j\omega})$$ where $^*$ denotes complex conjugation. If $D(e^{j\omega})$ is conjugate symmetric, i.e., if $D(e^{j\omega}) = D^*(e^{-j\omega})$, then the resulting filter will be real-valued. In sum, the method described below can be used to design linear and non-linear phase filters with real-valued or complex-valued coefficients.

Define the complex desired frequency $D(e^{j\omega})$ on a dense grid of frequencies $\omega_i$, $i=1\ldots M$, in the interval $(-\pi,\pi]$. Let $H(e^{j\omega})$ be the frequency response of the FIR filter with coefficients $h[n]$: $$H(e^{j\omega})=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}$$ where $N$ is the filter length (i.e., the filter order equals $N-1$). The filter design problem can now simply be stated as $$H(e^{j\omega_i})\overset{!}{=}D(e^{j\omega_i}),\quad i=1\ldots M\tag{1}$$ where $\overset{!}{=}$ means "should be as close as possible" in some sense that is yet to be specified. We assume that $(1)$ is an overdetermined system, i.e., there are more equations than unknowns, or, in other words, there are more frequency points $M$ than unknown filter coefficients $N$.

It can be desirable to modify $(1)$ by introducing a strictly positive weight function $W(\omega)$. For notational convenience when formulating the least squares problem we use the square root of $W(\omega)$ in the overdetermined system: $$\sqrt{W}(\omega)H(e^{j\omega_i})\overset{!}{=}\sqrt{W}(\omega)D(e^{j\omega_i}),\quad i=1\ldots M\tag{2}$$ When minimizing the error in the overdetermined system $(2)$, the resulting approximation error will be smaller in frequency regions with a larger weight.

For formulating the weighted least squares approximation problem it is advantageous to write $(2)$ in matrix notation. Let $\mathbf C$ be a complex $M\times N$ matrix with rows ${\mathbf c}_i$ given by $${\mathbf c}_i=\big[1,e^{-j\omega_i},e^{-2j\omega_i},\ldots,e^{-(N-1)j\omega_i}\big]$$ With the filter coefficients $h[n]$ arranged in the column vector $\mathbf h$, the filter's frequency response evaluated at frequency $\omega_i$ can be written as $H(e^{j\omega_i})={\mathbf c}_i{\mathbf h}$. Consequently, the $M\times 1$ vector with elements $H(e^{j\omega_i})$, $i=1,\ldots,M$, is given by ${\mathbf C}{\mathbf h}$. With the complex desired frequency response values at frequencies $\omega_i$ arranged in a vector $\mathbf d$, we can write $(2)$ as $${\mathbf W}^{1/2}{\mathbf C \mathbf h} \overset{!}{=} {\mathbf W}^{1/2}{\mathbf d}\tag{3}$$ where $\mathbf W$ is a diagonal matrix with elements $W(\omega_i)$.

Solving $(3)$ in a least squares sense is equivalent to minimizing the following error measure: $$\begin{align}\epsilon&=\sum_{i=1}^M W(\omega_i)\big|H(e^{j\omega_i})-D(e^{j\omega_i})\big|^2\\&=({\mathbf C \mathbf h}-{\mathbf d})^H{\mathbf W}({\mathbf C \mathbf h}-{\mathbf d})\\&={\mathbf h}^H{\mathbf C}^H{\mathbf W}{\mathbf C}{\mathbf h}-{\mathbf h}^H{\mathbf C}^H{\mathbf W}{\mathbf d}-{\mathbf d}^H{\mathbf W}{\mathbf C}{\mathbf h}+{\mathbf d}^H{\mathbf d}\tag{4}\end{align}$$ where we use $^H$ to denote complex conjugate transposition.

Minimizing $(4)$ results in the following $N\times N$ system of linear equations for the unknown filter coefficients: $${\mathbf C}^H{\mathbf W}{\mathbf C}{\mathbf h} = {\mathbf C}^H{\mathbf W}{\mathbf d} \tag{5}$$

Now we can either solve $(5)$ using some linear equation solver, or we can directly use the overdetermined system $(3)$ as input to a QR-solver. Matlab and Octave offer both possibilities with the backslash `\` operator. Solving the overdetermined system directly using a QR-solver has numerical advantages, whereas solving the square system $(5)$ has the advantage that no matrix needs to be stored because the system matrix ${\mathbf C}^H{\mathbf W}{\mathbf C}$ has a Toeplitz structure, and is therefore completely determined by its first row. Such a system can be solved efficiently using Levinson's algorithm.

In general, the solutions to $(3)$ and $(5)$ are complex-valued. If we want to design real-valued filters, then for every frequency point $\omega_i>0$ there must be a point $-\omega_i$ for which $D(e^{-j\omega_i})=D^*(e^{j\omega_i})$ and $W(-\omega_i)=W(\omega_i)$ are satisfied. Another more efficient way to guarantee real-valued solutions is to slightly modify the systems $(3)$ and $(5)$. With these modifications we only need to specify frequency points in the interval $[0,\pi]$. The overdetermined system $(3)$ becomes $$\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf C}_R \\{\mathbf W}^{1/2}{\mathbf C}_I \end{array} \right]{\mathbf h}\overset{!}{=} \left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf d}_R\\{\mathbf W}^{1/2}{\mathbf d}_I\end{array}\right]\tag{6}$$ where subscripts $R$ and $I$ denote the real and imaginary parts, respectively. The $N\times N$ system $(5)$ becomes $$\textrm{Re}\left\{{\mathbf C}^H{\mathbf W}{\mathbf C}\right\}{\mathbf h} = \textrm{Re}\left\{{\mathbf C}^H{\mathbf W}{\mathbf d}\right\} \tag{7}$$ I.e., if we want a real-valued solution, we just take the real parts of the system matrix and of the right-hand side vector.

I've implemented the method described above in the Matlab/Octave function cfirls.ml:

Linear phase FIR filters: 4 types

This post is about the 4 types of linear phase finite impulse response (FIR) filters. The possibility to have a causal filter with an exactly linear phase response can be an advantage of FIR filters over infinite impulse response (IIR) filters for certain applications.

1. Type I: symmetrical impulse response, odd length
2. Type II: symmetrical impulse response, even length
3. Type III: asymmetrical impulse response, odd length
4. Type IV: asymmetrical impulse response, even length

When choosing one of these 4 types of linear phase filters there are mainly 3 things to consider:

1. constraints on the zeros of the transfer function $H(z)$ at $z=1$ (DC) and $z=-1$ (Nyquist)
2.  integer/non-integer group delay
3.  phase shift (apart from the linear phase)

For type I filters (odd number of taps, even symmetry) there are no constraints on the zeros at $z=1$ and $z=-1$, the phase shift is zero (apart from the linear phase), and the group delay is an integer value.

Type II filters (even number of taps, even symmetry) always have a zero at $z=-1$ (i.e., half the sampling frequency), they have a zero phase shift, and they have a non-integer group delay.

Type III filters (odd number of taps, odd symmetry) always have zeros at $z=1$ and $z=-1$ (i.e. at $f=0$ and $f=f_s/2$), they have a 90 degrees phase shift, and an integer group delay.

Type IV filters (even number of taps, odd symmetry) always have a zero at $z=1$, a phase shift of 90 degrees, and a non-integer group delay.

This implies (among other things) the following:

• Type I filters are pretty universal, but they cannot be used whenever a 90 degrees phase shift is necessary, e.g. for differentiators or Hilbert transformers.
• Type II filters would normally not be used for high pass or band stop filters, due to the zero at $z=-1$, i.e. at $f=f_s/2$. Neither can they be used for applications where a 90 degrees phase shift is necessary.
• Type III filters cannot be used for standard frequency selective filters because in these cases the 90 degrees phase shift is usually undesirable. For Hilbert transformers, type III filters have a relatively bad magnitude approximation at very low and very high frequencies due to the zeros at $z=1$ and $z=-1$. On the other hand, a type III Hilbert transformer can be implemented more efficiently than a type IV Hilbert transformer because in this case every other tap is zero.
• Type IV filters cannot be used for standard frequency selective filters, for the same reasons as type III filters. They are well suited for differentiators and Hilbert transformers, and their magnitude approximation is usually better because, unlike type III filters, they have no zero at $z=-1$.
• In some applications an integer group delay is desirable. In these cases type I or type III filters are preferred.

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.