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:
No comments:
Post a Comment