Open access

# Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design

Written By

Muhammad Ahsan and Tapio Saramäki

Submitted: December 8th, 2011 Published: September 26th, 2012

DOI: 10.5772/46477

From the Edited Volume

## MATLAB

Edited by Vasilios N. Katsikis

Chapter metrics overview

View Full Metrics

## 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 nz1

In this contribution, nz1is chosen to be the number of adjustable parameters in both the FORTRAN and the MATLAB implementations of the PM algorithm because in this case nz stands for the number of extrema at which Alternation Theorem should be satisfied in order to guarantee the optimality of the solution.

, the number of unknowns a[n] forn=0,1,,nz2, and the grid points grid(k) included in the vector grid of lengthngrid, which contains values between 0 and 0.5

In the original PM algorithm, this range is the baseband for the so-called normalized frequencies from which the corresponding angular frequencies are obtained by multiplying these frequencies by 2π [17].,

along with the vectors des and wt of the same lengthngrid, the entries of which carry the information of the desired and weight values, respectively, at the corresponding grid points of the vector grid, find the unknowns a[n] to minimize the following quantity:

ε=max1kngrid|E(k)|,E1

where

E(k)=wt(k)[G(k)des(k)]E2

and

G(k)=n=0nz2a[n]cos[2πngrid(k)].E3

According to Alternation (characterization) Theorem [4, 12, 13], G(k)of the form of (3) is the best unique solution minimizing ε as given by (1) if and only if there exists a vector opt that contains (at least) nzentries opt(1),opt(2),,opt(nz) having the values of k within 1kngrid such that

opt(1)<opt(2)<<opt(nz1)<opt(nz)E[opt(m+1)]=E[opt(m)]form=1,2,,nz1|E[opt(m)]|=εform=1,2,,nz.E4

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:

• Initialize the element values of ”trial” vector trial of length nz for opt as trial(m)=1+(m1)(ngrid1)/nz for m=1,2,,nz1 andtrial(nz)=ngrid. Here, xstands for the integer part ofx.

• Initialize the iteration counter as niter=1 and set the maximum number of iterations asitrmax=250.

#### 3.1.2. Iteration phase

The Remez exchange loop iteratively locates the desired optimum vector trial having as its entries the nz values of k within1kngrid, at which Alternation Theorem is satisfied as follows. In the first loop, the vector trial found in the initialization phase is the first “trial” vector for being the desired optimum vectoropt. As this is extremely unlikely to happen, Segment 1, based on the vectoropt, generates the weighted error vector wei_err(k) for 1kngrid corresponding toE(k), as given by (2), in such a way that the nz1 unknowns a[0],a[1],,a[nz2] as well as dev

It should be noted that the value of dev is either positive or negative

are implicitly found so that the following nz equations

wei_err(trial(m))=(1)m+1devform=1,2,,nzE5

are satisfied. When concentrating only on the values of k being the enteries of the “trial” vectortrial, this solution is the best one according to Alternation Theorem. However, when considering all the values of k within1kngrid, this solution is not the best one.

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 realtrial or not. If this equivalence is established, then the best solution has been found as in this caseoptrealtrial. Otherwise, the whole process is repeated by using the “real” vector real of the present iteration as a “trial” vector for the next iteration. This exchange of the vectors is continued until real and trial coincide or the maximum allowable number of the iterations is exceeded, which is extremely unlikely to occur.

Segment 1: After knowing the “trial” vector trial that contains the nz trial values of k in the ascending order for 1kngrid in the present iteration, this first segment guarantees that Alternation Theorem is satisfied when concentrating only on those nz values of k being involved in the vectortrial. For this purpose, it generates the weighted error vector wei_err(k) for 1kngrid such that the following system of nz equations:

wt(trial(m))[m=0nz2a(n)cos(2πngrid(trial(m)))des(trial(m))]=(1)m+1devform=1,2,,nzE6

