Window Design by Linear Programming

This section, based on a class project by EE graduate student Tatsuki Kashitani, illustrates the use of linprog in Matlab for designing variations on the Chebyshev window3.10). In addition, some comparisons between standard linear programming and the Remez exchange algorithm (firpm) are noted.

Linear Programming (LP)

If we can get our filter or window design problems in the form

\begin{displaymath}\begin{array}[t]{ll} \mathrm{minimize} & f^{T}x\\ \mathrm{subject}\, \mathrm{to} & \begin{array}[t]{l} \mathbf{A}_{eq}x=b_{eq}\\ \mathbf{A}x\le b\end{array}, \end{array}\end{displaymath} (4.60)

where $ x,f\in{\bf R}^N$ , $ b\in{\bf R}^M$ , $ A$ is $ M\times N$ , etc., then we are done.

The linprog function in Matlab Optimization Toolbox solves LP problems. In Octave, one can use glpk instead (from the GNU GLPK library).


LP Formulation of Chebyshev Window Design

What we want:

  1. Symmetric zero-phase window.
  2. Window samples to be positive.

    $\displaystyle w\left(n\right)\geq 0\qquad \mathrm{for}\quad -\frac{M-1}{2}\leq n\leq \frac{M-1}{2}=L$ (4.61)

  3. Transform to be 1 at DC.

    $\displaystyle W\left(0\right)=1$ (4.62)

  4. Transform to be within $ \left[-\delta ,\delta \right]$ in the stop-band.

    $\displaystyle -\delta \leq W\left(\omega \right)\leq \delta \qquad \mathrm{for}\quad \omega _{sb}\leq \omega \leq \pi$ (4.63)

  5. And $ \delta $ to be small.


Symmetric Window Constraint

Because we are designing a zero-phase window, use only the positive-time part $ h\left(n\right)$ :

$\displaystyle h\left(n\right)=w\left(n\right)\qquad n\geq 0$ (4.64)

