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:

	function h = cfirls(N,f,d,w,hreal)
	% function h = cfirls(N,f,d,w,hreal)
	% Least-squares FIR filter design with prescribed magnitude and phase
	% responses and with possibly complex coefficients
	%
	% h real or complex-valued impulse response
	% N desired filter length
	% f frequency grid
	% for complex-valued filters ('hreal=0'): -1 <= f < 1
	% ( '1' corresponds to Nyquist) or 0 <= f < 2
	% for real-valued filters ('hreal=1'): 0 <= f <= 1
	% d desired complex-valued frequency response on the grid f:
	% d = ( desired_magnitude ) .* exp( 1i * desired_phase )
	% w (positive) weighting function on the grid f
	% hreal flag for constraining the filter coefficients h to be real-valued
	% hreal=1: h is real-valued
	% hreal=0: h can be complex-valued
	% h will be real-valued (up to numerical accuracy) if for
	% every grid point f[i] there is a point -f[i] and the
	% desired response at -f[i] is the complex conjugate of d
	% at f[i], and the weights at f[i] and -f[i] must be equal
	% differences with firls.m:
	% 1. the desired phase response can be non-linear
	% 2. the filter coefficients can be complex (i.e., the specified
	% frequency response does not need to be symmetrical)
	% 3. the desired magnitude doesn't need to be piecewise linear,
	% but it is specified as a vector defined on a frequency grid
	% 4. the weighting does not need to be piecewise constant
	% EXAMPLE 1: real-valued linear phase bandpass filter
	% N = 61;
	% tau = (N-1)/2;
	% f = [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
	% d = [zeros(1,230),ones(1,200),zeros(1,430)] .* exp(-j*tau*pi*f);
	% w = [10*ones(1,230),ones(1,200),10*ones(1,430)];
	% h = cfirls(N,f,d,w,1);
	% EXAMPLE 2: same as Ex. 1 but with reduced delay in the passband
	% N = 61;
	% tau = 25;
	% f = [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
	% d = [zeros(1,230),ones(1,200),zeros(1,430)] .* exp(-j*tau*pi*f);
	% w = [10*ones(1,230),ones(1,200),10*ones(1,430)];
	% h = cfirls(N,f,d,w,1);
	% EXAMPLE 3: linear phase complex low pass filter
	% N = 51;
	% tau = (N-1)/2;
	% f = [linspace(-1,-.18,328),linspace(-.1,.3,160),linspace(.38,1,248)];
	% d = [zeros(1,328),ones(1,160),zeros(1,248)].*exp(-1i*pi*f*tau);
	% w = [10*ones(1,328),ones(1,160),10*ones(1,248)];
	% h = cfirls(N,f,d,w);
	% EXAMPLE 4: same as Ex. 3 but with reduced delay in the passband
	% N = 51;
	% tau = 20;
	% f = [linspace(-1,-.18,328),linspace(-.1,.3,160),linspace(.38,1,248)];
	% d = [zeros(1,328),ones(1,160),zeros(1,248)].*exp(-1i*pi*f*tau);
	% w = [10*ones(1,328),ones(1,160),10*ones(1,248)];
	% h = cfirls(N,f,d,w);
	%
	% author: Mathias Lang, 2018-01-01
	% revision 2022-10-17: added option 'hreal' for designing exactly
	% real-valued filters
	% mattsdspblog@gmail.com
	% check arguments
	if ( nargin > 5 | nargin < 4 ),
	error('wrong number of input arguments');
	end
	if ( nargin == 4 ), hreal = 0; end
	L = length(f);
	if ( length(d) ~= L ), error('d and f must have the same lengths.'); end
	if ( length(w) ~= L ), error('w and f must have the same lengths.'); end
	if ( N < 1 ), error('N must be greater than 0.'); end
	if ~all( w ),
	error('All elements of the weight vector w must be positive.');
	end
	minf = min(f); maxf = max(f);
	if ( hreal ),
	if ( minf < 0 | maxf > 1 ),
	error('f must be in the range [0,1] for real-valued filters.');
	end
	else
	if ~( ( minf >= -1 & maxf <= 1 ) | ( minf >= 0 & maxf <= 2 ) ),
	error('f must be in the interval [-1,1) or [0,2).');
	end
	end
	f = f(:);
	d = d(:);
	w = sqrt(w(:));
	% solve overdetermined system in a least squares sense
	A = w( :, ones(1,N) ) .* exp( -1i * pi * f * (0:N-1) );
	b = w .* d;
	if ( hreal )
	h = [real(A); imag(A)] \ [real(b); imag(b)];
	else
	h = A \ b;
	end

view raw cfirls.m hosted with ❤ by GitHub