InterludeWorking with Mathematica: Optimal Filters
Because Signals and Systems is written in the form of Mathematica packages, it is automatically integrated with the rest of the capabilities of Mathematica. All of the computational and programming functions in Mathematica can be applied to problems in signal processing, while using Signals and Systems as a useful tool to speed development.
As an example of this procedure, we will investigate the design of optimal FIR filters by linear programming. Signals and Systems does not currently contain functions specifically for doing optimal FIR filter design, so functions must be designed in their entirety for this purpose. However, such functions can then be used in conjunction with the existing filter manipulation and analysis routines.
We begin by noting that the frequency response of a linearphase FIR filter of odd length is an even function. By symmetry considerations, we can write the response as
with
The function is called the amplitude. is defined as , where is the (odd) length of the filter. The coefficients are related in a fairly straightforward way to the impulse response.
In the design of an equiripple filter, we can apply the Chebyshev error criteria, defined as
where is the ideal (design) response. The error is only defined over the design frequencies in the passband and stopband, not over the transition band. This expression for the error can be written in terms of an inequality of the form
This, in turn, can be written as a pair of inequalities
Since is defined in the summation written earlier, the entire inequality can be nicely written in matrix form
where is a matrix of the cosines generated from the expression for by sampling at a number of frequencies, and computing the cosines at those same frequencies. The coefficients in the expression for appear as , the vector of . Subject to these constraints, we simply want to minimize .
This problem can be tackled by using the Mathematica function LinearProgramming.
Here is the usage message for LinearProgramming.
In[1]:=?LinearProgramming
We see that the equations we have are very close to this definition. The only remaining factor is that LinearProgramming requires that the variables being minimized be nonnegative, while we require that the coefficients be able to take on negative values. This situation is handled by writing each variable as the difference of two nonnegative variables (i.e., ). See any text on linear programming for a more detailed description of the technique.
To write the problem in a computable form, we need to know , which is an odd integer corresponding to the length of the filter, and we need samples of the ideal frequency response in the domain we care about. One convenient method for getting these samples is to sample the passband and the stopband separately at equispaced points. To do this, we need to know the edge frequencies, and the number of samples we want from the passband and from the stopband. It is convenient to work with normalized frequencies in radians per second, where the maximum frequency is Pi (we don't need to sample to 2 Pi because of the symmetry considerations used in defining the filter).
Begin by making certain Signals and Systems routines are available.
In[2]:=Needs["SignalProcessing`"]
Let passco be the cutoff frequency for the passband, and stopco be the cutoff frequency for the stopband.
In[3]:=passco = .6 Pi; stopco = .8 Pi;
Since we are designing a lowpass filter for this example, there is only one passband and one stopband.
Let npass be the number of samples in the passband, and nstop be the samples in the stopband.
In[5]:=npass = 8; nstop = 3;
Here are the samples of the ideal response in the area we care about.
In[7]:=ideal = Join[ Table[1, {npass}], Table[0, {nstop}] ]
Out[7]=
These are the frequencies at which the ideal samples were taken.
In[8]:=freqs = Join[ Range[0, passco, passco/(npass  1)], Range[stopco, Pi, (Pi  stopco)/(nstop  1)] ]//N
Out[8]=
Here is the filter length (which must be an odd integer), and the halflength, defined as in the equations above.
In[9]:=filtlen = 11; halflen = (filtlen  1)/2;
Given the value for and the sample frequencies, the matrix of cosines used in the expression for can be computed quite easily.
By the listability of Cos, we can readily create the matrix of products of frequencies and values of , and take the cosine of the entire matrix.
In[11]:=cosines = Cos[ Map[# Range[0, halflen] &, freqs ] ];
The number of rows of the matrix is the number of frequency samples, and the number of columns is .
In[12]:=Dimensions[cosines]
Out[12]=
We also need the column vector of coefficients corresponding to .
In[13]:=ones = Table[{1}, {Length[freqs]}]
Out[13]=
We now have all the elements we need to assemble the matrix and vectors for the linear program. A useful utility for creating a matrix from an array of matrices can be found in the package LinearAlgebra`MatrixManipulation`.
This loads some auxilliary matrix manipulation routines that come with Mathematica.
In[14]:=Needs["LinearAlgebra`MatrixManipulation`"]
Here is the routine we are interested in.
In[15]:=?BlockMatrix
The horizontal duplication of the arrays cosines and ones corresponds to the creation of dummy variables to meet the nonnegativity condition of LinearProgramming.
In[16]:=cc = BlockMatrix[ {{ cosines, ones, cosines, ones}, { cosines, ones, cosines, ones}} ];
The objective function is written as a weighting vector, equivalent to the coefficients of the variables. Again, the doubling of the vector allows us to meet the nonnegativity condition. Only the variables corresponding to have nonzero coefficients.
In[17]:=weights = Join[ Append[Table[0, {halflen + 1}], 1], Append[Table[0, {halflen + 1}], 1] ];
The linear program can now be set up and solved. Remember that this solution corresponds to the difference of pairs of nonnegative variables.
In[18]:=lp = LinearProgramming[ weights, cc, Join[ideal, ideal] ];
Here is the actual result.
In[19]:=lp = Take[lp, Length[lp]/2]  Drop[lp, Length[lp]/2]
Out[19]=
The vector that results from the linear program now contains the values for and . The impulse response of the digital filter can be assembled from the .
The symmetry considerations used in the definition of cause the values resulting from the linear program to be double that of the impulse response at all points but the (shifted) zero point. As a result, most values need to be divided
by .
In[20]:=impresp = Join[ Reverse[Take[lp, {2,2}]]/2, {First[lp]}, Take[lp,{2,2}]/2 ]
Out[20]=
We can now proceed with an analysis of the filter.
We first create an FIR filter object.
In[21]:=dfilt = DigitalFIRFilter[ impresp, n ];
Here is the frequency response of the filter. Note particularly the equiripple response and linear phase.
In[22]:=MagnitudePhasePlot[ DiscreteTimeFourierTransform[ dfilt, n, w ], {w, 0, Pi} ];
Here is a polezero plot of the filter object.
In[23]:=PoleZeroPlot[ ZTransform[dfilt, n, z] ]
Out[23]=
Success! We now have a method for designing equiripple linearphase FIR filters by linear programming. Once the necessary series of steps has been established, it is easy to place the code into a function that can be reused.
Needs["LinearAlgebra`MatrixManipulation`"]
FIRCoefficients[filterlength_Integer?OddQ,
passbandend_,
stopbandstart_,
{npass_Integer, nstop_Integer} ] :=
Module[{ideal, freqs, cc, weights, halflen = (filterlength  1)/2, lp},
ideal = Join[Table[1, {npass}], Table[0, {nstop}]];
freqs = Join[
Range[0, passbandend, passbandend/(npass  1)],
Range[stopbandstart, Pi, (Pi  stopbandstart)/(nstop  1)]
]//N;
cc = BlockMatrix[
{{ #1, #2, #1, #2},
{ #1, #2, #1, #2}}
]&[ Cos[ Map[# Range[0, halflen] &, freqs ] ],
Table[{1}, {Length[freqs]}]
];
weights = Join[
Append[#1, 1],
Append[#1, 1]
]&[ Table[0, {halflen + 1}] ];
lp = LinearProgramming[ weights, cc, Join[ideal, ideal]];
If[!VectorQ[lp, NumberQ], Return[$Failed]];
Join[
Reverse[Take[#, {2, 2}]]/2,
{First[#]},
Take[#,{2, 2}]/2
]&[ Take[lp, Length[lp]/2] 
Drop[lp, Length[lp]/2]
]
]
Here is the algorithm written as a single function, FIRCoefficients. Parts of the code have been modified to be more concise and to reduce the number of intermediate variables.
Note that care must be taken with selecting the number of frequency samples taken. Linear programming can be numerically unstable depending on the arrangements of the constraints; a check has been placed in the above program to verify that the call to LinearProgramming succeeded. Not all arrangements of frequency samples (which correspond to the constraints) will be successful, so you may need to experiment. A reasonable first choice is for the number of filter sample points to be the same as the length of the filter.
This form uses the same inputs that we specified at the start of this section. However, a more general design would require the band specification to be given as a FilterSpecification object, and place the number of samples for the frequencies in an option. It might also return a full filter object, instead of just the coefficients. These and other similar enhancements would make it more general and standardized for a commercial release. For personal use, however, the above routine is quite adequate.
Here is an example that uses the new FIRCoefficients function. A length lowpass filter is developed, with the transition band running from to . Constraint conditions are generated from samples taken in the passband and taken in the stopband.
In[26]:=MagnitudePhasePlot[ DiscreteTimeFourierTransform[ DigitalFIRFilter[ FIRCoefficients[ 19, 0.3 Pi, 0.5 Pi, {10, 9} ], n ], n, w ], {w, 0, Pi} ];
