Force Linear Phase for a FIR Filter Synthesized Using Berchin's FDLS?

Question

As a follow-on to this post, how would you force linear phase for a FIR filter you synthesized using Berchin's FDLS?

Answer 1

Remember that, as its name suggests, Berchin's frequency-domain least squares method (FDLS) is really just a least-squares fit of a desired frequency-domain response. Recall that for an arbitrary transfer function of the form:

$$ H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{\sum_{k=1}^{N} a_k z^{-k}} $$

you can write the system's frequency response by simply letting $z=e^{j\omega}$:

$$ H(\omega) = \frac{\sum_{k=0}^{M} b_k e^{-j\omega k}}{\sum_{k=1}^{N} a_k e^{-j\omega k}} $$

This sort of system can be implemented as an IIR filter with $b_k$ as the feedforward (numerator) coefficients and $a_k$ as the feedback (denominator) coefficients. The FDLS method simply requires you to define the desired (complex-valued) frequency response at a collection of frequencies $\{\omega_k\}$. Evaluating the above expression for each frequency $\omega_k$ yields a system of linear equations in the $b_k$ and $a_k$ variables. You can then solve this using least-squares methods to yield the best coefficients for the desired frequency response and filter order (in the LS sense, of course).

If you have the requirement for exact linear phase (i.e. better linear phase response than you might get by just applying the method itself and specifying a desired linear phase response), then you could reformulate the above a bit, if you're willing to accept an FIR-only filter structure.

• Recall that linear phase in the frequency domain corresponds to symmetry in the filter's time-domain impulse response. For FIR filters, this impulse response is the same as the filter's coefficients $b_k$.

• Let $a_k=0\ \forall\ k$ in the above expression, yielding: $$ H(\omega) = \sum_{k=0}^{M} b_k e^{-j\omega k} $$

• The symmetry constraint corresponds to a loss in the number of degrees of freedom in the linear system; instead of each of the coefficients being independent variables, the linear-phase constraint requires that coefficients that are symmetric about the impulse response's midpoint to be equal.
• Rewrite the above equation taking these constraints into account. For example, for a type I FIR filter (odd length, even symmetry), then you would have something like: $$ H(\omega) = \sum_{k=0}^{\frac{M}{2}-1} b_k e^{-j\omega k} + \sum_{k=\frac{M}{2}+1}^{M} b_{M-k} e^{-j\omega k} + b_{\frac{M}{2}} e^{-j\omega \frac{M}{2}} $$ where the filter has $M+1$ taps. You would then solve the system for the reduced set of coefficients, then reconstruct the overall filter impulse response using the symmetry constraint: $$ h[k] = \begin{cases} b_k, && 0 \le k \le \frac{M}{2} \\ b_{M-k}, && \frac{M}{2}+1 \le k \le M \end{cases} $$ The result is an FIR filter that will have linear phase.

Answer 2

Jason R's answer is correct, but a little more information is needed in order to know how to apply it specifically to FDLS. That information may be found here:

https://groups.google.com/forum/?hl=en#!original/comp.dsp/qgYsnf5aIXs/NSg6VSXhkqoJ

If you are using my MATLAB code for FDLS, then you will have to modify it as indicated in that link.

Greg Berchin

Answer 3

Given an FDLS solution as a transfer function with a closest fit to your desired magnitude response, all possible stable phase responses can be derived by trying all permutations of flipping the zeros of that transfer function inside and outside the unit circle. This set can be evaluated for the closest fit to your desired phase response, linear or otherwise. The closest fit in phase for that constellation of poles and zeros may not be an exact fit.