The design problem is to find the coefficients $a_n$ in $(2)$ such that the phase $\phi(\omega)$ approximates a given desired phase response $\phi_d(\omega)$. The filter's phase $\phi(\omega)$ is a highly non-linear function of the filter coefficients $a_n$, so it's difficult to directly minimize some norm of the phase error $\phi(\omega)-\phi_d(\omega)$. It is more straightforward to consider the complex error function $H(e^{j\omega})-e^{j\phi_d(\omega)}$. This error function still depends non-linearly on the filter coefficients. However, multiplying this error function by the denominator of $H(e^{j\omega})$ results in an error function which is a linear function of the filter coefficients: $$\epsilon(\omega)=e^{-jN\omega}A^*(e^{j\omega}) - A(e^{j\omega})e^{j\phi_d(\omega)}\tag{4}$$ This error function is usually called "equation error", and it is also used for designing general IIR filters. In this DSP stackexchange post I describe the basics of that method for the general case, and there are also some relevant links.
Minimization of the error function $(4)$ is equivalent to solving the following overdetermined system of equations: $$e^{-jN\omega_i}A^*(e^{j\omega_i})\overset{!}{=}A(e^{j\omega_i})e^{j\phi_d(\omega_i)},\quad i=1,\ldots,M\tag{5}$$ where we've chosen a dense frequency grid $\omega_i$ in the interval $[0,\pi]$ such that there are considerably more equations than unknown filter coefficients. Since the equations in $(5)$ are linear functions of the filter coefficients $a_n$, we can rewrite $(5)$ using matrix/vector notation: $${\mathbf C}{\mathbf a}\overset{!}{=}{\mathbf d}\tag{6}$$ where $\mathbf C$ is a complex valued $M\times N$ matrix with rows $${\mathbf c}_i=\left[e^{-j(N-1)\omega_i}-e^{j(\phi_d(\omega_i)-\omega_i)},e^{-j(N-2)\omega_i}-e^{j(\phi_d(\omega_i)-2\omega_i)},\ldots,1-e^{j(\phi_d(\omega_i)-N\omega_i)}\right]\tag{7}$$ $\mathbf a$ is the $N\times 1$ vector of unknown filter coefficients $a_n$, and $\mathbf d$ is a complex $M\times 1$ vector with elements $e^{j\phi_d(\omega_i)}-e^{-jN\omega_i}$.
We can solve $(6)$ in a least squares sense either by solving the corresponding normal equations or we can directly use the overdetermined system $(6)$ as the input to a QR-solver, as provided by the backslash operator in Matlab or Octave. We could compute the filter coefficients from $(6)$ but we'll add a minor generalization in terms of a weight function which allows to weigh the error differently in different frequency regions. We define a positive weight function $W(\omega)$ on the frequency grid and we define an $M\times M$ diagonal matrix $\mathbf W$ with elements $W(\omega_i)$ on its diagonal. The overdetermined system we are going to solve is then $${\mathbf W}^{1/2}{\mathbf C}{\mathbf a}\overset{!}{=}{\mathbf W}^{1/2}{\mathbf d}\tag{8}$$ We use the square root of the weights in $(8)$ just to be consistent with other weighted least squares methods. Since we're only interested in real-valued solutions, we rewrite $(8)$ as $$\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf C}_R\\{\mathbf W}^{1/2}{\mathbf C}_I\end{array}\right]{\mathbf a}\overset{!}{=}\left[\begin{array}{c}{\mathbf W}^{1/2}{\mathbf d}_R\\{\mathbf W}^{1/2}{\mathbf d}_I\end{array}\right]\tag{9}$$ where the subscripts $R$ and $I$ denote real and imaginary parts, respectively. The overdetermined system of linear equations $(9)$ is the one that we need to solve in a least squares sense in order to compute the filter coefficients that minimize the weighted equation error.
It should be noted that the method described here does not guarantee the stability of the designed filter. If the resulting filter turns out to be unstable, reflecting poles about the unit circle is not an option, because that would ruin the phase approximation. Instead, a new design has to be done whereby the desired phase response is modified by a linear phase term, i.e., by adding (or subtracting) a constant delay. In this way, it is always possible to obtain a stable filter. The value of that delay, however, needs to be determined experimentally.
I'll show a few design examples to illustrate the method described above. The first design is a fractional delay. The goal is to design an allpass filter with a fractional delay of $7.5$ samples. We choose a filter order $N=8$ and no weighting, i.e., $W(\omega)=1$ for all frequencies. The figures below shows the result of this design. The top figure shows the phase error in radians. The phase error is smaller than $0.05$ radians for more than $90$% of the bandwidth. Only close to the Nyquist frequency does the phase error become relatively large. The bottom figure shows the group delay, which approximates the desired $7.5$ samples quite well for about $85$% of the bandwidth.
In the next design example, the allpass filter is used to equalize the phase distortion of a digital Cauer lowpass filter. The desired phase response of the allpass filter is the negative phase of the Cauer filter with some linear phase subtracted, such that the total desired phase response can be approximated well by a causal filter. I used a $5^{th}$ order Cauer lowpass and a $12^{th}$ order allpass equalizer. A good result was achieved with a desired total delay of $10$ samples. The figure below shows the design result. The top figure shows the phase approximation error in the passband, and the bottom plot shows the group delay functions in the passband for the Cauer filter (blue), the allpass equalizer (orange), and the cascade of Cauer filter and allpass filter (yellow). The total group delay is close to the specified $10$ samples over a large range of the passband.
The final example is the design of a Hilbert transformer. A Hilbert transformer shifts the phase of the input signal by $-\pi/2$ while keeping the amplitude unchanged. An ideal Hilbert transformer is a non-causal system, that's why a delay has to be added to the desired phase such that the approximation error of a causal allpass becomes reasonably small. I chose the filter order $N=10$ and a delay of $9$ samples. The top figure below shows the phase error. As is to be expected, the approximation is more difficult close to DC and Nyquist. However, over a large range of the bandwidth, the resulting filter computes an accurate Hilbert transform. The bottom figure shows the group delay, which approximates the desired value of $9$ samples, with larger errors at low and high frequencies.
No comments:
Post a Comment