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 window (§3.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
(4.60) |
where , , is , 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:
- Symmetric zero-phase window.
- Window samples to be positive.
(4.61)
- Transform to be 1 at DC.
(4.62)
- Transform to be within
in the stop-band.
(4.63)
- And to be small.
Symmetric Window Constraint
Because we are designing a zero-phase window, use only the positive-time part :
(4.64) |
(4.65) |
Positive Window-Sample Constraint
For each window sample, , or,
(4.66) |
Stacking inequalities for all ,
(4.67) |
or
(4.68) |
DC Constraint
The DTFT at frequency is given by
(4.69) |
By zero-phase symmetry,
So can be expressed as
(4.70) |
Sidelobe Specification
Likewise, side-lobe specification can be enforced at frequencies in the stop-band.
(4.71) |
or
(4.72) |
where
(4.73) |
We need to obtain many frequency samples per side lobe. Stacking inequalities for all ,
I.e.,
(4.74) |
LP Standard Form
Now gather all of the constraints to form an LP problem:
(4.75) |
where the optimization variables are .
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.
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:
(4.76) |
where is the weighted ripple amplitude at frequency . ( is an arbitrary ripple weighting function.) Note that the desired frequency-response amplitude is also arbitrary at each frequency sample.
Convergence of Remez Exchange
According to a theorem of Remez, 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 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:
(4.77) |
In matrix form,
or,
(4.78) |
See Fig.3.38.
L-Infinity Norm of Derivative Objective
We can add a smoothness objective by adding -norm of the derivative to the objective function.
(4.79) |
The -norm only cares about the maximum derivative. Large 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 which bounds all derivatives.
(4.80) |
In matrix form,
Objective function becomes
(4.81) |
The result of adding the Chebyshev norm of diff(h) to the objective function to be minimized ( ) is shown in Fig.3.39. The result of increasing to 20 is shown in Fig.3.40.
L-One Norm of Derivative Objective
Another way to add smoothness constraint is to add -norm of the derivative to the objective:
(4.82) |
Note that the norm is sensitive to all the derivatives, not just the largest.
We can formulate an LP problem by adding a vector of optimization parameters which bound derivatives:
(4.83) |
In matrix form,
(4.84) |
The objective function becomes
(4.85) |
See Fig.3.41 and Fig.3.42 for example results.
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 method (§4.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