## 1. Introduction

One of the main advantages of the linear-phase FIR filters over their IIR counterparts is the fact that there exist efficient algorithms for optimizing the arbitrary-magnitude FIR filters in the minimax sense. In case of IIR filters, the design of arbitrary-magnitude filters is usually time-consuming and the convergence to the optimum solution is not always guaranteed. The most efficient method for designing optimum magnitude linear-phase FIR filters with arbitrary-magnitude specifications is the Remez algorithm and the most frequent method to implement this algorithm is the one originally proposed by Parks and McClellan. Initially, they came up with the design of conventional low-pass linear phase FIR filters, whose impulse response is symmetric and the order is even [5]. Later on, along with Rabiner they extended the original design technique in [7] such that the generalized algorithm is applicable to the design of all the four linear-phase FIR filter types with arbitrary specifications. Due to the initial work of Parks and McClellan for the extended algorithm, this algorithm is famously known as the Parks-McClellan (PM) algorithm.

The PM algorithm was generated in the beginning of 1970 by using FORTRAN. During that era, the computer resources were quite limited. When people applied this algorithm in practice for high-order filters, they failed to achieve the optimum results. This gave rise to two main doubts about the PM algorithm. The first doubt was that the algorithm is not properly constructed and is valid for only low-order filters design. However, after noticing that the maximum number of iterations in the algorithm implementation is set to only 25 which is quite low to achieve the optimum solutions for high order filters, it became clear that the first doubt is superfluous.

The second doubt was concerned with the implementation of the PM algorithm’s search strategy for the “real” extremal points of the weighted error function, which is formed based on the “trial” extremal points. While mimicking the search technique included in the PM algorithm in various Remez-type algorithms for designing recursive digital filters [15, 16, 18, 19], it was observed that some of the them are quite sensitive to the selection of the initial set of the “trial” extremal points. Moreover, they suffered from certain issues which prevented them to converge at the optimum solution in some cases. For algorithms described in [15], and [18], the convergence issue was solved by increasing the maximum number of iterations in the above-mentioned manner. However, the algorithms described in [16] and [19] have still remained sensitive to the selection of the initial set of the “trial” extremal points. This sensitivity motivated the authors of this contribution to figure out how the search for the “true” extremal points is really carried out in the core discrete Remez algorithm part in the FORTRAN implementation of the PM algorithm. After a thorough study of the FORTRAN code, it was observed that this code utilizes almost twenty interlaced “go to” statements which are quite redundant for locating the “real” extremal points. Meticulous investigation of the code revealed that it is possible to decompose the one large chunk of the search technique into two compact search techniques referred to as *Vicinity Search* and *Endpoint Search* in such a manner that the same optimum solution can be achieved in an efficient manner as follows.

In *Vicinity Search*, the candidate “real” extremal point is located in the vicinity of each “trial” extremal point, which is bounded by the preceding and the following “trial” extremal points with the exception of the first (last) point for which the lower (upper) bound is the first (last) grid point in use. *Endpoint Search*, in turn, checks whether before the first (after the last) local extremum found by *Vicinity Search* there is an additional first (last) local extremum of opposite sign. If one or both of such extrema exist, then their locations are considered as candidate “real” extremal points of the overall search consisting of *Vicinity Search* and *Endpoint Search*. In this case, there are one or two more candidate “real” extremal points as *Vicinity Search* already provides the desired number of “real” extremal points. In the PM algorithm, the desired final “real” extremal points are determined according to the following three options:

*Option 1:* The additional last extremum exists such that its absolute value is larger than or equal to those of the first extremum of *Vicinity Search* and the possible additional first extremum.

*Option 2:* The additional first extremum exists such that its absolute value is larger than or equal to that of the first extremum of the *Vicinity Search* and larger than that of the possible additional last extremum.

*Option 3:* The conditions in *Options 1* and *2* are not valid.

For *Option 3*, the final “real” extremal points are the ones obtained directly from the *Vicinity Search*, whereas for *Option 1* (*Option 2*), these points are obtained by omitting the first (last) point based on *Option 3* and by replacing the last (first) point with the one found by *Endpoint Search*. Based on the above-mentioned facts, an extremely compact translation of the original FORTRAN implementation into a corresponding MATLAB implementation has been reporeted in [2].