is implicitly solved for the nz1 unknowns a[0],a[1],,a[nz2] as well as fordev.

is worth emphasizing that implicit solutions for the calculation ofa[n]’s are not required for the intermediate iterations. The explicit solution for the calculation ofa[n]’s is needed only after achieving the convergence to the best solution.

For this purpose, similar to the function firpm, for a given “trial” vectortrial, the value of wei_err(k) atk=trial(1), denoted bydev, the corresponding abscissa vector x, the ordinate vector y, and the coefficient vector ad, each of which are of lengthnz, are determined. These vectors are required to express the zero-phase frequency response when using the Lagrange interpolation formula in the barycentric form at each value of k for1kngrid, thereby making the implementation of the Remez loop very accurate and efficient.

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

x(m)=cos[2πtrial(m)]form=1,2,,nzE7

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

x_all(k)=cos[2πgrid(k)]fork=1,2,,ngrid.E11

Step 4: Select the entries of the vectors err_num and err_den of length ngrid to be zero valued, setm=1, and go to the next step.

Step 5: Generate the entries of the vector aid as

and update err_num(k) and err_den(k) for k=1,2,,nz as err_num(k)=err_num(k)+y(m)aid(k) anderr_den(k)=err_den(k)+aid(k), respectively.

Step 6: Setm=m+1. Ifm>nz, then go to the next step. Otherwise, go to the previous step.

Step 7: Generate the entries of the weighted error function wei_err for k=1,2,,ngrid as

wei_err(k)=[err_num(k)/err_den(k)des(k)]wt(k).E13

The resulting wei_err contains undefined values at the entries of trial due to the use of the Lagrange interpolation formula in the barycentric form. The undefined values can be conveniently filled based on the fact that at trial(m) with m odd (even), the desired value isdev(dev), where dev is given by (8). Hence, the vector wei_err can be completed by using the following three steps:

Step 8: Set m=1 and go to the next step.

Step 9: Update the vector wei_err as

