Perhaps the finest paper written on filter design is Moorer's cleverly-titled opus, “The Manifold Joys of Conformal Mapping.” While primarily an invaluable study of the methods of transforming prototype filters into other filter types and to different frequencies, it also describes methods for the design of second-order shelf filters and higher-order shelf filters.
While Moorer's paper contains enough information to design these filters, the approach requires dealing with mathematics that might be unfamiliar to many.
Readily-available, modern tools such as Matlab can greatly simplify this task and make it at least as available to the average engineer as the nearest Matlab license. This paper presents a streamlined approach to designing these filters and contains a full Matlab implementation.
This paper closely follows Moorer's approach to the design of higher-order shelf filters. But since these filters are a special case of elliptic high and low pass filters, the difference being that the “stopband” of a shelf filter tends to some finite limit, this paper augments Moorer's development with constructs from Parks' and Burrus' method of elliptic filter design. Of further interest to those researching this area will certainly be a paper by Gray and Markel.
As with Moorer's approach, the approach outlined in this paper requires using numerical methods. However, rather than use of lesser-known, less readily available techniques, this paper relies upon Matlab's optimization functions for numerical solutions. The result is a straightforward Matlab function that is effective for designing shelf filters for a wide range of input parameters. This function and the supporting functions needed for a complete implementation are included in the appendix.
Following is an outline of Moorer's approach:
1. The user specifies a set of monotonically decreasing or increasing linear gains A, B, C, D which define the gains, respectively at 0, π/2, some intermediate frequency, θ, and π. Here π/2 and θ are the edges of the transition band for the first order filter to be designed in the next step, and θ > π/2. The values specified for B and C indicate the gain at these edge frequencies. However, in terms of the overall, higher-order filter being designed, they are equivalently thought of as the ripple specification, as shown in Figure 1.
Figure 1. Definitions of A, B, C, and D
2. The next step is the design of a first order digital filter that fits the specifications of Item 1. This filter has the following transfer function, which is Moorer's Eq. (17):
The parameters are derived by solving the equation at the four specified frequencies to get four equations and four unknowns, the three shown in Eq. (1) along with θ, which is allowed to float according to the parameters specified.
3. The first-order digital filter designed above is then mapped to the s-plane using the following equations, which are Eqs. (21b), (21c) from Moorer:
4. After mapping the prototype filter to the s-plane, it is mapped to the u-plane using Jacobi elliptic functions. The essence of this is that a parameter, k, is defined by Moorer's Eq. (24) to be
The u-plane roots un and ud are defined as follows in Moorers Eqs. (26a) and (26b):
where sc() denotes an elliptic function.
5. Now the desired order is finally achieved in the u-plane by defining a rectangle that appropriately encloses, in a symmetric pattern, N each of an infinite array of poles and zeros, where N is the desired order of the filter under design.
Moorer's paper provides alternative methods of defining this rectangle, but this paper focuses on the approach of specifying the filter order, and hence the size of the u-plane rectangle.
This requires finding the band edge, /θ , of the new Nth order filter. Moorer provides an approach, but this paper presents a simpler method in Section 2.
Moorer's Eq. (27) defines /k as a function of /θ
6. Moorer provides equations for mapping the N-dimensional u-plane filter to the s-plane, but Parks and Burrus's, Eqs. 7.84 and 7.85 are more straightforward:
7. Now these roots are mapped to the z-plane using the bilinear transform.
8. Finally, the gain is adjusted to the desired value and conformal mapping is used to map this filter to the desired cut frequency. Moorer's paper focuses on conformal mapping and explains the method thoroughly. His Eq. (2b) describes how to find mapping parameter a:
Here ρ corresponds to the frequency of the original filter that is to be mapped to frequency φ of the new filter. The smaller magnitude root, the one inside the unit circle, of Eq. 8 will be used. (As can be seen by rederiving it from Moorer's Eq. (2a), there is a typographical error in Moorer's Eq. (2b). This explains the slight difference between Moorer's Eq. (2b) and the corrected version, Eq. (8) of this paper.)
Once the value of a is found it is used in Moorer's Eq. (1) to make the transformation, which maps the prototype filter to the desired cut frequency:
The Matlab implementation mostly follows the above sketch while using the native functionality to do as much of the computation as possible. This section outlines the approach.
Step 1. Due to the way the overall scaling is achieved in the Matlab function, small changes are implemented in the definitions of A, B, C and D. These ensure that the tails of the designed filter exactly reach the specified gains. These modified equations are shown here in terms of the low-frequency gain specification, the low-frequency ripple specification, the high-frequency ripple specification and the high-frequency gain specification, all in dB:
Steps 2 and 3 are straightforward implementations of Moorer's equations, as shown in the Matlab code of the appendix. However, Step 4 is implemented using Matlab's native functionality, rather than the approach of Moorer.
Step 4. This step requires inverting Jacobi elliptic functions. While Matlab does have Jacobi elliptic functions built in, it does not provide their inverses. Moorer and other authors provide information on how to form this inverse, but, in every case, it is a numerical technique. Therefore, rather than constructing some such technique in Matlab, Matlab's own numerical abilities are used to perform a numerical inverse of these functions. This is accomplished by seeking the values of u that make the following equation equal to 0:
This is done using Matlab's fminsearch function, a numerical search algorithm. To solve Eq. (11), the following function is defined:
function e = ellmin(u, k, sigma)
[Sn, Cn, Dn] = ellipj(u, k^2);
e = abs(sigma – Sn/Cn);
This function is minimized using the following line from the main function:
un = fminsearch(@ellmin, 0.5, opt, kprime, sigman);
Performing this twice solves for un, and then ud (using inputs σn and σd, respectively).
Notice that function ellmin uses k2 rather than k as an input. This is simply a matter of definition. Matlab's Jacobi elliptic functions are defined in terms of m rather than k, where m = k2 .
Step 5 uses the newly-found values un and ud. The only tricky part of this step is the computation of /θ. There are different ways of looking at this, but this paper assumes pre-definition of the order, N. Now all that remains is to find the value of /θ that corresponds with this order. This is accomplished by optimizing equation thetaN, as described below, for /θ:
function e = thetaN(theta, N, K, Kprime)
khat = 1/tan(theta/2);
khatprime = sqrt(1 – khat^2);
Khat = ellipke(khat^2);
Khatprime = ellipke(khatprime^2);
e = abs(N – Khat*Kprime/Khatprime/K);
The optimization this time is accomplished with a different Matlab function, fminbnd, which looks for a solution between two bounds. Note the following line from the main function:
theta_N = fminbnd(@thetaN, pi/2+0.01, pi, opt, order, K, Kprime);
Notice that this syntax seeks the solution between something slightly higher than π/2, and π. The reason to limit the search to something greater than π/2 is that the fminbnd function will check the endpoints, and π/2 is not a valid solution.
Now that the value of /θ has been found, the rest of Step 5 can be carried out in a straightforward manner, using Matlab's ellipke function as needed to compute the complete elliptic integral.
Step 6. Substituting the Parks and Burrus equations for the Moorer equations, as described in the sketch above, and using Matlab's ellipj command to compute the needed Jacobi elliptic functions, this step is straightforward.
Step 7. A call to Matlab's bilinear function handles this step, as shown in the main function in the appendix.
Step 8. Moorer's method puts the beginning of the transition band at π/2, so Eqs. (8) and (9) can be used to map the beginning of the transition band to any desired frequency. However, this paper defines the corner frequency as that point 3 dB below (above) the maximum (minimum) amplitude shelf gain when the gain between shelves is 6 dB or higher. When the gain between shelves is less than 6 dB, the mid-gain point is used instead.
The simplest way to find the frequency at which this corner occurs in the higher-order prototype filter is to once again turn to Matlab command fminsearch to minimize function, magfind, which is minimum for the frequency corresponding to the commanded gain:
function e = magfind(omega, magdB, num, den)
e = abs(magdB – 20*log10(abs(N/D)));
Once found the corner frequency can readily be mapped to the desired frequency using Eqs. (8) and (9).
The final step in the design is to find the gain needed to scale the filter to the specified level.
From this point, depending upon the form needed for implementation, Matlab's zp2tf command can be applied to map this gain and the poles and zeros to an Nth-order transfer function.
The appendix includes a full Matlab function, ellipshelf.m, for designing higher-order shelf filters. This section shows results of using this function.
The following is the syntax for calling ellipshelf.m:
[B, A] = ellipshelf(order, lfgaindB, cutfreqHz, Fs, hfgaindB, lfrippledB, hfrippledB, figs)
B – returned numerator coefficient vector
A – returned denominator coefficient vector
order – desired filter order
lfgaindB – desired dB gain at low frequencies
cutfreqHz – desired 3 dB point
Fs – sample rate
hfgaindB – desired dB gain at high frequencies
lfrippledB – allowable low frequency shelf ripple in dB
hfrippledB – allowable high frequency shelf ripple in dB
figs – 1 for plots, 0 otherwise
Figs. 2 and 3 show respectively the magnitude and phase responses resulting from executing the function ellipshelf three times with the following syntax:
. ellipshelf(3, 6, 2000, 48000, 0, 0.001, 0.001, 1)
ellipshelf(6, 6, 2000, 48000, 0, 0.001, 0.001, 1)
ellipshelf(9, 6, 2000, 48000, 0, 0.001, 0.001, 1)
Figure 2. Magnitude Responses: Shelves of Increasing Orders
Figure 3. Phase Responses: Shelves of Increasing Orders
The solid line represents the 3rd order case, the dashed represents the 6th order case, and the dot-dashed line represents the 9th order case.
Figs. 4 and 5 show the magnitude and phase responses shown in Figs. 2 and 3 but zoomed in the frequency dimension to show detail.
Figure 4. Magnitude Responses: Shelves of Increasing Orders (zoomed)
Figure 5. Phase Responses: Shelves of Increasing Orders (zoomed)
Figs. 6 and 7 show respectively the (zoomed) magnitude and phase responses resulting from executing the function three times with the following syntax:
ellipshelf(7, 0, 10000, 96000, -5, 0.1, 0.1, 1);
ellipshelf(7, 0, 10000, 96000, -5, 0.01, 0.01, 1);
ellipshelf(7, 0, 10000, 96000, -5, 0.001, 0.001, 1);
Figure 6. Magnitude Responses: Shelves with Decreasing Ripple (zoomed)
The solid line is the first filter, with the ripple of 0.1 dB on both ends of the spectrum. This ripple is clearly visible in the magnitude response, as is the increased slope as compared to the other cases. This plot shows that it is possible to trade transition width for ripple, as would be expected.
Figure 7. Phase Responses: Shelves with Decreasing Ripples (zoomed)
The results above show that the presented function is highly effective for designing a wide variety of filters. However, as with many numerical techniques, there are those cases where this function will return something other than the desired filter. Fortunately it is easy to spot these problems from the plots. When such a problem occurs, reducing the order or changing the ripple specification will usually correct it.
The optimization routines used in the filter design have been adjusted for a good balance between design time and accuracy, but there might be value in adjusting them differently for certain applications.