The above study on how the PM algorithm performs the search for the “real” extremal points indicates that mimicking this search principle in the Remez-type algorithms proposed in [16], [19] does not give the flexibility to transfer two extremal points between the two consecutive bands, for example, from a passband to a stopband or vice versa, which is a necessary prerequisite for convergence to the optimum solution in certain cases. Most significantly, the search technique included in the PM algorithm does not follow the fundamental idea of the Remez multiple exchange (RME) algorithm when the approximation interval is a union of three or more disjoint intervals [8, 11, 20]. That is, if there are more candidate “real” extremal points than required, then the desired points should be selected in such a way that the ones corresponding to the largest absolute values of the weighted error functions are retained subject to the condition that the sign of the weighted error function alternates at the consecutive points. An efficient MATLAB based implementation following the above mentioned fundamental notion of the RME algorithm has been reported in [Ahsan Saramaki, 2011] and provides significant improvements in designing the multiband FIR filters.

In the beginning, the main purpose of the authors of this contribution was to modify the core discrete Remez part of the PM algorithm in FORTRAN such that it follows the fundamental principle of the RME algorithm and can ultimately be incorporated in the algorithms proposed in [16] and [19]. However, during the course of modifications, it was observed that a modified MATLAB implementation mimicking the original FORTRAN implementation is quite effective and superior to already available MATLAB implementation of the algorithm in the function firpm [21]. Based on the above discussion, this chapter describes two novel MATLAB based implementations of the Remez algorithm within the PM algorithm. Implementation I is an extremely fast and compact translation of the Remez algorithm part of the original FORTRAN code to the corresponding MATLAB code and is valid for general purpose linear-phase FIR filters design [2]. It is worth noting that Implementation I imitates the implementation idea of the Remez algorithm presented in PM algorithm. Implementation II is based on the fundamental notion of the Remez algorithm as described in [20] and provides significant improvements in designing the multiband FIR filters [3]. It is important to note that this chapter emphasizes on the practical MATLAB based implementation aspects of the Remez algorithm. In order to get an overview of the theoretical aspects of the Remez algorithm, please refer to [10, 17]. The organization of this chapter is as follows. Section (2) formally states the problem under consideration, Section (3) describes the Implementation I in detail, Section (4) discusses the Implementation II in detail, and finally, the concluding remarks are presented in Section (5).

## 2. Problem statement

After specifying the filter type, the filter order, and the filter specifications such that the problem is solvable using the RME algorithm, the essential problem in the PM algorithm is the following:

*Core Discrete Approximation Problem*: Given

, the number of unknowns

where

and

According to *Alternation (characterization) Theorem* [4, 12, 13],

(4) |

It is worth mentioning that the core discrete approximation problem is the same for both the implementations I and II as defined in the Introduction.

## 3. Implementation I

This section discusses Implementation I in detail as follows. First, the theoretical formulation of the algorithm is described so that the reader can grasp the very essence of the MATLAB code snippet provided later on in this section. Second, during this inclusion, it is emphasized that instead of using approximately 15 nested loops and around 300 lines of code, only 3 looping structures and approximately 100 lines of code are required by Implementation I. Third, it is shown, by means of four examples, that the overall CPU execution time required by the proposed implementation to arrive practically in the same manner at the optimum FIR filter designs is only around one third in comparison with the original implementation. Fourth, in the last two examples, there are unwanted peaks in the transition bands. In order to suppress these peaks to acceptable levels, two methods of including the transition bands in the original problem are introduced.

### 3.1. Theoretical formulation

The theoretical formulation of the proposed algorithm is roughly classified into the initialization phase and the iteration phase. The initialization phase performs the necessary initializations for the algorithm, whereas the iteration phase carries out the actual Remez exchange loop. In order to explain why Implementation I is a compact and efficient MATLAB based routine, the iteration phase is further decomposed into four well-defined primary segments. Each segment is constructed in such a way that before the start of the basic steps, there is a thorough explanation on the benefits of carrying out the segment under consideration with the aid of the proposed basic steps.