wei_err(trial(m))=(+devformodddevformeven.E14

Step 10: Setm=m+1. Ifm<nz+1, then go to the previous step. Otherwise, go to Step 1 under Segment 2.

Segment 2: This segment explains how to perform Vicinity Search based on the values of the weighted error function wei_err(k) for1kngrid, which has been generated at Segment 1, and the “trial” vectortrial, which is under consideration in the present iteration. The key idea in Vicinity Search is to determine the mth entry of the “real” vectorreal, denoted by real(m) form=1,2,,nz, to be the value of k in the close vicinity ofk=real(m), where a local extremum of wei_err(k) with the same sign occurs. The location of these nz entries are simplified as follows.

In the first phase, the search of both local minima and maxima is reduced to that of local maxima by multiplying the values of wei_err(k) for 1kngrid by sign[wei_err(trial(m))] as in this case the values of the resulting signed weighted function sign[wei_err(trial(m))]×wei_err(k) become positive at k=trial(m) as well as in its proximity. In the second phase, the proper location of each real(m) for m=1,2,,nz can be obtained conveniently based on the following facts. First Segment 1 guarantees that form>1[m<nz], the “signs” of wei_err(k) at both k==trial(m1) and k=trial(m+1) is opposite to that ofk=trial(m), or correspondingly, atk=real(m). Second, during the course of the present search, real(m1)is located before real(m) in such a way that the sign of wei_err(k) at k=real(m1) is opposite to that of the sign of wei_err(k) atk=real(m+1). The above mentioned facts together with the reasoning that the lowest value of k for locating real(1) is 1 and the highest value of k for locating real(nz) is ngrid inherently lead to carrying out Vicinity Search by using the following three steps:

Step 1: Set m=1 and go to the next step.

Step 2: Find the mth element, denoted by˜real(m), at which the vector

err_vicinity=sign[wei_err(trial(m))]wei_err(low:upp),E15

where

low=(1form=1max{trial(m1)+1,real(m1)+1}for2mnzE16

and

upp=(trial(m+1)1for1mnz1ngridform=nzE17

achieves the maximum value. Generatereal(m)=˜real(m)+low1. Ifm=1[m=nz], then store the corresponding maximum value as err_vic(1)[err_vic(nz)] to be used in Endpoint Search at Segment 3. Update real(m) asreal(m)=˜real(m)+low1.

Step 3: Setm=m+1. Ifm<nz+1, then go to the previous step. Otherwise, go to 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 wei_err(k) contains an additional local extremum before k=real(1) [afterk=real(nz)] such that its sign is opposite to that of occurring atk=real(1)[k=real(nz)]. It is worth emphasizing that in order to take into account all the candidate extrema, 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

Vicinity Search automatically guarantees that the sign of the weighted error function wei_err(k) is same at both k=trial(1) andk=real(1). Hence, uppendis the smallest value ofk, where the sign of wei_err(k) can be opposite before these values. Similarly, the sign of wei_err(k) is same at both k=trial(nz) and k=real(nz) and lowend as given by (18) is the largest value ofk, where the sign of wei_err(k) can be opposite after these values.

uppend=min{trial(1)1,real(1)1}E18

is larger than or equal to1. If this fact holds true, then the largest entry in the sub-vector sign[wei_err(real(1))]wei_err(1:uppend) should be positive. Similarly, the existence of the additional last local extremum implies that

lowend=max{trial(nz)+1,real(nz)+1}E19

is smaller than or equal tongrid. If this fact holds true, then the largest entry in the subvector sign[wei_err(real(nz))]wei_err(lowend:ngrid) should be positive.

If no additional extremum exists, then the “final” realis the one found in Vicinity Search. Otherwise (that is, one or both the additional extrema exist), the final real is constructed according to the following alternatives:

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 real is formed such that the first (last) entry of real of Vicinity Search is disregarded and the last (first) entry is the value of k for1kngrid, where the additional last (first) maximum of the signed weighted error function sign[wei_err(real(nz))]wei_err(k) [sign[wei_err(real(1))]wei_err(k)] occurs.

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: Setendsearch=0.

Step 2: Determineuppendaccordingto(17). Ifuppend=0, then seterr_end(1)=0. Otherwise, find the index, denoted byend_real(1), where the vector

err_endpoint=sign[wei_err(real(1))]wei_err(1:end(1))E20

achieves the corresponding maximum entry value. Store this maximum entry value aserr_vic(1).

Step 3: Iferr_end(1)<err_vic(nz), where err_vic(nz) has been saved at Step 2 under Segment 2, then go to the next step. Otherwise, setendsearch=1.

Step 4: Determine lowend according to (14). Iflowend=ngrid+1, then go to Step 6. Otherwise, find the index, denoted by˜end_real(nz), where the vector

err_endpoint=sign[wei_err(real(nz))]wei_err(lowend:ngrid)E21

achieves its maximum entry value. Set end_real(nz)=˜end_real(nz)+lowend1 and store the corresponding maximum entry value aserr_end(nz).

Step 5: Iferr_end(nz)<max{err_end(1),err_vic(1)}, where err_vic(1) has been saved at Step 2 under Segment 2, then go to the next step. Otherwise, setendsearch=2.

Step 6: Ifendsearch=0, then go to Step 1 under Segment 4. Otherwise, go to the next step.

Step 7: Ifendsearch=1, then set real(nz+1m)=real(nzm) for m=1,2,,nz1 and real(1)=end_real(1) and got to Step 1 under Segment 4. Otherwise, go to the next step.

Step 8: Ifendsearch=2, set real(m)=real(m+1) for m=1,2,,nz1 andreal(nz)=end_real(nz). Go to 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 trial and real are the same, then stop as in this case the ultimate goal optrealtrial has been achieved. Otherwise, use the present “real” vector real as the “trial” vector for the subsequent iteration by using the substitution trial=real and go to Step 1 under Segment 1. Continue the Remez loop until trial and real coincide or the value of the iteration counter niter exceedsitrmax=250, which is extremely unlikely to occur.

This segment requires only the following two basic steps:

Step 1: If trial and real coincide, then stop. Otherwise, go to the next step.

Step 2: Set trial=real andniter=niter+1. Ifniter>itrmax, then stop. Otherwise, go to Step 1 under Segment 1.

### 3.2. MATLAB code snippet

The code pasted below has been tested for realizing Implementation I by using MATLAB version7.11.0.584(R2010b). In order to embed this code snippet in the MATLAB function firpm, edit the function by taking away the code between lines 214 and514, introduce the function call (the first line of the function of remez_imp1) and copy the function at the end of the file. Remember to take the backup copy of the original function. It is worth emphasizing that the implementation of Remez algorithm in the function firpm uses approximately 15 nested loops and 300 lines of code, whereas the code snippet provided below requires only 3 looping structures and approximately 100 lines of code to achieve the same optimum solution.

3.3. Performance comparison

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:

ωp=0.05π,ωs=0.1π,δp=0.01,andδs=0.001.E22

The minimum order to meet these criteria is 108 and the relevant MATLAB commands are

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:

ωs=0.02π,ωp=0.05π,δp=0.01,andδs=0.001.E23

The minimum order to meet these criteria is 172 and the relevant MATLAB commands are

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:

ωs1=0.2π,ωp1=0.25π,ωp2=0.6π,ωs2=0.7π,δp=δs2=0.01,andδs1=0.001.E24

The minimum order required to meet the criteria is 102 and the relevant MATLAB commands are

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 16 dB. This is because of the fact that the approximation interval is a union of passband and stopband regions and transition bands are considered as “don’t care” regions. This assumption works perfectly for filters having bands less than three. However, in case of three or more bands, there is no guarantee that the response is well-behaved in the transition bands, even though it is optimal according to the approximation theory. This fact is especially prominent if any one of the following holds true [9]:

• 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 R interlaced passband and stopband regions as given by

Ωκ=[ωκ(low),ωκ(upp)]forκ=1,2,,R.E25

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

Lκ(low)=DκδκandLκ(upp)=Dκ+δκE26

Here, theDκ’s alternatingly achieve the values of zero and unity such that the first value is zero (unity) if the first band is a stopband (passband). For this original problem, the overall approximation region is

Ω=Ω1Ω2ΩR.E27

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

Ω^=Ω1Ω1TΩ2Ω2TΩR1TΩR,E28

where

ΩκT=[ωκ(upp)+α,ωκ(low)α]forκ=1,2,,R1.E29

In order to guarantee that Ω^κ is still a closed subset of[0,π], αshould be a small positive number.

The only condition for α is that it should be small enough to avoid the extra peaks between the adjacent passbands and the newly formed intervals in the transition band regions.

There are two natural ways to state the transition band constraints, referred to as Type A and Type B transition band constraints. For both types, the upper and lower limits for the zero-phase frequency response in the κth transition band ΩκT are specified as follows. For both Type A and Type B, the upper limit is

It is worth emphasizing that the use of max{Dκ+δκ,Dκ+1+δκ+1} implies that the maximum allowable value in the nearest passband is the upper limit. Similarly, in the following equation, min{Dκδκ,Dκ+1δκ+1}means that the lower limit is the one in the nearest stopband.

L^κ(upp)=max{Dκ+δκ,Dκ+1+δκ+1},E30

whereas the lower limits depend on the type as follows:

L^κ(low)=(min{Dκδκ,Dκ+1δκ+1},forTypeAL^κ(upp)forTypeBE31

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 κth transition band ΩκT for both types is the average of the corresponding lower and upper limits, whereas the admissible deviation is the difference between the upper limit and the above-mentioned desired value. Hence, in ΩκT for κ=1,2,,R1 the desired values, denoted byD^κ, and the admissible deviations, denoted byδ^κ, are as follows:

D^κ=([max{Dκ+δκ,Dκ+1+δκ+1}+min{Dκδκ,Dκ+1δκ+1}]/2forTypeA0forTypeB.E32

and

δ^κ=(max{Dκ+δκ,Dκ+1+δκ+1}D^κforTypeAmax{Dκ+δκ,Dκ+1+δκ+1}forTypeB.E33

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 2R edges of the R bands as a fraction of π as well as the desired values and the admissible deviations from these values in the R bands in the original specifications.. "alpha" corresponds directly to α which is used in (25), whereas 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α=0.0005, the ten band-edges of the five bands as fractions of π for both Type A and Type B transition band constraints become

Ω^κ=[0,0.2,0.2005,0.2495,0.25,0.6,0.6005,0.6995,0.7,1].E34

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

DA=[0,0.5045,1,0.5,0]E35

WA=[1000,1.9782,100,1.9608,100],E36
and

DB=[0,0,1,0,0]E37
WB=[1000,0.9901,100,0.9901,100].E38

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:

ωp1=0.15π,ωs1=0.3π,ωs2=0.6π,ωp2=0.65π,δp1=δp2=0.01,andδs=0.001.E39

The minimum order required to meet these criteria is 102 and the relevant MATLAB commands are

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 33 and 15 dB in the first transition band. By using the technique described above, the transition band peaks can be attenuated to an acceptable level. The relevant MATLAB commands are

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δs=0.001. Furthermore, the response in the first passband remains equiripple.

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.

 Example Original Implementation Implementation I Total Time* Self Time Total Time Self Time 1 0.152 0.125 0.032 0.032 2 0.235 0.136 0.064 0.064 3Τ 0.159 0.120 0.032 0.032 3⊥A 0.262 0.191 0.056 0.056 3⊥B 0.307 0.184 0.061 0.061 4Τ 0.198 0.138 0.047 0.047 4⊥A 0.272 0.169 0.051 0.051 4⊥B 0.318 0.186 0.055 0.055

### Table 1.

Performance Comparison of Original Implementation and Implementation I

## 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 wei_err(k) for 1kngrid generated at Segment 1, this main step generates the vector real(start) to include as many entries as possible in the ascending order subject to the following three conditions:

Condition 1: At each entry ofreal(start), wei_err(k), when regarded as a function ofk, achieves a local extremum whose absolute value is larger than or equal to|dev|, where |dev| is determined according to (8).

Condition 2: In case of several consecutive local extrema of wei_err(k) with the same sign fork(low)kk(upp), only one entry is included in real(start) and its value is the value of k fork(low)kk(upp), where |wei_err(k)| achieves its maximum value.

Condition 3: The sign of wei_err as a function of k alternates at the consecutive enteries ofreal(start).

This vector real(start) serves as a start-up vector for generating the “real” vector real at Segment 3. This segment is carried out using the following four steps:

Step 1: Find all the values of k in the range1kngrid, where a local extremum of wei_err(k) occurs. Store these values of k in the ascending order into the vectoraid1.

Step 2: Extract those entries fromaid1, where the absolute value of wei_err is larger than or equal to|dev|, and store these entries into the vectoraid2.

Step 3: Split the range1κnaid2, where naid2 is the length ofaid2, into the sub-ranges κ(low)(m)κκ(upp)(m) for m=1,2,,naid2 in such a way that the signs of wei_err(aid2(κ)) in the mth sub-range as given by κ(low)(m)κκ(upp)(m) are the same.

Step 4: Generate the vector real(start) of length nzstart such that its mth entry is the value among aid2(κ) forκ(low)(m)κκ(upp)(m), at which |wei_err(aid2(κ))| achieves its maximum value. Go to Step 1 under Segment 3.

Segment 3: Based on the vector real(start) generated at Segment 2, this main step extracts from its enteries the nz enteries to be included in the “real” vector real such that the largest absolute values ofwei_err(k), when regarded as a function ofk, are retained subject to the condition that the maxima and minima alternate at the consecutive extracted enteries. If the length of real(start) isnz, then realreal(start) and no extraction is required. Otherwise, the extraction is performed using the following steps:

STEP A: Denote the length of real(start) bynreal(start). If nreal(start)nz is an odd integer, then the first (last) entry is discarded from real(start) if the absolute value of wei_err(k) at k=real(start)(1) is less than or equal (larger) than the corresponding value atk=real(start)(k=nreal(start)). Go to the next step.

STEP B: Denote the remaining vector and its length by ˜real(start) andn˜real(start), respectively. Ifn˜real(start)nz=0, then stop. Otherwise go to the next step.

STEP C: Since n˜real(start)nz is an even integer, two entries of ˜real(start) should be simultaneously discarded. There are altogether n˜real(start) optional pairs such that the first pair consists of the first and last entries of ˜real(start) and the remaining ones are n˜real(start)1 consecutive entry pairs. The pair to be discarded is the pair, where the maximum of two absolute values of wei_err(k) is the smallest. Go to STEP C.

This segment is carried out using the following seven steps:

Step 1: Setreal(init)=real(start). Ifnzreal(init)nz, where nzreal(init) is the length ofreal(init), is zero, then set real=real(init) and go to Step 1 under Segment 4. Otherwise, go to the next step.

Step 2: If nzreal(init)nz is an even integer, then go to Step 4. Otherwise, go to the next step.

Step 3: If|wei_err(real(init)(1))||wei_err(real(init)(nzreal(init)))|, then discard the first entry fromreal(init). Otherwise, discard the last entry fromreal(init). Go to the next step.

Step 4: Ifnzreal(init)nz, where nzreal(init) is the length of the remaining vectorreal(init), is zero, then set real=real(init) and go to Step 1 under Segment 4. Otherwise, go to the next step.

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

wei_comp(m)=max(|wei_err(real(init)(m))|,|wei_err(real(init)(m+1))|)form=1,2,,nzreal(init)1.E40

Step 6: If max(|wei_err(real(init)(1))|,|wei_err(real(init)(nzreal(init)))|) is less than or equal to the largest entry of wei_comp, then remove the first and last entries from real(init) and go to Step 4. Otherwise, go to the next step.

Step 7: Find the entry of wei_comp with the smallest value. If this is the mth entry, then remove the mth and the (m+1)th entries fromreal(init). Go to Step 4.

### 4.2. MATLAB code snippet

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.

### 4.3. Performance comparison

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:

ωs1=0.17π,ωp1=0.23π,ωp2=0.47π,ωs2=0.53π,ωs3=0.67π,ωp3=0.73π,ωp4=0.82π,ωs4=0.88π,δs1=δs2=δs3=0.001,andδp1=δp2=0.01.E41

The minimum order to meet the criteria is 91 and the relevant MATLAB commands are

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:

ωp1=0.1π,ωs1=0.15π,ωs2=0.3π,ωp2=0.35π,ωp3=0.75π,ωs3=0.8π,ωs4=0.85π,ωp4=0.9π,δp1=δp2=δp3=0.01,andδs1=δs2=0.001.E42

The minimum order to meet the criteria is 106 and the relevant MATLAB commands are

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:

ωs1=0.15π,ωp1=0.2π,ωp2=0.45π,ωs2=0.55π,ωs3=0.7π,ωp3=0.8π,ωp4=0.85π,ωs4=0.93π,δs1=δs2=δs3=0.001,andδp1=δp2=0.01.E43

The minimum filter order required to meet these specifications is100. The magnitude response of the resulting filter designed without any constraints in the transition bands is shown by the solid blue line in Fig. 9. It is observed that in the second and third transition bands there are unwanted peaks of approximately 9 dB and 16 dB, respectively. These undesired peaks can be avoided by using Type A and Type B transition band constraints in the approximation problem according to the discussion of Subsection 3.3. When usingα=0.0005, the minimum order to meet the resulting criteria for both Type A and Type B constraints is101. The responses of the resulting filters meeting these additional constraints are depicted in Fig. 9 by using a dashed red line and a dot-dashed black line, respectively.

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:

ωp1=0.17π,ωs1=0.27π,ωs2=0.47π,ωp2=0.52π,ωp3=0.69π,ωs3=0.79π,ωs4=0.87π,ωp4=0.92π,δp1=δp2=δp3=0.01,andδs1=δs2=0.001.E44

The minimum filter order to meet these specifications is102. The magnitude response of the resulting filter designed without any constraints in the transition bands is shown by the solid blue line in Fig. 10. It is observed that in the first, second, and third transition bands, there are undesired peaks of approximately 1 dB, 2dB, and 1.7 dB, respectively. These undesired peaks can be avoided in a manner similar to that used in the previous example. In this example, for Type A and Type B constraints, the minimum filter orders are 104 and102, respectively, whereas the responses of the resulting filters are depicted in Fig. 10 using a dashed red line and a dot-dashed black line, respectively.

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

 Example Original Implementation Implementation II Reduction Percentage(%) Iterations Time* Iterations Time Iterations Time 15 34 0.223 7 0.024 79 89 26 35 0.269 16 0.040 54 85 3Τ 15 0.160 15 0.033 − 79 3⊥A 93 0.617 23 0.055 75 91 3⊥B 38 0.283 20 0.051 47 82 4Τ 11 0.146 11 0.029 − 80 4⊥A 49 0.360 37 0.082 24 77 4⊥B 66 0.398 23 0.057 65 86

### Table 2.

Performance Comparison of Original Implementation and Implementation II.

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.

## References

1. 1. AhsanM.2008Design of optimum linear phase FIR filters with a modified implementation of the Remez multiple exchange algorithm, In:Department of Signal Processing, Tampere University of Technology, Master’s thesis, (Sep. 2008), Tampere, Finland, 107 pages.
2. 2. AhsanM.SaramäkiT.2009Significant improvements in translating the Parks-McClellan algorithm from its FORTRAN code to its corresponding MATLAB codeIn:IEEE Symp. Circuits Syst., (May 2009), Taipei, Taiwan, 289292
3. 3. AhsanM.SaramäkiT.2011A MATLAB Based Optimum Multiband FIR Filters Design Program Following the Original Idea of the Remez Multiple Exchange AlgorithmIn: IEEE Symp. Circuits Syst., (May 2011), Rio de Janiero Brazil, 137140
4. 4. CheneyW. E.1966Introduction to Approximation TheoryAMS Chalsea Publishing.
5. 5. Mc ClellanJ. H.ParksT. W.1972A program for the design of linear phase finite impulse response digital filtersIn: IEEE Trans. Audio Electroacoust., AU-20Aug. 1972) 195199
6. 6. Mc ClellanJ. H.ParksT. W.1973A unified approach to the Design of Optimum FIR linear phase digital filters, In: IEEE Trans. Circuit Theory, 20Nov. 1973) 697701
7. 7. Mc ClellanJ. H.ParksT. W.RabinerL. R.1973A computer program for designing optimum FIR linear phase digital filtersIn: IEEE Trans. Audio Electroacoust., 21Dec. 1973) 506526
8. 8. NovodvorskiiE. P.PinskerI. S.1951The process of equating maxima, In: Uspekhi Mat. Nauk (USSR), 6174181Engl. transl. by A. Schenitzer).
9. 9. RabinerL. R.KaiserJ. F.SchaferR. W.1974Some considerations in the design of multiband finite-impulse-response digital filtersIn: IEEE Trans. Acoust. Speech, Signal Processing, 22Dec. 1974) 462472
10. 10. RabinerR. L.GoldG.1975Theory and Application of DIGITAL SIGNAL PROCESSINGPrentice Hall.
11. 11. RemezE.1934Sur le calcul effectifdes polynomes d‘approximation de Tchebychef, In: Compt. Rend. Acad. Sci, 199337340
12. 12. RiceR. J.1964The Approximation of Functions1Linear Theory, Addison-Wesley Pub (Sd).
13. 13. RivlinJ. T.2010An Introduction to the Approximation of FunctionsDover Publications.
14. 14. SaramäkiT.1981Design of digital filters requiring a small number of arithmetic operationsIn:Dept. of Electrical Engineering, Tampere University of Technology, Dr. Tech. Dissertation, Publ. 12, Tampere, Finland, 1981,226 pages.
15. 15. SaramäkiT.1987Efficient iterative algorithms for the design of optimum IIR filters with arbitrary specifications, In:Proc. Int. Conf. Digital Signal Process., Florence,Italy, (Sep. 1987) 3236
16. 16. SaramäkiT.1992An efficient Remez-type algorithm for the design of optimum IIR filters with arbitrary partially constrained specifications, In: IEEE Symp. Circuits Syst., 5May 1992) San Diego CA, 25772580
17. 17. SaramäkiT.1993Finite impulse response filter design, In: Handbook for Digital Signal Processing, Mitra, S. K. & Kaiser, J. F., (Eds.), Ch. 4, (1993) New York, NY: JohnWiley and Sons, 155277
18. 18. SaramäkiT.1994Generalizations of classical recursive digital filters and their design with the aid of a Remez-type Algorithm, In: IEEE Symp. Circuits Syst., 2May 1994), London UK, 549552
19. 19. SaramäkiT.RenforsM.1995A Remez-type algorithm for designing digital filters composed of all-pass sections based on phase approximationsIn: Proc. 38th Midwest Symp. Circuits Syst., 1Aug. 1995), Rio de Janiero Brazil, 571575
20. 20. TemesG. C.CalahanD. A.1967Computer-aided network optimization the state-of-the-artIn:Proc. IEEE, 55Nov. 1967), 18321863Nov. 1967.
21. 21. The MathWorks Inc.2009Filter design toolbox user’s guide, In: MATLAB Product Help, Version 4.6, (Sep. 2009), The MathWorks Inc., Natick, MA.

