FIR design for arbitrary magnitude and phase: Why can I separate real and imaginary parts like this?

Question

In this question about the design of a FIR filter with arbitrary magnitude and phase specifications user robert bristow-johnson suggested to split the desired complex frequency response $H$ into its real and imaginary parts $H_R$ and $H_I$, then use firpm or firls to design two linear phase filters for both, setting the ftypeoption to hilbert for the imaginary part. Apparently, the sum of the two resulting FIR filters would yield a filter that meets the specification.

I'm trying to understand why this works. Here's my attempt: $$H_R = |H|\cos(\arg\{H\})$$ $$H_I = |H|\sin(\arg\{H\})$$ Now we design two FIR filters with firpm or firls, setting the ftypeoption to hilbert for the imaginary part. To my understanding, the frequency responses of the two resulting filters $\tilde h_R$ and $\tilde h_I$ will then look like this: $$\tilde H_R = \Re\{H\}\ e^{jp(k)}$$ $$\tilde H_I = \Im\{H\}\ e^{jp(k)+j\frac\pi2}$$ Where $p(k)$ is the linear phase term and the second filter has the additional phase shift due to the hilbert option that makes it a type III or IV FIR filter. The sum of the coefficients $\tilde h[k] = \tilde h_R[k] + \tilde h_I[k]$ would then be a FIR with the frequency response $$\tilde H = \tilde H_R + \tilde H_I$$ because adding the coefficients is basically the same as adding the filter output signals for the same filter input signal.

\begin{align}\tilde H =& \Re\{H\}\cos(p(\omega)) + \Im\{H\}\cos(p(\omega) + \frac\pi2) \\&+ j\left[\Re\{H\}\sin(p(\omega)) + \Im\{H\}\sin(p(\omega) + \frac\pi2)\right]\\ =& \Re\{H\}\cos(p(\omega)) - \Im\{H\}\sin(p(\omega)) \\&+ j\left[\Re\{H\}\sin(p(\omega)) - \Im\{H\}\cos(p(\omega))\right] \end{align}

This looks nothing like the response $H$ that I actually want. And I don't see how this could work. Is the original approach wrong or am I?

Answer 1

Neither of you is wrong. Your final expression is just the same as $H(\omega)\cdot e^{-j\omega(N-1)/2}$, which is (an approximation of) the desired complex frequency response with an additional linear phase term that depends on the chosen filter length $N$. This means that you need to subtract that linear phase term from the desired phase response, because it will be added in the design process.

Let the desired complex frequency response be given by

$$D(\omega)=|D(\omega)|\cdot e^{j\phi(\omega)}\tag{1}$$

Now you have to compute a new complex response by subtracting the linear phase term $-\omega (N-1)/2$ from the phase of $(1)$:

$$\tilde{D}(\omega)=|D(\omega)|\cdot e^{j[\phi(\omega)+\omega(N-1)/2]}\tag{2}$$

You have to use the real and imaginary parts of $\tilde{D}(\omega)$ to design the filter.

What is meant by the option 'hilbert', is that you tell the design routine that the desired response (the imaginary part of $\tilde{D}(\omega)$) should actually be multiplied by $j$ and it is an odd function of frequency, because otherwise it will be interpreted as a purely real and even function. This will result in an asymmetric impulse response. In principle this can be done, but I don't know if it can actually be done for arbitrary prescribed responses with the current versions of Matlab's filter design routines.

Note that such an approach is actually not really necessary for the least squares design of FIR filters, because solving a complex least squares problem (for non-linear phase designs) is just as simple as solving a real least squares problem (for linear phase designs). As far as I know it's just not implemented in Matlab's Signal Processing Toolbox. I've written a few Matlab/Octave functions that solve the complex least squares approximation problem for FIR filter design: lslevin.m for complex least squares design of real-valued filters, and cfirls.m for filters with complex coefficients.

For equi-ripple designs, the complex (non-linear phase) approximation problem is indeed much harder than the real (linear phase) problem. For the latter we have the Remez exchange algorithm, which was adapted to the linear phase FIR filter design problem by Parks and McClellan and which is implemented in Matlab (firpm.m). There are algorithms for solving the complex equi-ripple design problem, but they are much less efficient than the Parks-McClellan algorithm. One of them is implemented in Matlab (cfirpm.m).

Note that for the least squares approximation, solving two real-valued approximation problems for the real and imaginary parts gives the same solution as solving one complex least squares approximation problem. This is easily seen as follows. The complex approximation error is given by

$$E(\omega)=D(\omega)-H(\omega)\tag{3}$$

where $D(\omega)$ and $H(\omega)$ are the desired and actual complex frequency responses. The complex least squares solution minimizes the integral over $|E(\omega)|^2$, whereas the individual real-valued least squares solutions minimize the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$, where $E_R(\omega)$ and $E_I(\omega)$ are the real and imaginary parts of $E(\omega)$, respectively. Since

$$|E(\omega)|^2=E_R^2(\omega)+E_I^2(\omega)\tag{4}$$

minimizing the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$ is equivalent to minimizing the integral over $|E(\omega)|^2$.

For the complex Chebyshev approximation problem the situation is different. The quantity that is minimized is the maximum of the magnitude of the complex error over the bands of interest $\max_{\omega}|E(\omega)|$, but since

$$\max_{\omega}|E(\omega)|\neq \max_{\omega}|E_R(\omega)|+\max_{\omega}|E_I(\omega)|\tag{5}$$

independent Chebyshev approximation of the real and imaginary parts is not equivalent to complex Chebyshev approximation. If for a given complex approximation problem the maximum error of the optimal solution equals $\delta$, then in the worst case the errors of both real-valued approximations will also be equal to $\delta$. This means that in the worst case, the maximum complex error of the solution obtained from independent optimization of the real and imaginary parts becomes

$$\max_{\omega}|E(\omega)|=\max_{\omega}\sqrt{E_R^2(\omega)+E_I^2(\omega)}\le\sqrt{\max_{\omega}E_R^2(\omega)+\max_{\omega}E_I^2(\omega)}\le\sqrt{2}\delta\tag{6}$$

This is equivalent to a maximum of $3$ dB increase in error compared to the optimal solution.

I'll show an example design to illustrate a few points that I made. I designed a bandpass filter with a desired phase response that is linear but with a smaller group delay than an FIR filter with perfectly linear phase. The filter length is $N=61$, so the delay of a linear phase FIR filter would be $(N-1)/2=30$ samples. I chose the desired delay to be $20$ samples instead. I used a complex least squares and a complex Chebyshev criterion to approximate the given complex frequency response. For each criterion I designed the filter in two different ways. First, directly in the complex domain, and second by using the real and imaginary parts of the desired frequency response to design two linear phase filters, which, after combining them, should approximate the original complex frequency response. The figure below shows the design results. As expected, for the least squares design, both solutions are equivalent (up to numerical accuracy) and only one of them is shown in the plots (in blue). For the Chebyshev criterion both solutions are slightly different, as predicted. The solution from splitting the response into real and imaginary parts (in red) has a larger error by about $1.2$dB (as can be seen at the first sidelobe in the stopband). In the passband all three solutions are very similar. Note that all three filters approximate the desired group delay value of $20$ samples relatively well over a large portion of the passband.