#### 3.1.1. Initialization phase

The overall implementation starts with the following initializations:

#### 3.1.2. Iteration phase

The Remez exchange loop iteratively locates the desired optimum vector *Alternation Theorem* is satisfied as follows. In the first loop, the vector

are satisfied. When concentrating only on the values of *Alternation Theorem*. However, when considering all the values of

The efficiency of Implementation I in comparison with the original MATLAB function firpm, in terms of significant reduction in the code compactness and a considerable reduction in the CPU execution time for obtaining practically in the same manner the best linear-phase FIR solutions are mostly based on the following two novel facts. First, the steps under Segment 1 are accomplished by employing the efficient MATLAB vectorization operations whenever possible and, most importantly, by avoiding the call for one subroutine by replacing this call with highly efficient matrix operations available in MATLAB. Second, as already mentioned in the introduction, the lengthy search technique involved in the function firpm for locating the true extremal points based on the weighted error function can be compressed into *Vicinity Search* and *Endpoint Search*. In the sequel, Segment 2 and Segment 3 will take care of *Vicinity Search* and *Endpoint Search*, respectively. More detail can be found in the actual implementations of Segments 1, 2, and 3.

Finally, Segment 4 checks whether

Segment 1: After knowing the “trial” vector *Alternation Theorem* is satisfied when concentrating only on those

is implicitly solved for the

In comparison with many scattered scalar operations in the original function firpm, the MATLAB code snippet, which is available in the following subsection, is computationally efficient and is highly compact due to the above-mentioned vectors. In addition to that, the time consuming subroutine of “remezdd” is replaced with simple and highly efficient matrix operations. Further improvements are obtained by using the vector grid, which contains the grid points under consideration, as well as des and wt, which carry information of the desired values and weights at these grid points, respectively. With the above-mentioned data, the weighted error function is generated only once during each iteration and is a single vector wei_err. This vector plays a pivotal role in the implementations of *Vicinity Search* in Segment 2 and *Endpoint Search* in Segment 3.

This segment is performed by using the following ten steps:

*Step 1:* Determine the entries of the vectors x and ad as

and

respectively, as well as the corresponding deviation value as

*Step 2:* Determine the entries of the vector y as

*Step 3:* Generate the entries of the abscissa vector x_all covering all the entries in the vector grid as

*Step 4:* Select the entries of the vectors err_num and err_den of length

*Step 5:* Generate the entries of the vector aid as

and update

*Step 6:* Set

*Step 7:* Generate the entries of the weighted error function wei_err for

The resulting wei_err contains undefined values at the entries of

*Step 8:* Set

*Step 9:* Update the vector wei_err as

*Step 10:* Set*Step 1* under Segment 2.

Segment 2: This segment explains how to perform *Vicinity Search* based on the values of the weighted error function *Vicinity Search* is to determine the

In the first phase, the search of both local minima and maxima is reduced to that of local maxima by multiplying the values of *Vicinity Search* by using the following three steps:

*Step 1:* Set

*Step 2:* Find the

where

and

achieves the maximum value. Generate*Endpoint Search* at Segment 3. Update

*Step 3:* Set*Step 1* under Segment 3.

Segment 3: This segment explains how to perform *Endpoint Search*. After *Vicinity Search*, the role of *Endpoint Search* is to check whether the weighted error function *Endpoint Search* is necessary to be used after *Vicinity Search* as *Vicinity Search* totally omits the existence of these possible additional extrema.

The appearance of the additional first local extremum implies that[5] -

is larger than or equal to

is smaller than or equal to

If no additional extremum exists, then the “final” *Vicinity Search*. Otherwise (that is, one or both the additional extrema exist), the final

*Alternative 1*: The additional last extremum exists such that its absolute value is larger than or equal to those of the first extremum found by *Vicinity Search* and the possible additional first extremum.

*Alternative 2*: The additional first extremum exists such that its absolute value is larger than or equal to that of the first extremum found by *Vicinity Search* and larger than that of the possible additional last extremum.

If *Alternative 1* (*Alternative 2*) is valid, then the final *Vicinity Search* is disregarded and the last (first) entry is the value of