## Notes

• In this contribution, nz−1is chosen to be the number of adjustable parameters in both the FORTRAN and the MATLAB implementations of the PM algorithm because in this case nz stands for the number of extrema at which Alternation Theorem should be satisfied in order to guarantee the optimality of the solution.
• In the original PM algorithm, this range is the baseband for the so-called normalized frequencies from which the corresponding angular frequencies are obtained by multiplying these frequencies by 2π [17].,
• It should be noted that the value of dev is either positive or negative
• is worth emphasizing that implicit solutions for the calculation ofa[n]’s are not required for the intermediate iterations. The explicit solution for the calculation ofa[n]’s is needed only after achieving the convergence to the best solution.
• Vicinity Search automatically guarantees that the sign of the weighted error function wei_err(k) is same at both k=ℓtrial(1) andk=ℓreal(1). Hence, uppendis the smallest value ofk, where the sign of wei_err(k) can be opposite before these values. Similarly, the sign of wei_err(k) is same at both k=ℓtrial(nz) and k=ℓreal(nz) and lowend as given by (18) is the largest value ofk, where the sign of wei_err(k) can be opposite after these values.
• The only condition for α is that it should be small enough to avoid the extra peaks between the adjacent passbands and the newly formed intervals in the transition band regions.
• It is worth emphasizing that the use of max{Dκ+δκ,Dκ+1+δκ+1} implies that the maximum allowable value in the nearest passband is the upper limit. Similarly, in the following equation, min{Dκ−δκ,Dκ+1−δκ+1}means that the lower limit is the one in the nearest stopband.

Written By

Muhammad Ahsan and Tapio Saramäki

Submitted: December 8th, 2011 Published: September 26th, 2012