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

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

Answer

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.