The above explanation is the key idea to perform *Endpoint Search* in the function firpm. However, the function firpm performs *Endpoint Search* in a lengthier manner and in order to exactly follow this strategy, it is carried out by using the following eight steps:

*Step 1*: Set

*Step 2*: Determine

achieves the corresponding maximum entry value. Store this maximum entry value as

*Step 3*: If*Step 2* under Segment 2, then go to the next step. Otherwise, set

*Step 4*: Determine *Step 6*. Otherwise, find the index, denoted by

achieves its maximum entry value. Set

*Step 5*: If*Step 2* under Segment 2, then go to the next step. Otherwise, set

*Step 6*: If*Step 1* under Segment 4. Otherwise, go to the next step.

*Step 7*: If*Step 1* under Segment 4. Otherwise, go to the next step.

*Step 8*: If*Step 1* under Segment 4.

Segment 4: This concluding segment check the convergence of the Remez exchange loop as follows. If the entries of the vectors *Step 1* under Segment 1. Continue the Remez loop until

This segment requires only the following two basic steps:

*Step 1*: If

*Step 2*: Set *Step 1* under Segment 1.

### 3.2. MATLAB code snippet

The code pasted below has been tested for realizing Implementation I by using MATLAB version

This section presents a performance comparison between the Implementation I and the implementation available in the MATLAB function firpm. The performance measurement criteria is the time taken by both the implementations to design a particular filter and is measured with the help of MATLAB built-in function profiler. This function indicates the CPU execution time taken by a function and provides the following information:

Calls - The number of times the function was called while profiling was on.

Total Time - The total time spent in a function, including all child functions called, in seconds.

Self Time - The total time taken by an individual function, not including the time for any child functions called, in seconds.

Time measurement was carried out on an IBM ThinkCentre machine equipped with Intel Core 2 Duo processor E6550 running at a speed of 2.33 GHz with a memory of 3 GB.

The following four filters were designed with both implementations. After the examples, the time taken by them will be tabulated in Table 1.

Example 1: It is desired to design a lowpass filter meeting the following criteria:

The minimum order to meet these criteria is

The magnitude response of the resulting filter is shown in Fig. 1.

Example 2: It is desired to design a highpass filter meeting the following criteria:

The minimum order to meet these criteria is

The magnitude response of the resulting filter is shown in Fig. 2.

Example 3: It is desired to synthesize a bandpass filter meeting the following criteria:

The minimum order required to meet the criteria is

Figure 3 illustrates the magnitude response of the resulting filter. Although the amplitude response is optimal according to *Alternation Theorem*, it is worth noting that this particular filter has an extra peak in the second transition band region of approximately

Transition bandwidths are very large compared to the passband and/or stopband widths.

Width of transition bands is different; the larger the difference, the greater the problem.

In order to conveniently avoid the appearance of the unwanted transition peaks, consider an original problem stated for linear-phase Type I and Type II filters as follows.

First, there are

Such that these regions do not overlap. The lower and upper limits for the zero-phase frequency response in these bands are, respectively, specified as

Here, the

This region can be extended to cover the transition bands as follows:

where

In order to guarantee that *Type A* and *Type B* transition band constraints. For both types, the upper and lower limits for the zero-phase frequency response in the *Type A* and *Type B*, the upper limit is [7] -

whereas the lower limits depend on the type as follows:

The above limits for *Type A* are determined such that if the filter meets the overall criteria, then the maximum (minimum) value of the zero-phase frequency response in each transition band is less than or equal to the stated upper limit in the nearest passband region (larger than or equal to the stated lower limit in the nearest stopband region). For *Type B*, in turn, the upper limit is the same, whereas the lower limit is obtained from the upper limit by changing its sign, thereby indicating that the magnitude response of the filter is less than or equal to the stated upper limit in the nearest passband region.

The desired value in the

and

The following MATLAB function converts the original design specifications into those ones including either *Type A* or *Type B* transition band constraints as follows. The first three input parameters F_ori, Des_ori, and Dev_ori contain the *itype=1* (*itype=2*) means that *Type A* (*Type B*) transition band constraints are in use. The output of this function consists of vectors F, Des, and Wt that are in the desired form when calling the MATLAB function firpm in its original form or its modifications referred to as Implementation I or II in this contribution.