$\displaystyle w\left(n\right)=\left\{ \begin{array}{ll} h\left(n\right) & n\geq 0\\ h\left(-n\right) & n<0\end{array} \right.$ (4.65)


Positive Window-Sample Constraint

For each window sample, $ h\left(n\right) \geq 0$ , or,

$\displaystyle -h\left(n\right) \leq 0.$ (4.66)

Stacking inequalities for all $ n$ ,

$\displaystyle \left[\begin{array}{ccccc} -1 & 0 & \cdots & 0 & 0\\ 0 & -1 & & & 0\\ \vdots & & \ddots & & \vdots \\ 0 & & & -1 & 0\\ 0 & 0 & \cdots & 0 & -1\end{array} \right]\left[\begin{array}{c} h\left(0\right)\\ h\left(1\right)\\ \vdots \\ h\left(L-1\right)\\ h\left(L\right)\end{array} \right] \le \left[\begin{array}{c} 0\\ 0\\ \vdots \\ 0\\ 0 \end{array} \right]$ (4.67)

or

$\displaystyle \zbox {-\mathbf{I}\, h \le 0.}$ (4.68)


DC Constraint

The DTFT at frequency $ \omega$ is given by

$\displaystyle W\left(\omega \right)=\sum _{n=-L}^{L}w\left(n\right)e^{-j\omega n}.$ (4.69)

By zero-phase symmetry,
$\displaystyle W\left(\omega \right)$ $\displaystyle =$ $\displaystyle h\left(0\right)+2\sum _{n=1}^{L}h\left(n\right)\cos \left(n\omega \right)$  
  $\displaystyle =$ $\displaystyle \left[\begin{array}{cccc}
1 & 2\cos \left(\omega \right) & \cdots & 2\cos \left(L\omega \right)\end{array}\right]\left[\begin{array}{c}
h\left(0\right)\\
h\left(1\right)\\
\vdots \\
h\left(L\right)\end{array}\right]$  
  $\displaystyle =$ $\displaystyle d\left(\omega \right)^{T}h.$  

So $ W\left(0\right)=1$ can be expressed as

$\displaystyle \zbox {d\left(0\right)^{T}h=1.}$ (4.70)


Sidelobe Specification

Likewise, side-lobe specification can be enforced at frequencies $ \omega_{i}$ in the stop-band.

$\displaystyle -\delta \leq d\left(\omega _{i}\right)^{T}h\leq \delta$ (4.71)

or

$\displaystyle \left\{ \begin{array}{ccc} -d\left(\omega _{i}\right)^{T}h-\delta & \leq & 0\\ \;d\left(\omega _{i}\right)^{T}h-\delta & \leq & 0\end{array} \right.$ (4.72)

where

$\displaystyle \omega _{sb}\leq \omega _1,\omega _{2}\ldots ,\omega _{K}\leq \pi .$ (4.73)

We need $ K\gg L$ to obtain many frequency samples per side lobe. Stacking inequalities for all $ \omega_{i}$ ,
$\displaystyle \left[\begin{array}{c}
-d\left(\omega _1\right)^{T}\\
\vdots \\
-d\left(\omega _{K}\right)^{T}\\
d\left(\omega _1\right)^{T}\\
\vdots \\
d\left(\omega _{K}\right)^{T}\end{array}\right]h+\left[\begin{array}{c}
-\delta \\
\vdots \\
-\delta \\
-\delta \\
\vdots \\
-\delta \end{array}\right]$ $\displaystyle \le$ $\displaystyle \mathbf{0}$  
$\displaystyle \left[\begin{array}{cc}
-d\left(\omega _1\right)^{T} & -1\\
\vdots & \vdots \\
-d\left(\omega _{K}\right)^{T} & -1\\
d\left(\omega _1\right)^{T} & -1\\
\vdots & \vdots \\
d\left(\omega _{K}\right)^{T} & -1
\end{array}\right]\left[\begin{array}{c}
h\\
\delta \end{array}\right]$ $\displaystyle \le$ $\displaystyle \mathbf{0}.$  

I.e.,

$\displaystyle \zbox {\mathbf{A}_{sb}\left[\begin{array}{c} h\\ \delta \end{array} \right] \le \mathbf{0}.}$ (4.74)


LP Standard Form

Now gather all of the constraints to form an LP problem:

\begin{displaymath}\begin{array}[t]{ll} \mathrm{minimize} & \left[\begin{array}{cccc} 0 & \cdots & 0 & 1\end{array} \right] \left[\begin{array}{c} h\\ \delta \end{array} \right]\\ [5pt] \mbox{subject to} & \begin{array}[t]{l} \left[\begin{array}{cc} d\left(0\right)^{T} & 0\end{array} \right]\left[\begin{array}{c} h\\ \delta \end{array} \right]=1\\ \left[\begin{array}{c} \left[\begin{array}{cc} -\mathbf{I} & \mathbf{0}\end{array} \right]\\ [5pt] \mathbf{A}_{sb}\end{array} \right]\left[\begin{array}{c} h\\ \delta \end{array} \right]\le \mathbf{0}\end{array} \end{array}\end{displaymath} (4.75)

where the optimization variables are $ [h, \delta]^T$ .

Solving this linear-programming problem should produce a window that is optimal in the Chebyshev sense over the chosen frequency samples, as shown in Fig.3.37. If the chosen frequency samples happen to include all of the extremal frequencies (frequencies of maximum error in the DTFT of the window), then the unique Chebyshev window for the specified main-lobe width must be obtained. Iterating to find the extremal frequencies is the heart of the Remez multiple exchange algorithm, discussed in the next section.

Figure 3.37: Normal Chebyshev Window
\includegraphics[width=\twidth,height=6.5in]{eps/print_normal_chebwin}


Remez Exchange Algorithm

The Remez multiple exchange algorithm works by moving the frequency samples each iteration to points of maximum error (on a denser grid). Remez iterations could be added to our formulation as well. The Remez multiple exchange algorithm (function firpm [formerly remez] in the Matlab Signal Processing Toolbox, and still remez in Octave) is normally faster than a linear programming formulation, which can be regarded as a single exchange method [224, p. 140]. Another reason for the speed of firpm is that it solves the following equations non-iteratively for the filter exhibiting the desired error alternation over the current set of extremal frequencies:

$\displaystyle \left[ \begin{array}{c} H(\omega_1) \\ H(\omega_2) \\ \vdots \\ H(\omega_{K}) \end{array} \right] = \left[ \begin{array}{cccccc} 1 & 2\cos(\omega_1) & \dots & 2\cos(\omega_1L) & \frac{1}{W(\omega_1)} \\ 1 & 2\cos(\omega_2) & \dots & 2\cos(\omega_2L) & \frac{-1}{W(\omega_2)} \\ \vdots & & & \\ 1 & 2\cos(\omega_{K}) & \dots & 2\cos(\omega_{K}L) & \frac{(-1)^{K}}{W(\omega_{K})} \end{array} \right] \left[ \begin{array}{c} h_0 \\ h_1 \\ \vdots \\ h_{L} \\ \delta \end{array} \right]$ (4.76)

where $ W(\omega_k)\delta$ is the weighted ripple amplitude at frequency $ \omega_k$ . ( $ W(\omega_k)$ is an arbitrary ripple weighting function.) Note that the desired frequency-response amplitude $ H(\omega_k)$ is also arbitrary at each frequency sample.

Convergence of Remez Exchange

According to a theorem of Remez, $ \delta $ is guaranteed to increase monotonically each iteration, ultimately converging to its optimal value. This value is reached when all the extremal frequencies are found. In practice, numerical round-off error may cause $ \delta $ not to increase monotonically. When this is detected, the algorithm normally halts and reports a failure to converge. Convergence failure is common in practice for FIR filters having more than 300 or so taps and stringent design specifications (such as very narrow pass-bands). Further details on Remez exchange are given in [224, p. 136].

As a result of the non-iterative internal LP solution on each iteration, firpm cannot be used when additional constraints are added, such as those to be discussed in the following sections. In such cases, a more general LP solver such as linprog must be used. Recent advances in convex optimization enable faster solution of much larger problems [22].


Monotonicity Constraint

We can constrain the positive-time part of the window to be monotonically decreasing:

$\displaystyle \Delta h_{i}=h\left(i+1\right)-h\left(i\right)\leq 0\qquad i=1,\ldots ,L-1$ (4.77)

In matrix form,
$\displaystyle \left[\begin{array}{ccccc}
-1 & 1 & & & 0\\
& -1 & 1 & & \\
& & \ddots & \ddots & \\
0 & & & -1 & 1\end{array}\right]h$ $\displaystyle \le$ $\displaystyle 0,$  

or,

$\displaystyle \zbox {\mathbf{D}\,h \le 0.}$ (4.78)

See Fig.3.38.

Figure 3.38: Monotonic Chebyshev Window
\includegraphics[width=\twidth,height=6.5in]{eps/print_monotonic_chebwin}


L-Infinity Norm of Derivative Objective

We can add a smoothness objective by adding $ L-infinity$ -norm of the derivative to the objective function.

$\displaystyle \mathrm{minimize}\quad \delta +\eta \left\Vert \Delta h\right\Vert _{\infty }.$ (4.79)

The $ L-infinity$ -norm only cares about the maximum derivative. Large $ \eta $ means we put more weight on the smoothness than the side-lobe level.

This can be formulated as an LP by adding one optimization parameter $ \sigma $ which bounds all derivatives.

$\displaystyle -\sigma \leq \Delta h_{i}\leq \sigma \qquad i=1,\ldots ,L-1.$ (4.80)

In matrix form,
$\displaystyle \left[\begin{array}{r}
-\mathbf{D}\\
\mathbf{D}\end{array}\right]h-\sigma \mathbf1$ $\displaystyle \le$ $\displaystyle 0.$  

Objective function becomes

$\displaystyle \mathrm{minimize}\quad \delta +\eta \sigma .$ (4.81)

The result of adding the Chebyshev norm of diff(h) to the objective function to be minimized ($ \eta =1$ ) is shown in Fig.3.39. The result of increasing $ \eta $ to 20 is shown in Fig.3.40.

Figure: Chebyshev norm of diff(h) added to the objective function to be minimized ($ \eta =1$ )
\includegraphics[width=\twidth,height=6.5in]{eps/print_linf_chebwin_1}

Figure: Twenty times the norm of diff(h) added to the objective function to be minimized ($ \eta =20$ )
\includegraphics[width=\twidth,height=6.5in]{eps/print_linf_chebwin_2}


L-One Norm of Derivative Objective

Another way to add smoothness constraint is to add $ L1$ -norm of the derivative to the objective:

$\displaystyle \mathrm{minimize}\quad \delta +\eta \left\Vert \Delta h\right\Vert _1$ (4.82)

Note that the $ L1$ norm is sensitive to all the derivatives, not just the largest.

We can formulate an LP problem by adding a vector of optimization parameters $ \tau$ which bound derivatives:

$\displaystyle -\tau _{i}\leq \Delta h_{i}\leq \tau _{i}\qquad i=1,\ldots ,L-1.$ (4.83)

In matrix form,

$\displaystyle \left[\begin{array}{r} -\mathbf{D}\\ \mathbf{D}\end{array} \right]h-\left[\begin{array}{c} -\tau \\ -\tau \end{array} \right]\le 0.$ (4.84)

The objective function becomes

$\displaystyle \mathrm{minimize}\quad \delta +\eta \mathbf1^{T}\tau .$ (4.85)

See Fig.3.41 and Fig.3.42 for example results.

Figure: $ L1$ norm of diff(h) added to the objective function ($ \eta =1$ )
\includegraphics[width=\twidth,height=6.5in]{eps/print_lone_chebwin_1}

Figure: Six times the $ L1$ norm of diff(h) added to the objective function ($ \eta =6$ )
\includegraphics[width=\twidth,height=6.5in]{eps/print_lone_chebwin_2}


Summary

This section illustrated the design of optimal spectrum-analysis windows made using linear-programming (linprog in matlab) or Remez multiple exchange algorithms (firpm in Matlab). After formulating the Chebyshev window as a linear programming problem, we found we could impose a monotonicity constraint on its shape in the time domain, or various derivative constraints. In Chapter 4, we will look at methods for FIR filter design, including the window method4.5) which designs FIR filters as a windowed ideal impulse response. The formulation introduced in this section can be used to design such windows, and it can be used to design optimal FIR filters. In such cases, the impulse response is designed directly (as the window was here) to minimize an error criterion subject to various equality and inequality constraints, as discussed above for window design.4.16


Next Section:
The Ideal Lowpass Filter
Previous Section:
Optimized Windows