4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Just curious, did you try Matt L's suggestion in the second answer? It appears he worked this all out as his PhD thesis, including Matlab code for FIR and IIR implementations (it looks like ccirls on page 108 will do what you want?) drive.google.com/file/d/0B_s50T2gL6HaUWZwazJvT3BZTmc/view Nice work @MattL! $\endgroup$ – Dan Boschen Nov 12 '18 at 12:27
  • 1
    $\begingroup$ @DanBoschen: Thanks Dan, looks like I have to read up on this in my own thesis after all those years ... :) $\endgroup$ – Matt L. Nov 12 '18 at 13:02
  • $\begingroup$ @MattL I know the feeling! That is why I am glad I write everything down but nothing I have has approached the detail of that thesis. The appendices are especially helpful also- well done! $\endgroup$ – Dan Boschen Nov 12 '18 at 13:20
  • $\begingroup$ Yes, I'm using the lslevin functions with great success (thank you a lot for making those available to everyone!). I also tried the LS IIR functions with mixed success but I may just have used the wrong metrics for comparison so far, so I can't really tell how good they work in practice. $\endgroup$ – TheSlowGrowth Nov 13 '18 at 8:34
  • $\begingroup$ i think Eric Jacobsen (a comp.dsp regular, but i don't see him much around here) has advocated doing the very same thing. $\endgroup$ – robert bristow-johnson Nov 13 '18 at 18:52
4
$\begingroup$

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.

enter image description here

$\endgroup$
  • $\begingroup$ Interesting, thank you a lot. I actually worked it out to almost to the same form but didn't consider that the additional linear phase term $e^{-j\omega(N-1)/2}$ would indeed remain in the final transfer function. Can you point me towards the reason for the additional 3 dB worst-case error? $\endgroup$ – TheSlowGrowth Nov 13 '18 at 8:48
  • 1
    $\begingroup$ I've updated the answer to include that information and a design example. $\endgroup$ – Matt L. Nov 13 '18 at 18:56
  • $\begingroup$ Thank you a lot for the time and effort, this is extremely helpful! $\endgroup$ – TheSlowGrowth Nov 14 '18 at 8:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.