When using the above MATLAB function with*Type A* and *Type B* transition band constraints become

The corresponding desired and weight values for *Type A* and *Type B* transition band constraints are, respectively,

The relevant MATLAB commands

As seen in Fig. 4, by increasing the original filter’s order from 102 to 103, the transition bands constraints guarantee that the overall response of the filters stays within the desired limits.

Example 4: It is required to synthesize a bandstop filter meeting the following criteria:

The minimum order required to meet these criteria is

The magnitude response of the resulting bandstop filter is shown in Fig. 5. This response is the best one according to *Alternation Theorem*, but contains two extra peaks of approximately

As seen in Fig. 6, the overall response of the filters of the same order as the original one, that is, 102, stay within the desired limits for both *Type A* and *Type B* transition band constraints. Among these two constrained designs, *Type A* constrained design is preferred as for it the undershoot of the zero-phase frequency response is limited to be larger than or equal to

Table 1 indicates the outcomes obtained from the original implementation and the Implementation I of the Remez algorithm, both of which work practically in the same manner. It is evident that the time required by the Implementation I is almost one third of the time taken by the original implementation and illustrates the superiority of the proposed MATLAB implementation of the algorithm.

## 4. Implementation II

This section discusses the Implementation II in detail. First, a theoretical formulation of the algorithm is presented and, then, the corresponding MATLAB code is specified. After that, a detailed comparison shows how this implementation is superior to the original implementation of the Remez algorithm in the function firpm, in terms of significant reductions in the number of iterations and the CPU execution time required to generate the same optimum solution, especially in multiband cases.

### 4.1. Theoretical formulation

As mentioned in the introduction, the key difference between Implementations I and II is the search strategies employed for the “real” extremal points of the weighted error function, which is formed based on the “trial” extremal points. Consequently, Segment 1 and Segment 4 are same for both the implementations. The remaining Segment 2 and Segment 3 along with the accompanying steps are as follows.

Segment 2: Based on the values of

*Condition 1*: At each entry of

*Condition 2*: In case of several consecutive local extrema of

*Condition 3*: The sign of

This vector

*Step 1:* Find all the values of

*Step 2:* Extract those entries from

*Step 3:* Split the range

*Step 4:* Generate the vector *Step 1* under Segment 3.

Segment 3: Based on the vector

*STEP A*: Denote the length of

*STEP B*: Denote the remaining vector and its length by

*STEP C*: Since *STEP C*.

This segment is carried out using the following seven steps:

*Step 1:* Set*Step 1* under Segment 4. Otherwise, go to the next step.

*Step 2:* If *Step 4*. Otherwise, go to the next step.

*Step 3:* If

*Step 4:* If*Step 1* under Segment 4. Otherwise, go to the next step.

*Step 5:* Determine the entries of the vector wei_comp as follows:

(40) |

*Step 6:* If *Step 4*. Otherwise, go to the next step.

*Step 7:* Find the entry of wei_comp with the smallest value. If this is the *Step 4*.

### 4.2. MATLAB code snippet

### 4.3. Performance comparison

The code pasted below has been tested and implemented with MATLAB version 7.11.0.584 (R2010b). It can be embedded into the MATLAB function firpm in a similar fashion to Implementation I.

This subsection shows how the proposed Implementation II, following the fundamental principle of the RME algorithms, outperforms the original implementation in terms of significant reductions in both the number of iterations and the CPU execution time required to arrive at the same optimum solution. For this purpose, four practical filter design examples are discussed. In all these examples, the problem is to design a filter having five interlaced passbands and stopbands. In order to achieve the accepted behavior in the transition band regions, the last two examples require the use of the *Type A* or *Type B* transition band constraints described in Subsection 3.3.

Example 5: It is desired to design a five-band filter with two passbands and three stopbands meeting the following specifications:

(41) |

The minimum order to meet the criteria is

The magnitude response of the resulting filter is shown in Fig. 7.

Example 6: It is desired to design a five-band filter with three passbands and two stopbands with following specifications:

(42) |

The minimum order to meet the criteria is

The magnitude response of the resulting filter is shown in Fig. 8.

Example 7: It is desired to design a five-band filter with two passbands and three stopbands with following specifications:

(43) |

The minimum filter order required to meet these specifications is*Type A* and *Type B* transition band constraints in the approximation problem according to the discussion of Subsection 3.3. When using*Type A* and *Type B* constraints is

The relevant MATLAB commands for designing the above-mentioned three filters are

Example 8: It is desired to design a five bands filter with three passbands and two stopbands with following specifications:

(44) |

The minimum filter order to meet these specifications is*Type A* and *Type B* constraints, the minimum filter orders are

The relevant MATLAB commands for designing the above-mentioned three filters are

The number of iterations as well as the CPU execution times required by the original implementation and the proposed second implementation are summarized in Table 2 for synthesizing all the eight filters considered in this section. The hardware and MATLAB versions are the same as used earlier during the comparison of the original and the proposed first implementations. It is seen that the reduction in the numbers of iterations is 79 and 54 percent, respectively, when synthesizing the filters in Example 5 and 6. In case of examples 7 and 8 with the transition band constraints in effect, the reduction in the numbers of iterations is 45 percent (*Type A*), 38 percent (*Type B*) and 18 percent (*Type A*), and 55 percent (*Type B*), respectively. The reduction percentage in the CPU execution time is between 71 and 86 percent for all the eight filter designs under consideration. Hence, the proposed second implementation is highly efficient in the design of multiband filters.

## 5. Conclusion

This chapter has introduced two novel MATLAB based Remez algorithms for the design of optimum arbitrary-magnitude linear-phase FIR filters. The first algorithm is a highly optimized and compact translation of the PM algorithm from its original FORTRAN code to its MATLAB counterpart in comparison with the existing MATLAB function firpm. These attractive properties have been achieved by first observing that the PM algorithm’s very lengthy search strategy for the “real” extremal points of the weighted error function, which is formed based on the “trial” extremal points, can be compressed into two very compact basic search techniques. Second, the MATLAB vectorization techniques are employed whenever possible. As a result, the CPU execution time is roughly one third to synthesize linear-phase FIR filters practically in the same manner in comparison with the function firpm being more or less a direct translation from the FORTRAN code. Moreover, the code complexity is reduced to a considerable extent. The original implementation utilizes approximately 15 nested loops and around 300 lines of code whereas the first proposed implementation uses only 3 looping structures and approximately 100 lines of code. Thus, same efficient results are achieved with one fifth of the looping structures and one third of the code lines in the original implementation.

It is, however, important to note that the first technique does not follow the fundamental idea of Remez algorithm as suggested in [20] as it tries to find the new “trial” extremal points in the vicinity of previously found points as well as in the surroundings of the first and last grid points under consideration.

The second implementation obeys the fundamental principle of the Remez multiple exchange algorithm. This means that while searching for the “real” set of extrema, there is a possibility to obtain more than the required points in intermediate calculations. In this situation, the idea is to keep as many extremal points as possible subject to the condition that the corresponding error values are the maximum absolute ones and they obey the sign alternation property. Another prominent feature is that the weighted error function is calculated over the entire grid. This, not only makes sure that no potential extremal frequency point is skipped during a particular iteration, but also enables to transfer the two extremal points between the consecutive bands which, in some cases, is a necessary prerequisite for the algorithms in [16] and [19] to converge. Furthermore, the number of iterations as well as the CPU execution times required by the proposed second implementation to design the linear-phase FIR filters in comparison with the existing MATLAB function firpm, especially in multi-band cases, are significantly lower. Examples have shown that in most five-band cases with some constraints in the transition bands, the reduction in the number of iteration is more than 50 percent, whereas the reduction in the CPU execution time is around 80 percent.

The quality of the filters designed with the proposed implementations is analogous to that of the PM algorithm with the added advantages of less number of iterations and CPU execution time.

The proposed two implementations have concentrated only on the core discrete Remez part of the PM algorithm. Future work is devoted to explore the possibilities of further improvements in the overall function firpm and reimplementing the other portions of this function efficiently.