**Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City**

Farhad E. Mahmood

18 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

quality and standard conforming.

software development process.

*Universidad de Valladolid, Spain* 

**Acknowledgement** 

**6. References** 

Mariano Raboso and Myriam Codes

*Universidad Pontificia de Salamanca (Facultad de Informática), Spain* 

Available from: http://www.vex.net/~cthuang/tcom/

Vol.17, No.1, (January 2009), pp. 122-135, ISSN 0718-3305

Available from: http://www.tcl.tk/doc/

*Documentation*, 23.03.2012, Available from:

María I. Jiménez, Lara del Val, Alberto Izquierdo and Juan J. Villacorta

**Author details** 

Engineers working with Matlab and other software can take advantage of CBSE using COM and .Net technologies from Microsoft. Furthermore, integrating different objects from different applications accelerates software development and reduces costs. Today, COTS (Commercial Off-The-Shelf) components are ready to be integrated into the engineering applications and vendors exploit these advantages developing components while assuring

After reading this chapter, I expect that the interested reader can take into account component technology on their future projects, and gain effectiveness on the overall

This research has been supported by projects: 10MLA-IN-S08EI-1 (Pontifical University of Salamanca), and PON323B11-2 (Junta de Castilla y León). I also want to thank the work made by professors Maribel and Domingo, who carefully made the technical review of the chapter.

Gunderloy, M. (2001). Calling COM Components from .NET clients, In: *MSDN Library*, 23.03.2012, Available from: http://msdn2.microsoft.com/en-us/library/ms973800.aspx Huang, C. (2006). Tcom, In: *Access and implement Windows COM objects with Tcl*, 23.03.2012,

Raboso, M.; Izquierdo, A. & Villacorta J.J. (2003). Beamforming Systems Modeling using Component Reusability with XML Language, Proceedings of the International Signal

Tcl community (2009). Tcl and Tk manual pages, In: *Tcl/Tk Documentation*, 23.03.2012,

The Mathworks (2012). MATLAB COM Automation Server Support, In: *Matlab R2012a* 

http://www.mathworks.es/help/techdoc/matlab\_external/brd0v3w.html

Processing Conference, ISPC 2003, Dallas, Texas, USA, March 31-April 3, 2003. Raboso, M. (2007). *Beamforming Systems Modeling Using XML Language, based on Software Component Reuse,* ProQuest Information and Learning, ISBN 978-0-549-26134-6, USA Raboso M. ; Izquierdo A. ; Villacorta J. ; del Val L. & Jiménez, M. (2009). Integración de componentes COM de MATLAB/SIMULINK en el entorno CASE XBDK para el modelado de sistemas de conformación de haz. *Ingeniare, Revista chilena de ingeniería*, Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46469

## **1. Introduction**

Since the nineties, a great progress observed in wireless and mobile communications field. Mobiles make many tasks simple and easily to be accomplished. This requires imperatively good and accurate knowing and describing radio propagation for the intended area. Propagation of mobile radio signals is complex, and depend on the applied environment (Theodore S Rappaport, 2002).

 In this chapter, two theoretical models considered for the prediction of path loss in two different districts at Mosul city, using MATLAB 7.4 program. The propagation of radio wave of built-up area strongly influenced by the nature of environment, Parameters like antennas height, frequency band, interference, polarization etc, have an influence upon channel behavior. Reflection, diffraction and scattering from building add multipath, fast fading and slow fading effects to the signal.

The Walfisch-Ikegami (W-I) model used for uniform heights and similar buildings in the Karama district to calculate the path loss. The other model used is Okumura-Hata (OH) model applied for irregular and dissimilar buildings (M.O. Kabaou 2007) in the Almajmoa'a district. The information buildings heights for two areas obtained from the Civil Engineering Department, in the University of Mosul.

## **2. Cellular concept**

The key idea of mobile communication system is to offer possibility for a communication link in the serving area of the mobile network regardless of the location of the user. These demand two requirements for service providers (operators) of mobile communication system. Firstly, the service area has to be geometrically large area in order to provide services for a large number of costumers. Secondly, the network has to be able to guarantee

mobility, i.e., customers have to be able to move even considerably long distances without breaking already establish connections. Achieving this requirement is technically challenging, and more over there is limited amount of air interface frequency band available (T. Rappaport 2002).

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 21

suburban areas comprise residential houses, parks and gardens. The term 'rural' defines open farmland with sparse buildings, forest and woodland. When the effects of environment have been conceder, six factors are useful in classifying land usage (J.D.Parson

Mobile radio channels needs a rigorous vision of channel characterization; the study of channel will be useful for mobile wireless communication to make the best decision between source and destination, and even advising recommendations and directives for dimensioning exploited network (antenna heights, carrier, network architecture, etc...).

Propagation channel is a linear system defining a transformation between an input and an output signal, for best understanding the channel and overcome its eventual imperfections. In fact, signal affronts many distortions like delays, spreading. These distortions induced by various reflections that signal faces up during the Emission- Reception path. Consequently, other additive signals will be perceived by the receiver besides the main transmitted signal which must be captured in ideal situations. These additional signals have followed various and different paths. That is what usually named: Multipath. With higher rates used in digital communications e.g., symbols have small duration in comparison with the delay

Multipath is a phenomenon that happens in the channel of mobile systems when the transmitted signal arrives at the receiver via different paths due to reflection, diffraction and scattering resulting in fading. There is only one transmitted signal, due to obstacles like buildings, hills, trees, and so on, in the signal paths that cause different signals to arrive at the receiver from various directions with different delays (M. Rahnema 2008). In such situations, receiver will collect light of seen's (LOS) rays but also other rays having arrived over indirect ways. They know as Non-light of Sight (NLOS) rays. More further, are Emission and the Reception one side from the other, more lower will be chosen the used frequency. Possible paths will consequently be higher. Propagation channel causes then Inter Symbol Interferences (ISI). For that, the parameter conditioning propagation delay spread values must imperatively specify. Fig. 2. illustrates three propagation

spreading scale and a really superposition between them will be observed.

Building density (percentage of area covered by building).

Building size (area covered by building)

**4. Radio mobile channel characterization** 

2000):

Building height.

 Building location. Vegetation location. Terrain density.

**4.1. Multipath propagation** 

components.

One solution to satisfy these requirements is utilization of *cellular concept* (Lioh. Wacker. 2006). in cellular concept, the basic approach is to divide a large service area into smaller sub areas called cells. The main advantage of cellular solution is increasing capacity. Since the limited radio frequency, range used only in small areas. The same radio frequency can utilized again after a certain physical distance. Once other advantage is that, the base station antenna heights can be much lower due to smaller serving areas. This enables the usage of lower transmit powers, which also saves the battery of the mobile. Furthermore, change capacity and coverage demands for different area can be easily adopted using smaller cells where the traffic density is higher.

In Fig.1, cells depicted with a shape of hexagon. In practice, the concept of hexagonal cells is conceptual due to non-homogenous environment. The morphology of the terrain introduces different propagation environments (e.g., open areas, forests, and waters) and topography fluctuations of the terrain. If the radio wave attenuation is enough, or the distances between the sites are long enough, it is possible to reuse frequencies after a certain distance. The frequencies used in all the grey areas of the figure could be the same, as well as those in the white and the lined areas, respectively, this increase the capacity of cellular communication system. The cellular concept introduces also some disadvantages. Due to the large service area covered by multiple cells, the system has to know the location of the users in order to, e.g., route an incoming call. Thus, some signaling capacity needed in order to control the communications. Another problem arises when continuous service provided for the mobiles on move. To solve these problems, the network management system becomes more complicated, expensive, and more difficult to handle.

**Figure 1.** Cellular Concept

## **3. An environment classification for a mobile radio channel**

The propagation of radio wave of built-up district strongly influenced by the nature of environment, in particular the density and size of buildings. In propagation studies for mobile radio, a qualitative description of the environment area is often employed using term such as dense urban, urban, sub urban and rural. 'Dense urban' area generally defined as being dominated by tall buildings, office blocks and other commercial buildings. Whereas suburban areas comprise residential houses, parks and gardens. The term 'rural' defines open farmland with sparse buildings, forest and woodland. When the effects of environment have been conceder, six factors are useful in classifying land usage (J.D.Parson 2000):


20 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

(T. Rappaport 2002).

where the traffic density is higher.

**Figure 1.** Cellular Concept

complicated, expensive, and more difficult to handle.

**3. An environment classification for a mobile radio channel** 

The propagation of radio wave of built-up district strongly influenced by the nature of environment, in particular the density and size of buildings. In propagation studies for mobile radio, a qualitative description of the environment area is often employed using term such as dense urban, urban, sub urban and rural. 'Dense urban' area generally defined as being dominated by tall buildings, office blocks and other commercial buildings. Whereas

mobility, i.e., customers have to be able to move even considerably long distances without breaking already establish connections. Achieving this requirement is technically challenging, and more over there is limited amount of air interface frequency band available

One solution to satisfy these requirements is utilization of *cellular concept* (Lioh. Wacker. 2006). in cellular concept, the basic approach is to divide a large service area into smaller sub areas called cells. The main advantage of cellular solution is increasing capacity. Since the limited radio frequency, range used only in small areas. The same radio frequency can utilized again after a certain physical distance. Once other advantage is that, the base station antenna heights can be much lower due to smaller serving areas. This enables the usage of lower transmit powers, which also saves the battery of the mobile. Furthermore, change capacity and coverage demands for different area can be easily adopted using smaller cells

In Fig.1, cells depicted with a shape of hexagon. In practice, the concept of hexagonal cells is conceptual due to non-homogenous environment. The morphology of the terrain introduces different propagation environments (e.g., open areas, forests, and waters) and topography fluctuations of the terrain. If the radio wave attenuation is enough, or the distances between the sites are long enough, it is possible to reuse frequencies after a certain distance. The frequencies used in all the grey areas of the figure could be the same, as well as those in the white and the lined areas, respectively, this increase the capacity of cellular communication system. The cellular concept introduces also some disadvantages. Due to the large service area covered by multiple cells, the system has to know the location of the users in order to, e.g., route an incoming call. Thus, some signaling capacity needed in order to control the communications. Another problem arises when continuous service provided for the mobiles on move. To solve these problems, the network management system becomes more


Mobile radio channels needs a rigorous vision of channel characterization; the study of channel will be useful for mobile wireless communication to make the best decision between source and destination, and even advising recommendations and directives for dimensioning exploited network (antenna heights, carrier, network architecture, etc...).

## **4. Radio mobile channel characterization**

Propagation channel is a linear system defining a transformation between an input and an output signal, for best understanding the channel and overcome its eventual imperfections. In fact, signal affronts many distortions like delays, spreading. These distortions induced by various reflections that signal faces up during the Emission- Reception path. Consequently, other additive signals will be perceived by the receiver besides the main transmitted signal which must be captured in ideal situations. These additional signals have followed various and different paths. That is what usually named: Multipath. With higher rates used in digital communications e.g., symbols have small duration in comparison with the delay spreading scale and a really superposition between them will be observed.

## **4.1. Multipath propagation**

Multipath is a phenomenon that happens in the channel of mobile systems when the transmitted signal arrives at the receiver via different paths due to reflection, diffraction and scattering resulting in fading. There is only one transmitted signal, due to obstacles like buildings, hills, trees, and so on, in the signal paths that cause different signals to arrive at the receiver from various directions with different delays (M. Rahnema 2008). In such situations, receiver will collect light of seen's (LOS) rays but also other rays having arrived over indirect ways. They know as Non-light of Sight (NLOS) rays. More further, are Emission and the Reception one side from the other, more lower will be chosen the used frequency. Possible paths will consequently be higher. Propagation channel causes then Inter Symbol Interferences (ISI). For that, the parameter conditioning propagation delay spread values must imperatively specify. Fig. 2. illustrates three propagation components.

**Figure 2.** Propagation Component between Base Station and Mobile Station

Where (L. Poole 2006):


Signals coming from the same source will arrive to at the receiver side affected by different delays. Delay spread evaluated as (S. Tabbane 2000):

$$\text{Multi path spread} = \frac{\text{longest path} - \text{shortest path}}{c} \tag{1}$$

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 23

**Figure 3.** Channel characterization versus different vision

**4.2. Slow and fast fading** 

also called log-normal fading.

**4.3. Rayleigh fading** 

Table 1 for different areas and situations (J.D.Parson 2000).

Direction and distance between reflectors (buildings, mountains, walls, cars) influence channel awkward and followed channel. The IR of a dispersive delay channel depends then on number of physical factors. The delay spread value varies from ten nanoseconds for indoor propagation up to some microseconds in outdoor propagation case as indicated in

In addition to distance dependent path loss, the received signal level experiences also fluctuations that called slow and fast fading. When the receiver placed in, a coverage area that either omits the LOS component or other dominant component, the receiver is considered shadowed. Usually these shadowing obstacles are big trees, buildings, or in rural environment the hills. Slow fading is shown to follow log-normal distribution, and thus it is

Fast fading is a consequence of the multipath propagation phenomena. As described earlier, the multipath propagated signal components can received either constructively or deconstructively depending on the relative phase difference between the components. This causes very rapid fluctuations in the received signal level. Slow and fast fading components

After reflection, the phases and magnitude of the radio wave might be attaint attenuation or fading phenomena cause essentially temporal variations in phases. These come from additive multiple signals In this case; the resulting received signal in the receiver could be null or a very weak signal. Multiple received signals could also add in a constructive manner and the resulting received signal will have a magnitude greater than signal issued along directly path. Taking account of propagation conditions's variation between kinds of

of received signal illustrated in Fig.4. in respect to receiver-transmitter distance.

domains, we must note fading depends on the following (J. S. Seybold 2005):

Where c commonly known light celerity. Table below shows and gives better and more accurate idea.


**Table 1.** Comparison of delay spread in different environments.

We can particularly note that delays caused indoor are fewer than other referred to abroad situations (some tens of s) (M. J. Nawrocki 2006). For that, a suitable modeling must accomplished for each situation often named Indoor/Outdoor. Therefore, the channel has many dimensions, giving a particular description as indicated in the fig. 3. (V. Erceg 2005)

Parameters like antennas heights, frequency band, interference, polarization etc, and have an influence upon channel behavior indicated above in the fig. 3.

**Figure 3.** Channel characterization versus different vision

Direction and distance between reflectors (buildings, mountains, walls, cars) influence channel awkward and followed channel. The IR of a dispersive delay channel depends then on number of physical factors. The delay spread value varies from ten nanoseconds for indoor propagation up to some microseconds in outdoor propagation case as indicated in Table 1 for different areas and situations (J.D.Parson 2000).

### **4.2. Slow and fast fading**

� (1)

22 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 2.** Propagation Component between Base Station and Mobile Station

Fast fading: due to the multipath reflection of the transmitted signal.

Slow fading (Shadowing): is due to obstruction of the signal by natural obstacles.

Signals coming from the same source will arrive to at the receiver side affected by different

Multi path spread = ������� ���� � �������� ����

Where c commonly known light celerity. Table below shows and gives better and more

Environment type Delay spread (s)

We can particularly note that delays caused indoor are fewer than other referred to abroad situations (some tens of s) (M. J. Nawrocki 2006). For that, a suitable modeling must accomplished for each situation often named Indoor/Outdoor. Therefore, the channel has many dimensions, giving a particular description as indicated in the fig. 3. (V. Erceg 2005)

Parameters like antennas heights, frequency band, interference, polarization etc, and have

Path loss: decay of the signal strength with distance.

delays. Delay spread evaluated as (S. Tabbane 2000):

Free space < 0.2 Rural area 1

Mountainous area 30 to 50 Suburban area 0.5 Urban area 3 Indoor area 0.1

**Table 1.** Comparison of delay spread in different environments.

an influence upon channel behavior indicated above in the fig. 3.

Where (L. Poole 2006):

accurate idea.

In addition to distance dependent path loss, the received signal level experiences also fluctuations that called slow and fast fading. When the receiver placed in, a coverage area that either omits the LOS component or other dominant component, the receiver is considered shadowed. Usually these shadowing obstacles are big trees, buildings, or in rural environment the hills. Slow fading is shown to follow log-normal distribution, and thus it is also called log-normal fading.

Fast fading is a consequence of the multipath propagation phenomena. As described earlier, the multipath propagated signal components can received either constructively or deconstructively depending on the relative phase difference between the components. This causes very rapid fluctuations in the received signal level. Slow and fast fading components of received signal illustrated in Fig.4. in respect to receiver-transmitter distance.

## **4.3. Rayleigh fading**

After reflection, the phases and magnitude of the radio wave might be attaint attenuation or fading phenomena cause essentially temporal variations in phases. These come from additive multiple signals In this case; the resulting received signal in the receiver could be null or a very weak signal. Multiple received signals could also add in a constructive manner and the resulting received signal will have a magnitude greater than signal issued along directly path. Taking account of propagation conditions's variation between kinds of domains, we must note fading depends on the following (J. S. Seybold 2005):

	- Transmitted Signal's bandwidth,
	- Delay spread of the received signal,
	- Random phase and magnitude of Multipath components,
	- Transmitter, receiver and around existing objects position cause CIR's temporal fluctuations.

**Figure 4.** Slow and fast fading of propagated signal in respect to receiver-transmitter distance

### **4.4. Doppler shift**

This phenomena result from Displacement and mobility of the MS in relation to BS. This introduces frequency shift of the received signal. Frequency shift depends essentially on:


If we note:


The frequency received by a receiver moving at elative speed ѵ relatively to the emitter is given by the following (F.I . Mahmoud 2008).

$$f' = f - \frac{\upsilon}{\lambda} \tag{2}$$

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 25

���� = (4����)�2 (3)

are several situations in which it can be used and it also useful as the basis for

The free space path loss formula or free space path loss equation is quite simple to use. Not only is the path loss proportional to the square of the distance between the transmitter and receiver, but the signal level is also proportional to the square of the frequency

due to multipath effect and fading it cannot depend on free space formula to find path loss

Okumura-Hata Model is an Empirical model; which based on a Statistical Analysis of a great number of experimental measures that takes in account many parameters such as buildings, Base stations (BS) and Mobile stations heights (S.A.Mawjoud 2008). The model based on the field measurements made by Okumura in Tokyo in 1968. Okumura's measurements were later fitted in a mathematical formula by Hata . The Okumura-Hata model is valid for frequency range of 150–1500 MHz, the BS height is 30-200m, the MS height is 1-10m and the distance is 1- 20km. Because of the frequency band limitation of the Okumura-Hata model, the original model was later extended by the European CO-operation in the field of Scientific and Technical research (COST 231) to the frequency band of 1500–2000 MHz , which covers the

� (�����)(��) = 46.3 + 33.9 ���(����) − 13.82 log(ℎ�) − �(ℎ�) + (44.9 − 6.55 log(ℎ�)) log(���) (4)

bands allocated to the 3G networks (E. Damasso 1999), is given in equation (4)

Methods to model the mobile radio propagation channel are (M.O. Kabaou 2007):

understanding many real life radio propagation situations.

d is the distance of the receiver from the transmitter (metres)

c is the speed of light in a vacuum (metres per second)

**5.1. Okumura-Hata Model (OH), statistical models** 

there are empirical model used to find path loss.

 Statistical models ,like (Okumura-Hata) , Physical Models ,like (Walfisch-Ikegami),

In addition, References Models.

We study only the first two cases.

**Calculation of free-space path loss:** 

FSPL is the Free space path loss

λ is the signal wavelength (metres) f is the signal frequency (Hertz)

(T. Rappaport, 2002).

Mixed Models,

Where:

Doppler shift causes a random frequency modulation for the signal, and could introduce Multipath Signals. Positive and negative shifts both introduced.

### **5. Path loss models**

#### **Free space:**

The free space path loss used in many areas for predicting radio signal strengths that may expected in a radio system. Although it does not hold for most terrestrial situations as there are several situations in which it can be used and it also useful as the basis for understanding many real life radio propagation situations.

#### **Calculation of free-space path loss:**

The free space path loss formula or free space path loss equation is quite simple to use. Not only is the path loss proportional to the square of the distance between the transmitter and receiver, but the signal level is also proportional to the square of the frequency (T. Rappaport, 2002).

$$FSPL = \langle 4\pi d/\lambda \rangle^{\wedge} 2 \tag{3}$$

Where:

24 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

Transmitter, receiver and around existing objects position cause CIR's temporal

**Figure 4.** Slow and fast fading of propagated signal in respect to receiver-transmitter distance

This phenomena result from Displacement and mobility of the MS in relation to BS. This introduces frequency shift of the received signal. Frequency shift depends essentially on:

The frequency received by a receiver moving at elative speed ѵ relatively to the emitter is

�� ��� <sup>ѵ</sup>

Doppler shift causes a random frequency modulation for the signal, and could introduce

The free space path loss used in many areas for predicting radio signal strengths that may expected in a radio system. Although it does not hold for most terrestrial situations as there

�

(2)

Random phase and magnitude of Multipath components,

 Transmitted Signal's bandwidth, Delay spread of the received signal,

fluctuations.

**4.4. Doppler shift** 

Mobile position,

*λ* is the wavelength¸

**5. Path loss models** 

**Free space:** 

*f* the carrier used frequency.

given by the following (F.I . Mahmoud 2008).

If we note:

The speed of the receiver according to the emitter.

Multipath Signals. Positive and negative shifts both introduced.

FSPL is the Free space path loss

d is the distance of the receiver from the transmitter (metres) λ is the signal wavelength (metres) f is the signal frequency (Hertz) c is the speed of light in a vacuum (metres per second)

due to multipath effect and fading it cannot depend on free space formula to find path loss there are empirical model used to find path loss.

Methods to model the mobile radio propagation channel are (M.O. Kabaou 2007):


We study only the first two cases.

#### **5.1. Okumura-Hata Model (OH), statistical models**

Okumura-Hata Model is an Empirical model; which based on a Statistical Analysis of a great number of experimental measures that takes in account many parameters such as buildings, Base stations (BS) and Mobile stations heights (S.A.Mawjoud 2008). The model based on the field measurements made by Okumura in Tokyo in 1968. Okumura's measurements were later fitted in a mathematical formula by Hata . The Okumura-Hata model is valid for frequency range of 150–1500 MHz, the BS height is 30-200m, the MS height is 1-10m and the distance is 1- 20km. Because of the frequency band limitation of the Okumura-Hata model, the original model was later extended by the European CO-operation in the field of Scientific and Technical research (COST 231) to the frequency band of 1500–2000 MHz , which covers the bands allocated to the 3G networks (E. Damasso 1999), is given in equation (4)

� (�����)(��) = 46.3 + 33.9 ���(����) − 13.82 log(ℎ�) − �(ℎ�) + (44.9 − 6.55 log(ℎ�)) log(���) (4)

Where:

���� is the operating frequency in MHz (1500 MHz -2000 MHz).

ℎ� the base station(BS) antenna height in meter (30m- 200m).

ℎ� the mobile(MS) antenna height in meter (1m – 10m) .

��� the distance between BS and MS in km.

�(ℎ�) the correction factor for mobile antenna height which is given by:

$$a(h\_m)\_{dB} = (1.1\log(f\_{MHz}) - 0.7)h\_m - 1.56\log(f\_{MHz}) - 0.8\tag{5}$$

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 27

*L=LO +Lrts+Lmsd* (6)

�� � ���� � �� ��� ��� � ���������� (7)

(9)

���� � −���� − �� ��� � � ��������� � �� ����ℎ���� − ℎ��−���� (8)

10 0.354 for 0 35

2.5 0.075 35 for 35 <55 deg

4.0-0.114 55 for 55 <90

deg

deg

4. The orientation of the streets relative to the line of sight (LOS) and non line of sight (NLOS). 5. The distance between the Tx. and receiver Rx, the heights of Tx. and Rx. antennas and

COST 231 Walfisch–Ikegami (WI) model is the most widely used empirical model today, being an extension of the models from J. Walfisch and F. Ikegami. To work up to 2GHz by The European comparative for scientific and technical, research (EURO-COST). It has been adopted as a standard model for 3G IMT 2000/UMTS systems. It is valid within the following constraints:

3. The distance between buildings.

the frequency of operation.

 BS height: 4–50m. MS height: 1–3m.

Where:

Where:

*LO* : free space loss.

BS to MS separation : 0.02–5 km.

ℎ����: The height of the buildings.

The Walfisch – Ikegami equation is given by:

*Lrts* : rooftop to street diffraction and scatter loss.

*ori L*

*Lmsd* : multi screen diffraction loss due to the rows of building .

**Figure 5.** Propagation geometry for Walfisch– Ikegami model

Where: 1 m ≤ hm ≤ 10 m

Equation (3) is an extension version of original Hata model to work up to 2GHz by The European comparative for scientific and technical research (EURO-COST) (M. J. Nawrocki 2006).

The BS antenna height must be above the rooftop level of the buildings adjacent to the BS. Thus, the model is proposed to be used in propagation studies of macro-cells. There are several weaknesses in the empirical or semi-empirical models for propagation studies in micro-cellular environments. If the BS antenna height is below the level of the rooftop of surrounding buildings, the nature of the propagation phenomena will changes. This situation cannot be analyzed with statistical methods because the individual buildings are too large compared with the cell size and the exact geometrical properties of the buildings can no longer be ignored as they can in macro cellular models.

### **5.2***.* **Walfisch-Ikegami Model (W-I) physical models**

This model based on the assumption that the transmitted wave propagates over the rooftops by a process of multiple diffraction in regular area, regular area has rows of buildings, which are nearly of equal and uniform height and are located on flat terrain. The building are organized along the street grid with little or no side-to-side spacing and nearly equal front-to-front and back-to-back spacing. The propagation takes place primarily over building (T. Rappaport. 2002). The model considers the impact on rooftops and building heights using diffraction to predict average signal strength at street level. The models consider the path loss (L) to be the product of three factors (J.D.Parson 2000):


The geometry used in the Walfisch model shown in fig.5.

From fig. 5, it can be noted that the diffracted signals arrive at the receiver from the 1st, 2nd and 3rd as if the model show that the signals travels from rooftops.

This model depends on:


3. The distance between buildings.

26 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

���� is the operating frequency in MHz (1500 MHz -2000 MHz).

�(ℎ�) the correction factor for mobile antenna height which is given by:

�(ℎ�)�� = (1.1 log(���� ) − 0.7)ℎ� − 1.56 log(����) − 0.8 (5)

Equation (3) is an extension version of original Hata model to work up to 2GHz by The European

The BS antenna height must be above the rooftop level of the buildings adjacent to the BS. Thus, the model is proposed to be used in propagation studies of macro-cells. There are several weaknesses in the empirical or semi-empirical models for propagation studies in micro-cellular environments. If the BS antenna height is below the level of the rooftop of surrounding buildings, the nature of the propagation phenomena will changes. This situation cannot be analyzed with statistical methods because the individual buildings are too large compared with the cell size and the exact geometrical properties of the buildings

This model based on the assumption that the transmitted wave propagates over the rooftops by a process of multiple diffraction in regular area, regular area has rows of buildings, which are nearly of equal and uniform height and are located on flat terrain. The building are organized along the street grid with little or no side-to-side spacing and nearly equal front-to-front and back-to-back spacing. The propagation takes place primarily over building (T. Rappaport. 2002). The model considers the impact on rooftops and building heights using diffraction to predict average signal strength at street level. The models

From fig. 5, it can be noted that the diffracted signals arrive at the receiver from the 1st, 2nd

consider the path loss (L) to be the product of three factors (J.D.Parson 2000):

comparative for scientific and technical research (EURO-COST) (M. J. Nawrocki 2006).

ℎ� the base station(BS) antenna height in meter (30m- 200m). ℎ� the mobile(MS) antenna height in meter (1m – 10m) .

can no longer be ignored as they can in macro cellular models.

**5.2***.* **Walfisch-Ikegami Model (W-I) physical models**

 Losses introduced by rooftops near the studied area. The geometry used in the Walfisch model shown in fig.5.

and 3rd as if the model show that the signals travels from rooftops.

2. The width of the streets and the width of the buildings.

��� the distance between BS and MS in km.

Where: 1 m ≤ hm ≤ 10 m

Free space loss.

This model depends on:

1. The height of the buildings.

Losses added by diffraction.

Where:


**Figure 5.** Propagation geometry for Walfisch– Ikegami model

COST 231 Walfisch–Ikegami (WI) model is the most widely used empirical model today, being an extension of the models from J. Walfisch and F. Ikegami. To work up to 2GHz by The European comparative for scientific and technical, research (EURO-COST). It has been adopted as a standard model for 3G IMT 2000/UMTS systems. It is valid within the following constraints:


The Walfisch – Ikegami equation is given by:

$$L = L\_{\cdot} + L\_{\cdot rts} + L\_{\cdot msd} \tag{6}$$

Where:

*LO* : free space loss.

*Lrts* : rooftop to street diffraction and scatter loss.

*Lmsd* : multi screen diffraction loss due to the rows of building .

$$L\_o = 32.4 + 20\log d\_{km} + 20\log f\_{\text{MHz}} \tag{7}$$

$$L\_{\rm rts} = -16.9 - 10 \log w + 10 \log f\_{\rm MHz} + 20 \log \left( h\_{\rm rnof} - h\_m \right) - L\_{\rm orl} \tag{8}$$

Where:

ℎ����: The height of the buildings.

$$L\_{ori} = \begin{cases} -10 + 0.354 \frac{\text{\(\rho\)}}{\text{deg}} & \text{for } 0^\circ \le \rho < 35^\circ\\ 2.5 + 0.075 \left(\frac{\rho}{\text{deg}} - 35\right) \text{for } 35^\circ \le \rho < 55^\circ\\ 4.0 \text{-} 0.114 \left(\frac{\rho}{\text{deg}} - 55 & \text{for } 0^\circ 55 \le \rho < 90^\circ \end{cases} \right) \tag{9}$$

is the angle between incident wave and street oriented.

$$L\_{msd} = L\_{bsh} + k\_a + k\_d \log(d\_{Km}) + k\_f \log(f\_{MHz}) - 9 \log(b) \tag{10}$$

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 29

Fig.9. Shows the variation of path loss with distance, which shows that as the coverage of cell increases the path loss increases, using equation 4 & 5 for operation frequency of 2000MHz, 1.5m mobile antenna height and 25m base station antenna height. That is mean

when mobile move far from base station the received signal decreases gradually.

i. Study of coverage range in OH model:

**Figure 6.** Mosul city map by Google Earth

**Figure 7.** Almajmoa'a region in Mosul city.

*Where b:* The distance between buildings.

$$L\_{bsh} = \begin{cases} -18\log\left[1 + \left(h\_b - h\_{R\alpha f}\right)\right] & \text{for } |h\_b > h\_{R\alpha f}| \\ 0 & \text{for } |h\_b \le h\_{R\alpha f}| \end{cases} \tag{11}$$

$$k\_a = \begin{cases} 54 & \text{for } h\_b > h\_{\text{Roff}} \\ 54 - 0.8 \left( h\_b - h\_{\text{Roff}} \right) & \text{for } \text{d} \ge 0.5 \text{km and } h\_b \le h\_{\text{Roff}} \\ 54 - 0.8 \left( h\_b - h\_{\text{Roff}} \frac{d\_{\text{km}}}{0.5} \right) & \text{for } \text{d} < 0.5 \text{km and } h\_b \le h\_{\text{Roff}} \end{cases} \tag{12}$$

$$k\_d = \begin{cases} 18 & \text{for } h\_b > h\_{\text{Rowf}} \\ 18 - 15 \frac{\left(h\_b - h\_{\text{Rowf}}\right)}{h\_{\text{Rowf}}} & \text{for } \quad h\_b \le h\_{\text{Rowf}} \end{cases} \tag{13}$$

$$k\_f = -4 + 1.5\left(\frac{f\_{MHz}}{925} - 1\right) \tag{14}$$

The term *<sup>a</sup> k* represents the increase in path loss when the base station antenna is below rooftop height. *<sup>d</sup> k* and *<sup>f</sup> k* allow for the dependence of the diffraction loss on range and frequency, respectively .

#### **6. Radio network planning**

Mosul city have large variation in topography there are terrain and water, shown in fig.6, two regions have been investigated, for terrain there is nonhomogenous area like Almajmoa'a region, also there is the Karama region that have regular area with equal height for all building. Usually the OH model used by mobile companies in all Mosul city because it is assumed irregular region, but regular region can be seen in Karama district (building structure is nearly equal), W-I propagation model is applicable to district of uniform area buildings, which is appropriate for Karama region.

#### **6.1. Path loss study for Almajmoa'a region using Okumura-Hata Model**

In Almajmoa'a region, shown in fig.7. , all buildings are irregular and not of uniform shape, for these reasons W-I model cannot be applied but OH is more suitable.

Simulations for *OH Model* in Almajmoa'a:

### i. Study of coverage range in OH model:

28 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

18 1 for 0 for

 

54 for 54 0.8 for d>0.5km and

 

*k hh h h*

54 0.8 for d 0.5km and 0.5

*Roof*

*h*

*d b Roof*

*k h h*

*f*

18 15 for

4 1.5 ( 1) <sup>925</sup> *MHz*

The term *<sup>a</sup> k* represents the increase in path loss when the base station antenna is below rooftop height. *<sup>d</sup> k* and *<sup>f</sup> k* allow for the dependence of the diffraction loss on range and

Mosul city have large variation in topography there are terrain and water, shown in fig.6, two regions have been investigated, for terrain there is nonhomogenous area like Almajmoa'a region, also there is the Karama region that have regular area with equal height for all building. Usually the OH model used by mobile companies in all Mosul city because it is assumed irregular region, but regular region can be seen in Karama district (building structure is nearly equal), W-I propagation model is applicable to district of uniform area

In Almajmoa'a region, shown in fig.7. , all buildings are irregular and not of uniform shape,

**6.1. Path loss study for Almajmoa'a region using Okumura-Hata Model** 

for these reasons W-I model cannot be applied but OH is more suitable.

18 for

*km*

*a b Roof b Roof*

*Log h h h h <sup>L</sup>*

���� � ���� � �� � �� log(���) � �� log(����) � �log (�) (10)

*b Roof b Roof*

*b Roof b Roof*

*<sup>d</sup> h h h h*

*b Roof*

*b Roof*

*h h*

(11)

(13)

(12)

*h h*

*b Roof*

*h h*

*h h*

*b Roof*

*<sup>f</sup> <sup>k</sup>* (14)

is the angle between incident wave and street oriented.

*Where b:* The distance between buildings.

*bsh*

frequency, respectively .

**6. Radio network planning** 

buildings, which is appropriate for Karama region.

Simulations for *OH Model* in Almajmoa'a:

Fig.9. Shows the variation of path loss with distance, which shows that as the coverage of cell increases the path loss increases, using equation 4 & 5 for operation frequency of 2000MHz, 1.5m mobile antenna height and 25m base station antenna height. That is mean when mobile move far from base station the received signal decreases gradually.

**Figure 6.** Mosul city map by Google Earth

**Figure 7.** Almajmoa'a region in Mosul city.

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 31

**hb=30 hb=50 hb=70 hb=100**

The result in Fig.10., shows the variation of path loss with distance for four antennas height of base station (30m,50m,70m and 100m ), which shows that as the antenna height of base station increases, path loss decreases. Therefore, the coverage of cell will increases, which mean for BS height 30m the coverage of cell, could be 5km, where for 100m the coverage may arrive 35km. Therefore, inside the city where there are large number of building and MS, this demand large number of BS, if the coverage of these BSs is large, that mean the interference among these cells will increase, then the blocking probability will increase. Therefore, it must decrease BSs height inside cities to reduce this interference. However, for outside cities, where there are large area, small building and small number of MSs, it recommended to increase BSs height to increase coverage to cover as large area as can.

**relation between hb&pl**

**Figure 10.** Path loss versus distance for four BS heights in OH model

given from Civil Engineering Department in Mosul university).

Simulations for W-I in karama:

90

100

110

120

**loss in dB**

130

140

150

160

mobile unit at cell boarder.

a. Study of coverage range in W-I model:

**6.2. Path loss study for karama region using Walfisch-Ikegami (W-I) model** 

In karama region shown in fig. 8, buildings are regular and of uniform shape, for these reasons W-I model can be applied. (karama's information which is used in simulation, are

0 0.5 1 1.5 2 2.5 3

**distance in km**

Buildings in this region have three floors with 3m height for each floor, i.e. the whole building's height (hroof) is 9m. Base station's height (hb) is 25m. buildings distance (b) are 6m , street width (w) is 4m, operating frequency 2000MHz, mobile antenna's height (hm) is 1.5 m . Fig.11. shows the relation between path loss and distance in Karama district with W-I model, this figure show that as the coverage of cell increases the path loss increases in

ii. Study of BS height in OH model:

**Figure 9.** Path loss versus distance for O-H in Almajmoa'a region.

The result in Fig.10., shows the variation of path loss with distance for four antennas height of base station (30m,50m,70m and 100m ), which shows that as the antenna height of base station increases, path loss decreases. Therefore, the coverage of cell will increases, which mean for BS height 30m the coverage of cell, could be 5km, where for 100m the coverage may arrive 35km. Therefore, inside the city where there are large number of building and MS, this demand large number of BS, if the coverage of these BSs is large, that mean the interference among these cells will increase, then the blocking probability will increase. Therefore, it must decrease BSs height inside cities to reduce this interference. However, for outside cities, where there are large area, small building and small number of MSs, it recommended to increase BSs height to increase coverage to cover as large area as can.

**Figure 10.** Path loss versus distance for four BS heights in OH model

## **6.2. Path loss study for karama region using Walfisch-Ikegami (W-I) model**

In karama region shown in fig. 8, buildings are regular and of uniform shape, for these reasons W-I model can be applied. (karama's information which is used in simulation, are given from Civil Engineering Department in Mosul university).

Simulations for W-I in karama:

30 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 8.** Karama region in Mosul city.

100

110

120

130

**pathloss in dB**

140

150

160

ii. Study of BS height in OH model:

**Figure 9.** Path loss versus distance for O-H in Almajmoa'a region.

0 0.5 1 1.5 2 2.5 3

**pathloss with distance in Almajmoaa region**

**distance in km**

a. Study of coverage range in W-I model:

Buildings in this region have three floors with 3m height for each floor, i.e. the whole building's height (hroof) is 9m. Base station's height (hb) is 25m. buildings distance (b) are 6m , street width (w) is 4m, operating frequency 2000MHz, mobile antenna's height (hm) is 1.5 m . Fig.11. shows the relation between path loss and distance in Karama district with W-I model, this figure show that as the coverage of cell increases the path loss increases in mobile unit at cell boarder.

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 33

The BS height larger than rooftop case, it can be seen that there are difference by about 20 dB

However, for BS height less than rooftop the difference in bath loss increase gradually to be

The simulation results of the two different districts considered are shown in fig.13. (for same BS and MS antenna height and same operating frequency),which shows that Karama region using W-I model have larger path loss than Almajmoa'a region using OH model about 7 - 10dB, because of rooftop and multi-screen effects in karama region which makes more

**Distance effect on Almajmoaa and karama regions**

0 0.5 1 1.5 2 2.5 3

**karama district Almajmoaa district**

**distance in km**

from rooftop BS height, regardless of distance between BS and MS.

10 dB at 3km, where in 600 meter the difference could be zero.

**Figure 13.** Path loss versus distance in karama and Almajmoa'a regions

**path loss in dB**

**7. Comparison between O-H and W-I** 

diffraction for signal.

32 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 11.** Path loss versus distance for W-I in Karama

b. Study of BS height in W-I model :

Fig.12. shows the variation of path loss with distance for three BSs height (hb=25m>hroof , hb=9m= hroof & hb=5m<hroof) this figure shows that to reduce the path loss they must increase BS height to become higher than buildings height, to avoid reflection from buildings .

**Figure 12.** Path loss versus distance for different BS heights.

The BS height larger than rooftop case, it can be seen that there are difference by about 20 dB from rooftop BS height, regardless of distance between BS and MS.

However, for BS height less than rooftop the difference in bath loss increase gradually to be 10 dB at 3km, where in 600 meter the difference could be zero.

## **7. Comparison between O-H and W-I**

32 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

Fig.12. shows the variation of path loss with distance for three BSs height (hb=25m>hroof , hb=9m= hroof & hb=5m<hroof) this figure shows that to reduce the path loss they must increase

**pathloss vs. distance in three case[hb>hr,hb=hr,hb<hr]**

0 0.5 1 1.5 2 2.5 3

**hBS < hroof hBS = hroof hBS > hroof**

**distance in km**

0 0.5 1 1.5 2 2.5 3

**pathloss vs. distance in Karama region**

**distance in km**

BS height to become higher than buildings height, to avoid reflection from buildings .

**Figure 11.** Path loss versus distance for W-I in Karama

**Figure 12.** Path loss versus distance for different BS heights.

**pathloss in dB**

b. Study of BS height in W-I model :

**path loss in dB**

The simulation results of the two different districts considered are shown in fig.13. (for same BS and MS antenna height and same operating frequency),which shows that Karama region using W-I model have larger path loss than Almajmoa'a region using OH model about 7 - 10dB, because of rooftop and multi-screen effects in karama region which makes more diffraction for signal.

**Figure 13.** Path loss versus distance in karama and Almajmoa'a regions

The path loss verses base stations height for O-H model and W-I model is shown in fig. 14., this result show that the path loss verses BS antenna height in karama district is higher than that for Almajmoa'a district, which can be attributed to high signal diffraction on roof top in karama region, and this figure shows that the decrease path loss in karama region greater than that in Almajmoa'a region by 7dB, due to the signals in karama region are release from rooftop and multi-screen effects for increasing BS height over the rooftop of karama regular buildings.

Mobile Radio Propagation Prediction for Two Different Districts in Mosul-City 35

In this study, two districts in Mosul city were investigated, one (Karama district) which has large number of similar and uniform building, Walfisch-Ikegami model is used which is suitable for radio network planning. The other is the Almajmoa'a district, which has dissimilar and irregular building, and less buildings density than Karama, Okumura-Hata model is applicable. Statistical Models are easier to be used than the physical ones. They do not need e.g. geographic databases. However, validity domain is often limited: Okumura-

Theodore S Rappaport, "Wireless communications principles and practice". John

J.D.Parson '' The mobile radio propagation channel '' second edition . John Wiely&Sons,

M.O. Kabaou."Multipath propagation models for radio mobile channel " Fourth international multi-conference on systems ,signals & devices .vol .3 2007 . Tunisia . S.A.Mawjoud "Estimation of design parameters for cellular WCDMA network" .AL-

M. J. Nawrocki M. Dohler "Understanding UMTS Radio Network Modelling, Planning and

F.I . Mahmoud & S.A.Mawjoud " planning and design of a WCDMA network compatible with existing GSM system in Mosul city " 5th international multi-conference on systems

M. Rahnema, "UMTS Network Planning Optimization, and Inter-Operation with GSM",

L. Poole, "Cellular Communications Explained from Basics to 3G", First Edition, Elsevier

E. Damasso and L. M. Correia, Eds., Digital Mobile Radio towards Future Generation

H. Holma and A. Toskala, WCDMA for UMTS: HSPA Evolution and LTE, 4th ed. John

S. Tabbane. *Handbook of Mobile Radio Networks*. Artech House Boston London, 2000.

V. Erceg. *Channels Modeling suitable for MBWA*. Activity report, Junary 2003.

Hata model cannot be used for distances less than 1 km.

*Electrical Department, College of Engineering, University of Mosul, Iraq* 

Rafidain engineering journal .vol 16. No .4 . Oct. 2008.

Automated Optimisation" John Wiley & Sons Ltd . 2006.

First Edition, John Wiley and Sons (Asia) Ltd, India, 2008.

J. S. Seybold, Introduction to RF Propagation. Wiley-IEEE, 2005.

,signals & devices. IEEE . 2008. Jordan .

Systems, COST 231 Final Report, 1999.

**8. Conclusions** 

**Author details** 

**9. References** 

2000.

Wiely&Sons, 2002.

Ltd, England, 2006.

Wiley & Sons, 2007.

Farhad E. Mahmood

**Figure 14.** Path loss versus BS antenna height.

## **8. Conclusions**

34 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

buildings.

**Figure 14.** Path loss versus BS antenna height.

**pathloss in dB**

The path loss verses base stations height for O-H model and W-I model is shown in fig. 14., this result show that the path loss verses BS antenna height in karama district is higher than that for Almajmoa'a district, which can be attributed to high signal diffraction on roof top in karama region, and this figure shows that the decrease path loss in karama region greater than that in Almajmoa'a region by 7dB, due to the signals in karama region are release from rooftop and multi-screen effects for increasing BS height over the rooftop of karama regular

15 20 25 30 35 40 45 50

**pathloss vs. BS antenna height for Almajmoaa and karama regions**

**karama district Almajmoaa district**

**BS antenna height (m)**

In this study, two districts in Mosul city were investigated, one (Karama district) which has large number of similar and uniform building, Walfisch-Ikegami model is used which is suitable for radio network planning. The other is the Almajmoa'a district, which has dissimilar and irregular building, and less buildings density than Karama, Okumura-Hata model is applicable. Statistical Models are easier to be used than the physical ones. They do not need e.g. geographic databases. However, validity domain is often limited: Okumura-Hata model cannot be used for distances less than 1 km.

## **Author details**

Farhad E. Mahmood *Electrical Department, College of Engineering, University of Mosul, Iraq* 

## **9. References**

	- J. Laiho, A. Wacker, and T. Novosad, Radio Network Planning and Optimisation for UMTS, 2nd ed. John Wiley & Sons, 2006.

**Chapter 0**

**Chapter 3**

**Two Novel Implementations of the Remez Multiple**

**Exchange Algorithm for Optimum FIR Filter Design**

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

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

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

> ©2012 Ahsan and Saramäki, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is

©2012 Ahsan and Saramäki, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

Muhammad Ahsan and Tapio Saramäki

known as the Parks-McClellan (PM) algorithm.

properly cited.

http://dx.doi.org/10.5772/46477

**1. Introduction**

doubt is superfluous.

Additional information is available at the end of the chapter

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

Muhammad Ahsan and Tapio Saramäki

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46477

## **1. Introduction**

36 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

2nd ed. John Wiley & Sons, 2006.

J. Laiho, A. Wacker, and T. Novosad, Radio Network Planning and Optimisation for UMTS,

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

©2012 Ahsan and Saramäki, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Ahsan and Saramäki, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 38 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>3</sup>

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.

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 [3] and provides

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

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

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

*Core Discrete Approximation Problem***:** Given *nz* <sup>−</sup> <sup>1</sup> 1, the number of unknowns *<sup>a</sup>*[*n*] for *<sup>n</sup>* <sup>=</sup> 0, 1, . . . , *nz* − 2, and the grid points *grid*(*k*) included in the vector **grid** of length *ngrid*, which contains values between 0 and 0.5 2, along with the vectors **des** and **wt** of the same length *ngrid*, 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

<sup>1</sup> In this contribution, *nz* <sup>−</sup> 1 is 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

<sup>2</sup> In the original PM algorithm, this range is the baseband for the so-called normalized frequencies from which the

<sup>1</sup>≤*k*≤*ngrid*|*E*(*k*)|, (1a)

*<sup>ε</sup>* <sup>=</sup> max

*Alternation Theorem* should be satisfied in order to guarantee the optimality of the solution.

corresponding angular frequencies are obtained by multiplying these frequencies by 2*π* [17].

significant improvements in designing the multiband FIR filters.

presented in Section (5).

**2. Problem statement**

the following:

following quantity:

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 [3] 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**

2 Will-be-set-by-IN-TECH

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:

extremum.

reporeted in [2].

additional last extremum.

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

*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

*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

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

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

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 *nz* <sup>−</sup> <sup>1</sup> 1, the number of unknowns *<sup>a</sup>*[*n*] for *<sup>n</sup>* <sup>=</sup> 0, 1, . . . , *nz* − 2, and the grid points *grid*(*k*) included in the vector **grid** of length *ngrid*, which contains values between 0 and 0.5 2, along with the vectors **des** and **wt** of the same length *ngrid*, 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:

$$\varepsilon = \max\_{1 \le k \le \eta \text{grid}} |E(k)|\,\tag{1a}$$

<sup>1</sup> In this contribution, *nz* <sup>−</sup> 1 is 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.

<sup>2</sup> 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].

#### 4 Will-be-set-by-IN-TECH 40 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>5</sup>

where

$$E(k) = \mathbf{wt}(k) \left[ G(k) - \mathbf{des}(k) \right] \tag{1b}$$

*3.1.1. Initialization phase*

*itrmax* = 250.

*3.1.2. Iteration phase*

for the integer part of *x*.

implementations of Segments 1, 2, and 3.

<sup>3</sup> It should be noted that the value of *dev* is either positive or negative.

to occur.

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 + (*m* − 1)�(*ngrid* − 1)/*nz*� for *m* = 1, 2, . . . , *nz* − 1 and **trial**(*nz*) = *ngrid*. Here, �*x*� stands

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

• Initialize the iteration counter as *niter* = 1 and set the maximum number of iterations as

The Remez exchange loop iteratively locates the desired optimum vector **trial** having as its entries the *nz* values of *k* within 1 ≤ *k* ≤ *ngrid*, 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 vector **opt**. As this is extremely unlikely to happen, **Segment 1**, based on the vector **opt**, generates the weighted error vector **wei\_err**(*k*) for 1 ≤ *k* ≤ *ngrid* corresponding to *E*(*k*), as given by (1b), in such a way that the *nz* − 1 unknowns *<sup>a</sup>*[0], *<sup>a</sup>*[1],..., *<sup>a</sup>*[*nz* <sup>−</sup> <sup>2</sup>] as well as *dev* <sup>3</sup> are implicitly found so that the following *nz* equations

are satisfied. When concentrating only on the values of *k* being the enteries of the "trial" vector **trial**, this solution is the best one according to *Alternation Theorem*. However, when

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

Finally, **Segment 4** checks whether **real** ≡ **trial** or not. If this equivalence is established, then the best solution has been found as in this case **opt** ≡ **real** ≡ **trial**. 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

**Segment 1:** After knowing the "trial" vector **trial** that contains the *nz* trial values of *k* in the ascending order for 1 ≤ *k* ≤ *ngrid* in the present iteration, this first segment guarantees that

considering all the values of *k* within 1 ≤ *k* ≤ *ngrid*, this solution is not the best one.

**wei\_err**(**trial**(*m*)) = (−1)*m*+1*dev* for *<sup>m</sup>* <sup>=</sup> 1, 2, . . . , *nz* (2)

and

$$G(k) = \sum\_{n=0}^{nz-2} a[n] \cos[2\pi n \cdot gri d(k)].\tag{1c}$$

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

$$\begin{aligned} \ell\_{\mathbf{opt}}(1) &< \ell\_{\mathbf{opt}}(2) < \dots < \ell\_{\mathbf{opt}}(nz - 1) < \ell\_{\mathbf{opt}}(nz) \\ \mathrm{E}[\ell\_{\mathbf{opt}}(m + 1)] &= -\mathrm{E}[\ell\_{\mathbf{opt}}(m)] \text{ for } m = 1, 2, \dots, nz - 1 \\ \left| \mathrm{E}[\ell\_{\mathbf{opt}}(m)] \right| &= \varepsilon \text{ for } m = 1, 2, \dots, nz . \end{aligned}$$

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*

4 Will-be-set-by-IN-TECH

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

**opt**(1) < **opt**(2) < ... < **opt**(*nz* − 1) < **opt**(*nz*)

= *ε* for *m* = 1, 2, . . . , *nz*.

*E*[**opt**(*m* + 1)] = −*E*[**opt**(*m*)] for *m* = 1, 2, . . . , *nz* − 1

It is worth mentioning that the core discrete approximation problem is the same for both the

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

*G*(*k*) =

*nz*−2 ∑ *n*=0

*E*(*k*) = **wt**(*k*)[*G*(*k*) − **des**(*k*)] (1b)

*a*[*n*] cos[2*πn* · *grid*(*k*)]. (1c)

where

and

1 ≤ *k* ≤ *ngrid* such that

**3. Implementation I**

the original problem are introduced.

**3.1. Theoretical formulation**

*E*[**opt**(*m*)]

implementations I and II as defined in the Introduction.

consideration with the aid of the proposed basic steps.

The overall implementation starts with the following initializations:


#### *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* within 1 ≤ *k* ≤ *ngrid*, 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 vector **opt**. As this is extremely unlikely to happen, **Segment 1**, based on the vector **opt**, generates the weighted error vector **wei\_err**(*k*) for 1 ≤ *k* ≤ *ngrid* corresponding to *E*(*k*), as given by (1b), in such a way that the *nz* − 1 unknowns *<sup>a</sup>*[0], *<sup>a</sup>*[1],..., *<sup>a</sup>*[*nz* <sup>−</sup> <sup>2</sup>] as well as *dev* <sup>3</sup> are implicitly found so that the following *nz* equations

$$\mathbf{wei\\_err}(\ell\_{\text{trial}}(m)) = (-1)^{m+1}dev \text{ for } m = 1, 2, \dots, nz \tag{2}$$

are satisfied. When concentrating only on the values of *k* being the enteries of the "trial" vector **trial**, this solution is the best one according to *Alternation Theorem*. However, when considering all the values of *k* within 1 ≤ *k* ≤ *ngrid*, 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 **real** ≡ **trial** or not. If this equivalence is established, then the best solution has been found as in this case **opt** ≡ **real** ≡ **trial**. 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 1 ≤ *k* ≤ *ngrid* in the present iteration, this first segment guarantees that

<sup>3</sup> It should be noted that the value of *dev* is either positive or negative.

#### 6 Will-be-set-by-IN-TECH 42 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>7</sup>

*Alternation Theorem* is satisfied when concentrating only on those *nz* values of *k* being involved in the vector **trial**. For this purpose, it generates the weighted error vector **wei\_err**(*k*) for 1 ≤ *k* ≤ *ngrid* such that the following system of *nz* equations:

$$\begin{aligned} \mathbf{wt}(\mathfrak{C}\_{\texttt{trial}}(m)) \left[ \sum\_{m=0}^{nz-2} a(n) \cos \left( 2\pi n \cdot \operatorname{grid}(\mathfrak{C}\_{\texttt{trial}}(m)) \right) - \mathbf{des}(\mathfrak{C}\_{\texttt{trial}}(m)) \right] \\ &= (-1)^{m+1} \operatorname{dev} \, \operatorname{for } m = 1, 2, \ldots, nz \end{aligned} \tag{3}$$

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

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

**aid**(*k*) = **ad**(*m*)

**y**(*m*)**aid**(*k*) and **err\_den**(*k*) = **err\_den**(*k*) + **aid**(*k*), respectively.

**wei\_err**(*k*)=[ **err\_num**(*k*)

**wei\_err**(**trial**(*m*)) =

set *m* = 1, and go to the next step.

the following three steps:

under **Segment 2**.

simplified as follows.

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

*Step 9:* Update the vector **wei\_err** as

**grid** as

**<sup>y</sup>**(*m*) = **des**[**trial**(*m*)] + (−1)*m*−1**ad**[**trial**(*m*)]/**wt**[**trial**(*m*)] for *<sup>m</sup>* <sup>=</sup> 1, 2, . . . , *nz*. (7)

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

**x\_all**(*k*) = cos[2*π* · **grid**(*k*)] for *k* = 1, 2, . . . , *ngrid*. (8)

+*dev* for *m* odd −*dev* for *m* even.

[**x\_all**(*k*) − **x**(*m*)] for *k* = 1, 2, . . . , *ngrid* (9)

**err\_den**(*k*) − **des**(*k*)]**wt**(*k*). (10)

(11)

*Step 3:* Generate the entries of the abscissa vector **x\_all** covering all the entries in the vector

*Step 4:* Select the entries of the vectors **err\_num** and **err\_den** of length *ngrid* to be zero valued,

and update **err\_num**(*k*) and **err\_den**(*k*) for *k* = 1, 2, . . . , *nz* as **err\_num**(*k*) = **err\_num**(*k*) +

*Step 6:* Set *m* = *m* + 1. If *m* > *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

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 is *dev* (−*dev*), where *dev* is given by (6). Hence, the vector **wei\_err** can be completed by using

*Step 10:* Set *m* = *m* + 1. If *m* < *nz* + 1, then go to the previous step. Otherwise, go to *Step 1*

**Segment 2:** This segment explains how to perform *Vicinity Search* based on the values of the weighted error function **wei\_err**(*k*) for 1 ≤ *k* ≤ *ngrid*, which has been generated at **Segment 1**, and the "trial" vector **trial**, which is under consideration in the present iteration. The key idea in *Vicinity Search* is to determine the *m*th entry of the "real" vector **real**, denoted by **real**(*m*) for *m* = 1, 2, . . . , *nz*, to be the value of *k* in the close vicinity of *k* = **real**(*m*), where a local extremum of **wei\_err**(*k*) with the same sign occurs. The location of these *nz* entries are

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 1 ≤ *k* ≤ *ngrid* by sign[**wei\_err**(**trial**(*m*))] as in this case the values of the resulting signed weighted function sign[**wei\_err**(**trial**(*m*))] ×

is implicitly solved for the *nz* <sup>−</sup> 1 unknowns *<sup>a</sup>*[0], *<sup>a</sup>*[1],..., *<sup>a</sup>*[*nz* <sup>−</sup> <sup>2</sup>] as well as for *dev*. <sup>4</sup> For this purpose, similar to the function **firpm**, for a given "trial" vector **trial**, the value of **wei\_err**(*k*) at *k* = **trial**(1), denoted by *dev*, the corresponding abscissa vector **x**, the ordinate vector **y**, and the coefficient vector **ad**, each of which are of length *nz*, 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* for 1 ≤ *k* ≤ *ngrid*, 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

$$\mathbf{x}(m) = \cos[2\pi \cdot \mathfrak{L}\_{\mathbf{trial}}(m)] \text{ for } m = 1, 2, \dots, nz \tag{4}$$

and

$$\mathbf{ad}(m) = 1 \Big/ \prod\_{\substack{k=1 \\ k \neq m}}^{nz} [\mathbf{x}(m) - \mathbf{x}(k)] \quad \text{for} \ m = 1, 2, \dots, nz,\tag{5}$$

respectively, as well as the corresponding deviation value as

$$dev = -\frac{\sum\_{m=1}^{nz} \mathbf{ad}(m) \mathbf{des}[\boldsymbol{\ell}\_{\text{trial}}(m)]}{\sum\_{m=1}^{nz} (-1)^{m-1} \mathbf{ad}(m) \mathbf{wt}[\boldsymbol{\ell}\_{\text{trial}}(m)]}.\tag{6}$$

<sup>4</sup> It is worth emphasizing that implicit solutions for the calculation of *a*[*n*]'s are not required for the intermediate iterations. The explicit solution for the calculation of *a*[*n*]'s is needed only after achieving the convergence to the best solution.

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

6 Will-be-set-by-IN-TECH

*Alternation Theorem* is satisfied when concentrating only on those *nz* values of *k* being involved in the vector **trial**. For this purpose, it generates the weighted error vector **wei\_err**(*k*) for

is implicitly solved for the *nz* <sup>−</sup> 1 unknowns *<sup>a</sup>*[0], *<sup>a</sup>*[1],..., *<sup>a</sup>*[*nz* <sup>−</sup> <sup>2</sup>] as well as for *dev*. <sup>4</sup> For this purpose, similar to the function **firpm**, for a given "trial" vector **trial**, the value of **wei\_err**(*k*) at *k* = **trial**(1), denoted by *dev*, the corresponding abscissa vector **x**, the ordinate vector **y**, and the coefficient vector **ad**, each of which are of length *nz*, 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* for 1 ≤ *k* ≤ *ngrid*, thereby making 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**

*a*(*n*) cos (2*πn* · *grid*(**trial**(*m*))) − **des**(**trial**(*m*))

**x**(*m*) = cos[2*π* · **trial**(*m*)] for *m* = 1, 2, . . . , *nz* (4)

*<sup>m</sup>*=<sup>1</sup> **ad**(*m*)**des**[**trial**(*m*)]

<sup>4</sup> It is worth emphasizing that implicit solutions for the calculation of *a*[*n*]'s are not required for the intermediate iterations. The explicit solution for the calculation of *a*[*n*]'s is needed only after achieving the convergence to the

[**x**(*m*) − **x**(*k*)] for *m* = 1, 2, . . . , *nz*, (5)

*<sup>m</sup>*=1(−1)*m*−1**ad**(*m*)**wt**[**trial**(*m*)] . (6)

= (−1)*m*+1*dev* for *<sup>m</sup>* <sup>=</sup> 1, 2, . . . , *nz* (3)

1 ≤ *k* ≤ *ngrid* such that the following system of *nz* equations:

implementation of the Remez loop very accurate and efficient.

This segment is performed by using the following ten steps:

 *nz* ∏ *k*=1 *k*�=*m*

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

**ad**(*m*) = 1

respectively, as well as the corresponding deviation value as

*dev* <sup>=</sup> <sup>−</sup> <sup>∑</sup>*nz*

∑*nz*

 *nz*−2 ∑ *m*=0

**wt**(**trial**(*m*))

**2** and *Endpoint Search* in **Segment 3**.

and

best solution.

$$\mathbf{y}(m) = \mathbf{des}[\mathfrak{E}\_{\text{trial}}(m)] + (-1)^{m-1} \mathbf{ad}[\mathfrak{E}\_{\text{trial}}(m)] / \mathbf{wt}[\mathfrak{E}\_{\text{trial}}(m)] \text{ for } m = 1, 2, \dots, n \text{z}. \tag{7}$$

*Step 3:* Generate the entries of the abscissa vector **x\_all** covering all the entries in the vector **grid** as

$$\mathbf{x}\_{-}\mathbf{all}(k) = \cos[2\pi \cdot \mathbf{grid}(k)] \text{ for } k = 1, 2, \dots, n \text{grid}. \tag{8}$$

*Step 4:* Select the entries of the vectors **err\_num** and **err\_den** of length *ngrid* to be zero valued, set *m* = 1, and go to the next step.

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

$$\mathbf{a}\mathbf{id}(k) = \mathbf{ad}(m) / [\mathbf{x}\_{-}\mathbf{all}(k) - \mathbf{x}(m)] \quad \text{for} \ k = 1, 2, \dots, n \text{grid} \tag{9}$$

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*) and **err\_den**(*k*) = **err\_den**(*k*) + **aid**(*k*), respectively.

*Step 6:* Set *m* = *m* + 1. If *m* > *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

$$\mathbf{\dot{\bf w}}\\_\mathbf{err}(k) = [\mathbf{\dot{\bf r}}\\_\mathbf{num}(k) / \mathbf{\dot{\bf r}}\\_\mathbf{den}(k) \ -\mathbf{des}(k)]\mathbf{wt}(k). \tag{10}$$

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 is *dev* (−*dev*), where *dev* is given by (6). 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

$$\text{fewi\\_err}(\ell\_{\text{trial}}(m)) = \begin{cases} +dev & \text{for } m \text{ odd} \\ -dev & \text{for } m \text{ even} \end{cases} \tag{11}$$

*Step 10:* Set *m* = *m* + 1. If *m* < *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*) for 1 ≤ *k* ≤ *ngrid*, which has been generated at **Segment 1**, and the "trial" vector **trial**, which is under consideration in the present iteration. The key idea in *Vicinity Search* is to determine the *m*th entry of the "real" vector **real**, denoted by **real**(*m*) for *m* = 1, 2, . . . , *nz*, to be the value of *k* in the close vicinity of *k* = **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 1 ≤ *k* ≤ *ngrid* 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 for *m* > 1 [*m* < *nz*], the "signs" of **wei\_err**(*k*) at both *k* = **trial**(*m* − 1) and *k* = **trial**(*m* + 1) is opposite to that of *k* = **trial**(*m*), or correspondingly, at *k* = **real**(*m*). Second, during the course of the present search, **real**(*m* − 1) is located before **real**(*m*) in such a way that the sign of **wei\_err**(*k*) at *k* = **real**(*m* − 1) is opposite to that of the sign of **wei\_err**(*k*) at *k* = **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 *m*th element, denoted by **real**(*m*), at which the vector

$$\mathbf{err\\_vicinity} = \text{sign}[\mathbf{wei\\_err}(\ell\_{\text{trial}}(m))] \mathbf{wei\\_err}(low\\_upper),\tag{12a}$$

is larger than or equal to 1. If this fact holds true, then the largest entry in the sub-vector −sign[**wei\_err**(**real**(1))] · **wei\_err**(1 : *upp***end**) should be positive. Similarly, the existence of

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

is smaller than or equal to *ngrid*. If this fact holds true, then the largest entry in the subvector

If no additional extremum exists, then the "final" **real** is the one found in *Vicinity Search*. Otherwise (that is, one or both the additional extrema exist), the final **real** is constructed

*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

*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

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* for 1 ≤ *k* ≤ *ngrid*, where the additional last (first) maximum of the signed weighted error function

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

*Step 2*: Determine *upp***end** *according to* (13). If *upp***end** = 0, then set **err\_end**(1) = 0. Otherwise,

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

*Step 3*: If **err\_end**(1) < **err\_vic**(*nz*), where **err\_vic**(*nz*) has been saved at *Step 2* under **Segment**

*Step 4*: Determine *low***end** according to (14). If *low***end** = *ngrid* + 1, then go to *Step 6*. Otherwise,

*Step 5*: If **err\_end**(*nz*) < max{**err\_end**(1), **err\_vic**(1)}, where **err\_vic**(1) has been saved at

�**end\_real**(*nz*), where the vector

*Step 2* under **Segment 2**, then go to the next step. Otherwise, set *endsearch* = 2.

**err\_endpoint** = −sign[**wei\_err**(**real**(1))] · **wei\_err**(1 : **end**(1)) (15)

**err\_endpoint** = −sign[**wei\_err**(**real**(*nz*))] · **wei\_err**(*low***end** : *ngrid*) (16)

**end\_real**(*nz*) + *low***end** − 1 and store

−sign[**wei\_err**(**real**(*nz*))] · **wei\_err**(*k*) [−sign[**wei\_err**(**real**(1))] · **wei\_err**(*k*)] occurs.

to exactly follow this strategy, it is carried out by using the following eight steps:

find the index, denoted by **end\_real**(1), where the vector

**2**, then go to the next step. Otherwise, set *endsearch* = 1.

achieves its maximum entry value. Set **end\_real**(*nz*) = ˜

the corresponding maximum entry value as **err\_end**(*nz*).

−sign[**wei\_err**(**real**(*nz*))] · **wei\_err**(*low***end** : *ngrid*) should be positive.

*low***end** = max{**trial**(*nz*) + 1, **real**(*nz*) + 1} (14)

the additional last local extremum implies that

according to the following alternatives:

extremum.

additional last extremum.

*Step 1*: Set *endsearch* = 0.

find the index, denoted by ˜

**err\_end**(1).

where

$$low = \begin{cases} 1 & \text{for } m = 1\\ \max\{\ell\_{\text{trial}}(m-1) + 1, \ell\_{\text{real}}(m-1) + 1\} & \text{for } 2 \le m \le n \end{cases} \tag{12b}$$

and

$$
\mu pp = \begin{cases}
\ell\_{\text{trial}}(m+1) - 1 & \text{for } 1 \le m \le nz - 1 \\
\eta grid & \text{for } m = nz
\end{cases}
\tag{12c}
$$

achieves the maximum value. Generate **real**(*m*) = **real**(*m*) + *low* − 1. If *m* = 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*) as **real**(*m*) = **real**(*m*) + *low* − 1.

*Step 3:* Set *m* = *m* + 1. If *m* < *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) [after *k* = **real**(*nz*)] such that its sign is opposite to that of occurring at *k* = **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 <sup>5</sup>

$$\mu p p\_{\mathbf{end}} = \min \{ \ell\_{\text{trial}}(1) - 1, \ell\_{\text{real}}(1) - 1 \} \tag{13}$$

<sup>5</sup> *Vicinity Search* automatically guarantees that the sign of the weighted error function **wei\_err**(*k*) is same at both *k* = **trial**(1) and *k* = **real**(1). Hence, *upp***end** is the smallest value of *k*, 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 *low***end** as given by (14) is the largest value of *k*, where the sign of **wei\_err**(*k*) can be opposite after these values.

is larger than or equal to 1. If this fact holds true, then the largest entry in the sub-vector −sign[**wei\_err**(**real**(1))] · **wei\_err**(1 : *upp***end**) should be positive. Similarly, the existence of the additional last local extremum implies that

$$low\_{\texttt{end}} = \max\{\ell\_{\texttt{trial}}(nz) + 1, \ell\_{\texttt{real}}(nz) + 1\}\tag{14}$$

is smaller than or equal to *ngrid*. If this fact holds true, then the largest entry in the subvector −sign[**wei\_err**(**real**(*nz*))] · **wei\_err**(*low***end** : *ngrid*) should be positive.

If no additional extremum exists, then the "final" **real** is 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* for 1 ≤ *k* ≤ *ngrid*, 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*: Set *endsearch* = 0.

8 Will-be-set-by-IN-TECH

**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 for *m* > 1 [*m* < *nz*], the "signs" of **wei\_err**(*k*) at both *k* = **trial**(*m* − 1) and *k* = **trial**(*m* + 1) is opposite to that of *k* = **trial**(*m*), or correspondingly, at *k* = **real**(*m*). Second, during the course of the present search, **real**(*m* − 1) is located before **real**(*m*) in such a way that the sign of **wei\_err**(*k*) at *k* = **real**(*m* − 1) is opposite to that of the sign of **wei\_err**(*k*) at *k* = **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

1 for *m* = 1 max{**trial**(*m* − 1) + 1, **real**(*m* − 1) + 1} for 2 ≤ *m* ≤ *nz*

*ngrid* for *m* = *nz*

store the corresponding maximum value as **err\_vic**(1) [**err\_vic**(*nz*)] to be used in *Endpoint*

*Step 3:* Set *m* = *m* + 1. If *m* < *nz* + 1, then go to the previous step. Otherwise, go to *Step 1*

**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) [after *k* = **real**(*nz*)] such that its sign is opposite to that of occurring at *k* = **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.

<sup>5</sup> *Vicinity Search* automatically guarantees that the sign of the weighted error function **wei\_err**(*k*) is same at both *k* = **trial**(1) and *k* = **real**(1). Hence, *upp***end** is the smallest value of *k*, 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 *low***end** as given by

**trial**(*m* + 1) − 1 for 1 ≤ *m* ≤ *nz* − 1

**real**(*m*) + *low* − 1.

*upp***end** = min{**trial**(1) − 1, **real**(1) − 1} (13)

**real**(*m*), at which the vector

**err\_vicinity** = sign[**wei\_err**(**trial**(*m*))]**wei\_err**(*low* : *upp*), (12a)

(12b)

(12c)

**real**(*m*) + *low* − 1. If *m* = 1 [*m* = *nz*], then

the following three steps:

where

and

under **Segment 3**.

*Step 1:* Set *m* = 1 and go to the next step. *Step 2:* Find the *m*th element, denoted by

*upp* =

achieves the maximum value. Generate **real**(*m*) =

*Search* at **Segment 3** . Update **real**(*m*) as **real**(*m*) =

The appearance of the additional first local extremum implies that <sup>5</sup>

(14) is the largest value of *k*, where the sign of **wei\_err**(*k*) can be opposite after these values.

*low* =

*Step 2*: Determine *upp***end** *according to* (13). If *upp***end** = 0, then set **err\_end**(1) = 0. Otherwise, find the index, denoted by **end\_real**(1), where the vector

$$\mathbf{err\\_endpoint} = -\text{sign}[\mathbf{wei\\_err}(\ell\_{\text{real}}(1))] \cdot \mathbf{wei\\_err}(1:\ell\_{\text{end}}(1)) \tag{15}$$

achieves the corresponding maximum entry value. Store this maximum entry value as **err\_end**(1).

*Step 3*: If **err\_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, set *endsearch* = 1.

*Step 4*: Determine *low***end** according to (14). If *low***end** = *ngrid* + 1, then go to *Step 6*. Otherwise, find the index, denoted by ˜ �**end\_real**(*nz*), where the vector

$$\mathbf{err\\_endpoint} = -\text{sign}[\mathbf{wei\\_err}(\ell\_{\text{real}}(nz))] \cdot \mathbf{wei\\_err}(low\_{\text{end}} : \eta grid) \tag{16}$$

achieves its maximum entry value. Set **end\_real**(*nz*) = ˜ **end\_real**(*nz*) + *low***end** − 1 and store the corresponding maximum entry value as **err\_end**(*nz*).

*Step 5*: If **err\_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, set *endsearch* = 2.

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

*Step 7*: If *endsearch = 1*, then set **real**(*nz* + 1 − *m*) = **real**(*nz* − *m*) for *m* = 1, 2, . . . , *nz* − 1 and **real**(1) = **end\_real**(1) and got to *Step 1* under **Segment 4**. Otherwise, go to the next step.

% WEIGHTED ERROR FUNCTION wei\_err(k) AT ALL THE GRID POINTS

dev = -dnum/dden; % Step 1: Current value of deviation.

A = x'\*ones(1,nz)-ones(nz,1)\*x;

 % Step 2: Lagrange ordinate vector y y = des(l\_trial) + dev\*add./wt(l\_trial); % Step 3: Overall abscissa vector x\_all x\_all = cos(2\*pi\*grid(1:ngrid));

err\_den = err\_num; % and err\_den.

err\_num = err\_num + y(jj)\*aid;

% SEGMENT 2: PERFORM VICINITY SEARCH

wei\_err(low:l\_trial(2)-1);

wei\_err(low:l\_trial(k+1)-1);

l\_real(k) = ind\_vicinity+low-1;

% SEGMENT 3: PERFORM ENDPOINT SEARCH

err\_vic(k) = wei\_err (l\_real(k));

[~,ind\_endpoint]=max(err\_endpoint);

endsearch=0; % Step 1: Start Endpoint search.

err\_cy = err\_num./err\_den;

 A(eye(nz)==1) = 1; ad = prod(A);

 add = ones(size(ad)); add(2:2:nz) = -add(2:2:nz); dnum = ad\*des(l\_trial)'; dden = add\*(ad./wt(l\_trial))';

end

 if k==1 low = 1;

else

end

elseif k==nz

x = cos(2\*pi\*grid(l\_trial)); % Step 1: Lagrange abscissa vector x.

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

 ad = ad \* (-2)^(nz-1); % Step 1: Lagrange coefficient vector ad... ad = 1./ad; % found efficiently without using the function remezdd.

err\_num = zeros(1,ngrid); % Step 4: Initializations of err\_num...

 wei\_err = (err\_cy - des).\*wt; % Step 7: Generate the vector wei\_err. dev\_vect = ones(size(l\_trial)); % Steps 8-10: Fill in the undefined dev\_vect(2:2:length(l\_trial))= -dev\_vect(2:2:length(l\_trial)); dev\_vect = dev\_vect \* dev; % entries of wei\_err at l\_trial(1:nz)... wei\_err(l\_trial)=dev\_vect; % by using the values of dev (-dev).

for k=1:nz % Steps 1,2, and 3: Start of Vicinity search.

err\_vicinity = sign(wei\_err(l\_trial(k)))\*wei\_err(low:ngrid);

 [~,ind\_vicinity]=max(err\_vicinity); % tilde operator does not... % work with older MATLAB releases If you are running an older... % MATLAB version, considering replacing it with a dummy value.

if k==1 || k==nz % Step 3: Find err\_vic(1)=wei\_err(l\_real(1))...

 end % and err\_vic(nz)=wei\_err(l\_real(nz)) for use at STEP III. end % Steps 1, 2, and 3: End of Vicinity search.

 err\_end(1) = 0; % Step 2: Needed for the case, where upp\_end = 0. if l\_real(1)>1 && l\_trial(1)> 1 % Step 2: Find l\_end\_true(1)... upp\_end = min(l\_real(1)-1,l\_trial(1)-1); % and err\_end(1). err\_endpoint = -sign(wei\_err(l\_real(1)))\*wei\_err(1: upp\_end);

err\_vicinity = sign(wei\_err(l\_trial(k)))\*...

low = max(l\_trial(k-1)+1,l\_real(k-1)+1);

 low = max(l\_trial(k-1)+1,l\_real(k-1)+1); err\_vicinity = sign(wei\_err(l\_trial(k)))\* ...

 for jj = 1:nz % Steps 5 and 6: Intermediate evaluations for... aid = ad(jj)./(x\_all - x(jj)); % obtaining the weighted error... err\_den = err\_den + aid; % wei\_err(k) at all the grid points.

*Step 8*: If *endsearch* = 2, set **real**(*m*) = **real**(*m* + 1) for *m* = 1, 2, . . . , *nz* − 1 and **real**(*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 **opt** ≡ **real** ≡ **trial** 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* exceeds *itrmax* = 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** and *niter* = *niter* + 1. If *niter* > *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 version 7.11.0.584(*R*2010*b*). In order to embed this code snippet in the MATLAB function **firpm**, edit the function by taking away the code between lines 214 and 514, 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.

```
1 function [x,y,ad,dev] = remez_imp1(nz,iext,ngrid,grid,des,wt)
2 % remez_imp1 implements the Segments 1 - 4 described in the preceding
3 % section, the function needs to be inserted within the MATLAB function
4 % firpm. The input argument values come directly from the function firpm
5 % and the output arguments are required to perform the Inverse Fourier
6 % transform in order to calculate the filter coefficients. In case of
7 % any issues send an e-mail to muhammad"dot"ahsan"at"tut "dot" fi.
8 % Last updated 04.15.2012 4:15 AM (UTC/GMT+2)
9
10 % INITIALIZATIONS PHASE
11 niter = 1; % Initialize the iteration counter.
12 itrmax = 250; % Maximum number of iterations.
13 l_trial = iext(1:nz)'; % Startup value of l_trial.
14
15 % ITERATION PHASE
16 % REMEZ LOOP FOR LOCATING DESIRED nz INDICES AMONG THE GRID POINTS
17 while (niter < itrmax)
18
19 % SEGMENT 1: BASED ON THE PRESENT 'TRIAL' VECTOR l_trial, GENERATE THE
```
Will-be-set-by-IN-TECH

*Step 7*: If *endsearch = 1*, then set **real**(*nz* + 1 − *m*) = **real**(*nz* − *m*) for *m* = 1, 2, . . . , *nz* − 1 and **real**(1) = **end\_real**(1) and got to *Step 1* under **Segment 4**. Otherwise, go to the next step.

*Step 8*: If *endsearch* = 2, set **real**(*m*) = **real**(*m* + 1) for *m* = 1, 2, . . . , *nz* − 1 and **real**(*nz*) =

**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 **opt** ≡ **real** ≡ **trial** 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* exceeds *itrmax* = 250, which is extremely unlikely to

*Step 2*: Set **trial** = **real** and *niter* = *niter* + 1. If *niter* > *itrmax*, then stop. Otherwise, go to

The code pasted below has been tested for realizing Implementation I by using MATLAB version 7.11.0.584(*R*2010*b*). In order to embed this code snippet in the MATLAB function **firpm**, edit the function by taking away the code between lines 214 and 514, 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.

**end\_real**(*nz*). Go to *Step 1* under **Segment 4**.

This segment requires only the following two basic steps:

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

function [x,y,ad,dev] = remez\_imp1(nz,iext,ngrid,grid,des,wt)

% Last updated 04.15.2012 4:15 AM (UTC/GMT+2)

 niter = 1; % Initialize the iteration counter. itrmax = 250; % Maximum number of iterations. l\_trial = iext(1:nz)'; % Startup value of l\_trial.

% REMEZ LOOP FOR LOCATING DESIRED nz INDICES AMONG THE GRID POINTS

% SEGMENT 1: BASED ON THE PRESENT 'TRIAL' VECTOR l\_trial, GENERATE THE

 % remez\_imp1 implements the Segments 1 - 4 described in the preceding % section, the function needs to be inserted within the MATLAB function % firpm. The input argument values come directly from the function firpm % and the output arguments are required to perform the Inverse Fourier % transform in order to calculate the filter coefficients. In case of % any issues send an e-mail to muhammad"dot"ahsan"at"tut "dot" fi.

occur.

*Step 1* under **Segment 1**.

**3.2. MATLAB code snippet**

% INITIALIZATIONS PHASE

% ITERATION PHASE

while (niter < itrmax)

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

```
20 % WEIGHTED ERROR FUNCTION wei_err(k) AT ALL THE GRID POINTS
21 x = cos(2*pi*grid(l_trial)); % Step 1: Lagrange abscissa vector x.
22 A = x'*ones(1,nz)-ones(nz,1)*x;
23 A(eye(nz)==1) = 1;
24 ad = prod(A);
25 ad = ad * (-2)^(nz-1); % Step 1: Lagrange coefficient vector ad...
26 ad = 1./ad; % found efficiently without using the function remezdd.
27 add = ones(size(ad));
28 add(2:2:nz) = -add(2:2:nz);
29 dnum = ad*des(l_trial)';
30 dden = add*(ad./wt(l_trial))';
31 dev = -dnum/dden; % Step 1: Current value of deviation.
32 % Step 2: Lagrange ordinate vector y
33 y = des(l_trial) + dev*add./wt(l_trial);
34 % Step 3: Overall abscissa vector x_all
35 x_all = cos(2*pi*grid(1:ngrid));
36 err_num = zeros(1,ngrid); % Step 4: Initializations of err_num...
37 err_den = err_num; % and err_den.
38 for jj = 1:nz % Steps 5 and 6: Intermediate evaluations for...
39 aid = ad(jj)./(x_all - x(jj)); % obtaining the weighted error...
40 err_den = err_den + aid; % wei_err(k) at all the grid points.
41 err_num = err_num + y(jj)*aid;
42 end
43 err_cy = err_num./err_den;
44 wei_err = (err_cy - des).*wt; % Step 7: Generate the vector wei_err.
45 dev_vect = ones(size(l_trial)); % Steps 8-10: Fill in the undefined
46 dev_vect(2:2:length(l_trial))= -dev_vect(2:2:length(l_trial));
47 dev_vect = dev_vect * dev; % entries of wei_err at l_trial(1:nz)...
48 wei_err(l_trial)=dev_vect; % by using the values of dev (-dev).
49
50 % SEGMENT 2: PERFORM VICINITY SEARCH
51 for k=1:nz % Steps 1,2, and 3: Start of Vicinity search.
52 if k==1
53 low = 1;
54 err_vicinity = sign(wei_err(l_trial(k)))*...
55 wei_err(low:l_trial(2)-1);
56 elseif k==nz
57 low = max(l_trial(k-1)+1,l_real(k-1)+1);
58 err_vicinity = sign(wei_err(l_trial(k)))*wei_err(low:ngrid);
59 else
60 low = max(l_trial(k-1)+1,l_real(k-1)+1);
61 err_vicinity = sign(wei_err(l_trial(k)))* ...
62 wei_err(low:l_trial(k+1)-1);
63 end
64 [~,ind_vicinity]=max(err_vicinity); % tilde operator does not...
65 % work with older MATLAB releases If you are running an older...
66 % MATLAB version, considering replacing it with a dummy value.
67 l_real(k) = ind_vicinity+low-1;
68 if k==1 || k==nz % Step 3: Find err_vic(1)=wei_err(l_real(1))...
69 err_vic(k) = wei_err (l_real(k));
70 end % and err_vic(nz)=wei_err(l_real(nz)) for use at STEP III.
71 end % Steps 1, 2, and 3: End of Vicinity search.
72
73 % SEGMENT 3: PERFORM ENDPOINT SEARCH
74 endsearch=0; % Step 1: Start Endpoint search.
75 err_end(1) = 0; % Step 2: Needed for the case, where upp_end = 0.
76 if l_real(1)>1 && l_trial(1)> 1 % Step 2: Find l_end_true(1)...
77 upp_end = min(l_real(1)-1,l_trial(1)-1); % and err_end(1).
78 err_endpoint = -sign(wei_err(l_real(1)))*wei_err(1: upp_end);
79 [~,ind_endpoint]=max(err_endpoint);
```
The minimum order to meet these criteria is 108 and the relevant MATLAB commands are

Angular frequency

*ω<sup>s</sup>* = 0.02*π*, *ω<sup>p</sup>* = 0.05*π*, *δ<sup>p</sup>* = 0.01, and *δ<sup>s</sup>* = 0.001.

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

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

1 >> [N,F,A,W] = firpmord([0.2 0.25 0.6 0.7],[0 1 0],[0.001 .01 .01]);

*ωs*<sup>1</sup> = 0.2*π*, *ωp*<sup>1</sup> = 0.25*π*, *ωp*<sup>2</sup> = 0.6*π*, *ωs*<sup>2</sup> = 0.7*π*, *δ<sup>p</sup>* = *δs*<sup>2</sup> = 0.01, and *δs*<sup>1</sup> = 0.001.

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

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

Passband details (linear scale)

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

0 0.01π 0.02π 0.03π 0.04π 0.05π

1 >> [N,F,A,W] = firpmord([0.05 0.1],[1 0],[0.01 0.001]);

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

0.99 1 1.01

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

3 >> fvtool(firr\_coeff); % filter visualization tool

2 >> firr\_coeff = firremez\_imp1(N+6,F,A,W);

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

**Figure 1.** Magnitude response of the lowpass filter of Example 1.

1 >> [N,F,A,W] = firpmord([0.02 0.05],[0 1],[0.001 0.01]);

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

2 >> firr\_coeff = firremez\_imp1(N-4,F,A,W);

2 >> firr\_coeff = firremez\_imp1(N,F,A,W);

are

Magnitude (dB)

```
80 l_end_real(1) = ind_endpoint;
81 err_end(1) = -sign(wei_err(l_real(1)))*wei_err(l_end_real(1));
82 if err_end(1) > abs(err_vic(nz)) % Step 3:Use 'endsearch=1'...
83 endsearch=1; % or not?
84 end
85 end
86 if l_real(nz) < ngrid & l_trial(nz) < ngrid % Step 4: Find...
87 low_end = max(l_real(nz)+1,l_trial(nz)+1); % l_end_real(nz)...
88 err_endpoint = -sign(wei_err(l_real(nz)))*wei_err(low_end:ngrid);
89 [~,ind_endpoint]=max(err_endpoint); % and err_end(nz).
90 l_end_real(nz) = ind_endpoint+low_end-1;
91 err_end(nz) = -sign(wei_err(l_real(nz)))*wei_err(l_end_real(nz));
92 if err_end(nz) > max(abs(err_vic(1)), err_end(1)) % Step 5:...
93 endsearch=2; % Use 'endsearch=2' or not?
94 end
95 end
96 if endsearch == 1 % Step 7: 'endsearch=1' is valid. Form...
97 l_real=[l_end_real(1) l_real(1:nz-1)]; % l_real accordingly.
98 elseif endsearch == 2 % Step 8: 'endsearch=2' is true. Form...
99 l_real=[l_real(2:nz) l_end_real(nz)]; % l_real accordingly.
100 end % End of Endpoint search.
101
102 % SEGMENT 4: TEST CONVERGENCE
103 if (l_real == l_trial) % Step 1: The real and trial vectors...
104 break; % coincide. Hence, stop. Remez loop ended successfully.
105 else
106 l_trial = l_real; % Step 2: Otherwise, replace the values of...
107 niter = niter + 1; % l_trial with the values of l_real and...
108 end % continue.
109 end % END OF THE OVERALL REMEZ LOOP
```
#### **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:


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:

*ω<sup>p</sup>* = 0.05*π*, *ω<sup>s</sup>* = 0.1*π*, *δ<sup>p</sup>* = 0.01, and *δ<sup>s</sup>* = 0.001.

48 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>13</sup> Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design 49

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

```
1 >> [N,F,A,W] = firpmord([0.05 0.1],[1 0],[0.01 0.001]);
2 >> firr_coeff = firremez_imp1(N+6,F,A,W);
3 >> fvtool(firr_coeff); % filter visualization tool
```
12 Will-be-set-by-IN-TECH

<sup>81</sup> err\_end(1) = -sign(wei\_err(l\_real(1)))\*wei\_err(l\_end\_real(1)); 82 if err\_end(1) > abs(err\_vic(nz)) % Step 3:Use 'endsearch=1'...

86 if l\_real(nz) < ngrid & l\_trial(nz) < ngrid % Step 4: Find... 87 low\_end = max(l\_real(nz)+1,l\_trial(nz)+1); % l\_end\_real(nz)... <sup>88</sup> err\_endpoint = -sign(wei\_err(l\_real(nz)))\*wei\_err(low\_end:ngrid);

<sup>91</sup> err\_end(nz) = -sign(wei\_err(l\_real(nz)))\*wei\_err(l\_end\_real(nz)); 92 if err\_end(nz) > max(abs(err\_vic(1)), err\_end(1)) % Step 5:...

89 [~,ind\_endpoint]=max(err\_endpoint); % and err\_end(nz).

 if endsearch == 1 % Step 7: 'endsearch=1' is valid. Form... l\_real=[l\_end\_real(1) l\_real(1:nz-1)]; % l\_real accordingly. elseif endsearch == 2 % Step 8: 'endsearch=2' is true. Form... l\_real=[l\_real(2:nz) l\_end\_real(nz)]; % l\_real accordingly. end % End of Endpoint search.

103 if (l\_real == l\_trial) % Step 1: The real and trial vectors... 104 break; % coincide. Hence, stop. Remez loop ended successfully.

106 l\_trial = l\_real; % Step 2: Otherwise, replace the values of... 107 niter = niter + 1; % l\_trial with the values of l\_real and...

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

• Total Time − The total time spent in a function, including all child functions called, in

• Self Time − The total time taken by an individual function, not including the time for any

Time measurement was carried out on an IBM ThinkCentre machine equipped with Intel Core

The following four filters were designed with both implementations. After the examples, the

*ω<sup>p</sup>* = 0.05*π*, *ω<sup>s</sup>* = 0.1*π*, *δ<sup>p</sup>* = 0.01, and *δ<sup>s</sup>* = 0.001.

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.

2 Duo processor E6550 running at a speed of 2.33 GHz with a memory of 3 GB.

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

90 l\_end\_real(nz) = ind\_endpoint+low\_end-1;

93 endsearch=2; % Use 'endsearch=2' or not?

80 l\_end\_real(1) = ind\_endpoint;

83 endsearch=1; % or not?

84 end 85 end

94 end 95 end

105 else

seconds.

108 end % continue.

102 % SEGMENT 4: TEST CONVERGENCE

**3.3. Performance comparison**

109 end % END OF THE OVERALL REMEZ LOOP

child functions called, in seconds.

time taken by them will be tabulated in Table 1.

101

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

**Figure 1.** Magnitude response of the lowpass filter of Example 1.

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

*ω<sup>s</sup>* = 0.02*π*, *ω<sup>p</sup>* = 0.05*π*, *δ<sup>p</sup>* = 0.01, and *δ<sup>s</sup>* = 0.001.

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

1 >> [N,F,A,W] = firpmord([0.02 0.05],[0 1],[0.001 0.01]); 2 >> firr\_coeff = firremez\_imp1(N-4,F,A,W);

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:

*ωs*<sup>1</sup> = 0.2*π*, *ωp*<sup>1</sup> = 0.25*π*, *ωp*<sup>2</sup> = 0.6*π*, *ωs*<sup>2</sup> = 0.7*π*, *δ<sup>p</sup>* = *δs*<sup>2</sup> = 0.01, and *δs*<sup>1</sup> = 0.001.

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

1 >> [N,F,A,W] = firpmord([0.2 0.25 0.6 0.7],[0 1 0],[0.001 .01 .01]); 2 >> firr\_coeff = firremez\_imp1(N,F,A,W);

• 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.

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

In order to conveniently avoid the appearance of the unwanted transition peaks, consider an

Such that these regions do not overlap. The lower and upper limits for the zero-phase

Here, the *Dκ*'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

<sup>1</sup> <sup>∪</sup> <sup>Ω</sup><sup>2</sup> <sup>∪</sup> <sup>Ω</sup>*<sup>T</sup>*

In order to guarantee that Ω *<sup>κ</sup>* is still a closed subset of [0, *π*], *α* should be a small positive number. <sup>6</sup> 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

min{*D<sup>κ</sup>* − *δκ*, *Dκ*+<sup>1</sup> − *δκ*+1}, for *Type A*

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

<sup>6</sup> The only condition for *α* is that it should be small enough to avoid the extra peaks between the adjacent passbands

<sup>7</sup> It is worth emphasizing that the use of max{*D<sup>κ</sup>* <sup>+</sup> *δκ*, *<sup>D</sup>κ*+<sup>1</sup> <sup>+</sup> *δκ*<sup>+</sup>1} implies that the maximum allowable value in the nearest passband is the upper limit. Similarly, in the following equation, min{*D<sup>κ</sup>* − *δκ*, *Dκ*+<sup>1</sup> − *δκ*<sup>+</sup>1} means that the

<sup>2</sup> <sup>∪</sup> ... <sup>Ω</sup>*<sup>T</sup>*

*<sup>κ</sup>* = [*ω*(*upp*) *<sup>κ</sup>* <sup>+</sup> *<sup>α</sup>*, *<sup>ω</sup>*(*low*) *<sup>κ</sup>* <sup>−</sup> *<sup>α</sup>*] for *<sup>κ</sup>* <sup>=</sup> 1, 2, . . . , *<sup>R</sup>* <sup>−</sup> 1. (17e)

*<sup>L</sup>*(*upp*) *<sup>κ</sup>* <sup>=</sup> max{*D<sup>κ</sup>* <sup>+</sup> *δκ*, *<sup>D</sup>κ*+<sup>1</sup> <sup>+</sup> *δκ*+1}, (17f)

*<sup>L</sup>*(*upp*) *<sup>κ</sup>* for *Type B* (17g)

for *κ* = 1, 2, . . . , *R*. (17a)

*<sup>R</sup>*−<sup>1</sup> <sup>∪</sup> <sup>Ω</sup>*R*, (17d)

*<sup>κ</sup>* are specified as follows. For

*<sup>L</sup>*(*low*) *<sup>κ</sup>* <sup>=</sup> *<sup>D</sup><sup>κ</sup>* <sup>−</sup> *δκ* and *<sup>L</sup>*(*upp*) *<sup>κ</sup>* <sup>=</sup> *<sup>D</sup><sup>κ</sup>* <sup>+</sup> *δκ* (17b)

Ω = Ω<sup>1</sup> ∪ Ω<sup>2</sup> ∪ ... ∪ Ω*R*. (17c)

original problem stated for linear-phase Type I and Type II filters as follows.

*<sup>ω</sup>*(*low*) *<sup>κ</sup>* , *<sup>ω</sup>*(*upp*) *<sup>κ</sup>*

First, there are *R* interlaced passband and stopband regions as given by

Ω*<sup>κ</sup>* = 

approximation region is

where

frequency response in these bands are, respectively, specified as

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

<sup>Ω</sup> <sup>=</sup> <sup>Ω</sup><sup>1</sup> <sup>∪</sup> <sup>Ω</sup>*<sup>T</sup>*

the zero-phase frequency response in the *κ*th transition band Ω*<sup>T</sup>*

whereas the lower limits depend on the type as follows:

and the newly formed intervals in the transition band regions.

−

Ω*<sup>T</sup>*

both *Type A* and *Type B*, the upper limit is <sup>7</sup>

 *<sup>L</sup>*(*low*) *<sup>κ</sup>* <sup>=</sup>

lower limit is the one in the nearest stopband.

**Figure 2.** Magnitude response of the highpass filter of Example 2.

**Figure 3.** Magnitude response and the highest transition band peak of the bandpass filter of Example 3.

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]:

50 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>15</sup> Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design 51


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

$$\Omega\_{\mathbb{K}} = \left[ \omega\_{\mathbb{K}}^{(low)}, \omega\_{\mathbb{K}}^{(up)} \right] \text{ for } \mathbb{x} = 1, 2, \dots, R. \tag{17a}$$

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\_{\mathbf{x}}^{(low)} = D\_{\mathbf{x}} - \delta\_{\mathbf{x}} \text{ and } L\_{\mathbf{x}}^{(upp)} = D\_{\mathbf{x}} + \delta\_{\mathbf{x}} \tag{17b}$$

Here, the *Dκ*'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

$$
\Omega = \Omega\_1 \cup \Omega\_2 \cup \dots \cup \Omega\_R. \tag{17c}
$$

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

$$
\widehat{\Omega} = \Omega\_1 \cup \Omega\_1^T \cup \Omega\_2 \cup \Omega\_2^T \cup \dots \cup \Omega\_{R-1}^T \cup \Omega\_{R'} \tag{17d}
$$

where

14 Will-be-set-by-IN-TECH

Passband details (linear scale)

0.2π 0.4π 0.6π 0.8π π

Angular frequency

Angular frequency

**Figure 3.** Magnitude response and the highest transition band peak of the bandpass filter of Example 3.

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

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

0.6π 0.65π 0.7π

Transition band peak (linear scale)

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

−100

of the following holds true [9]:

−80

−60

−40

Magnitude (dB)

−20

0

20

0.99 0.995 1 1.005 1.01

**Figure 2.** Magnitude response of the highpass filter of Example 2.

Magnitude (dB)

$$
\Omega\_{\mathbf{x}}^{T} = \left[ \omega\_{\mathbf{x}}^{(upp)} + \mathfrak{a}, \omega\_{\mathbf{x}}^{(low)} - \mathfrak{a} \right] \text{ for } \mathbf{x} = \mathbf{1}, \mathbf{2}, \dots, \mathbf{R} - \mathbf{1}. \tag{17e}
$$

In order to guarantee that Ω *<sup>κ</sup>* is still a closed subset of [0, *π*], *α* should be a small positive number. <sup>6</sup> 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 Ω*<sup>T</sup> <sup>κ</sup>* are specified as follows. For both *Type A* and *Type B*, the upper limit is <sup>7</sup>

$$\hat{L}\_{\mathbf{x}}^{(upper)} = \max \{ D\_{\mathbf{x}} + \delta\_{\mathbf{x}}, D\_{\mathbf{x}+1} + \delta\_{\mathbf{x}+1} \} \, \tag{17f}$$

whereas the lower limits depend on the type as follows:

$$
\hat{L}\_{\mathbf{x}}^{(low)} = \begin{cases}
\min\{D\_{\mathbf{x}} - \delta\_{\mathbf{x}\prime} D\_{\mathbf{x}+1} - \delta\_{\mathbf{x}+1}\}, & \text{for Type } A \\
\end{cases}
\tag{17g}
$$

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

<sup>6</sup> 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.

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

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 Ω*<sup>T</sup> <sup>κ</sup>* 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 Ω*<sup>T</sup> <sup>κ</sup>* for *κ* = 1, 2, . . . , *R* − 1 the desired values, denoted by *D<sup>κ</sup>*, and the admissible deviations, denoted by *δ <sup>κ</sup>*, are as follows:

$$\hat{D}\_{\mathbf{k}} = \begin{cases} \left[ \max\{ D\_{\mathbf{k}} + \delta\_{\mathbf{k}\prime} D\_{\mathbf{k}+1} + \delta\_{\mathbf{k}+1} \} + \min\{ D\_{\mathbf{k}} - \delta\_{\mathbf{k}\prime} D\_{\mathbf{k}+1} - \delta\_{\mathbf{k}+1} \} \right] / 2 & \text{for Type } \mathbf{A} \\ 0 & \text{for Type } \mathbf{B} . \end{cases} \tag{17b}$$

and

$$
\widehat{\delta}\_{\mathbb{K}} = \begin{cases}
\max\{D\_{\mathbb{K}} + \delta\_{\mathbb{K}}, D\_{\mathbb{K}+1} + \delta\_{\mathbb{K}+1}\} - \widehat{D}\_{\mathbb{K}} & \text{for Type } A \\
\max\{D\_{\mathbb{K}} + \delta\_{\mathbb{K}}, D\_{\mathbb{K}+1} + \delta\_{\mathbb{K}+1}\} & \text{for Type } B.
\end{cases}
\tag{17i}
$$

23 % Generate the output parameters

<sup>33</sup> Des(2\*k) = (aid1+aid2)/2; <sup>34</sup> Dev(2\*k) = aid1-Des(2\*k);

39 error('Type should be either 1 or 2.');

<sup>25</sup> F(4\*(k-1)+1)=F\_ori(2\*k-1); F(4\*(k-1)+2)=F\_ori(2\*k); <sup>26</sup> Des(2\*(k-1)+1)=Des\_ori(k); Dev(2\*(k-1)+1)=Dev\_ori(k);

<sup>29</sup> F(4\*(k-1)+3)=F\_ori(2\*k)+alpha;F(4\*(k-1)+4)=F\_ori(2\*k+1)-alpha; 30 aid1=max(Des\_ori(k)+Dev\_ori(k),Des\_ori(k+1)+Dev\_ori(k+1)); 31 aid2=min(Des\_ori(k)-Dev\_ori(k),Des\_ori(k+1)-Dev\_ori(k+1));

When using the above MATLAB function with *α* = 0.0005, the ten band-edges of the five

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

Ω *<sup>κ</sup>* = [0, 0.2, 0.2005, 0.2495, 0.25, 0.6, 0.6005, 0.6995, 0.7, 1] . The corresponding desired and weight values for *Type A* and *Type B* transition band

*WA* = [1000, 1.9782, 100, 1.9608, 100],

*WB* = [1000, 0.9901, 100, 0.9901, 100].

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.

bands as fractions of *π* for both *Type A* and *Type B* transition band constraints become

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

*DB* = [0, 0, 1, 0, 0]

2 >> [F1,A1,W1] = convert2constrt([0 0.2 0.25 0.6 0.7 1],[0 1 0],...

6 >> [F2,A2,W2] = convert2constrt([0 0.2 0.25 0.6 0.7 1],[0 1 0],... 7 [.001 .01 .01],0.0005,2); firr\_coeff2 = firremez\_imp1(N1,F2,A2,W2);

4 >> N1 = 103; firr\_coeff1 = firremez\_imp1(N1,F1,A1,W1);

24 for k=1:R

32 if itype==1

38 else

40 end 41 end

43 F = F'; 44 Des = temp(:); 45 Wt = (1./Dev)';

and

35 elseif itype==2 <sup>36</sup> Des(2\*k) = 0; <sup>37</sup> Dev(2\*k) = aid1;

42 temp = Des(ones(2,1),:);

constraints are, respectively,

The relevant MATLAB commands are

1 >> % Data for Type A design

5 >> % Data for Type B design

3 [.001 .01 .01],0.0005,1);

27 end 28 for k=1:R-1

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 2*R* 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 (17e), 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.

```
1 function [F,Des,Wt]=convert2constrt(F_ori, Des_ori,Dev_ori,alpha,itype)
2 % This function converts the original deisgn specifications into Type A
3 % or Type B transition band constraints compatible specifications.
4 %
5 % Input parameters:
6 % - F_ori contains the edges of the R bands, where R = length(Des_s).
7 % - Des_ori contains the desired values in the R bands.
8 % - Dev_ori contains the admissiable deviations in the R bands.
9 % - alpha is a small positive constant.
10 % - type=1 (2) generates Type A (B) transition band constraints.
11 %
12 % Output parameters
13 % - F contains the edges of the 2R-1 bands.
14 % - Des contains the desired values on all the edges of 2R-1 bands.
15 % - Wt contains the weights in the 2R-1 bands.
16
17 % Check if the input data is correct
18 if (alpha < 0)
19 error('alpha should be a small positive number.');
20 end
21 R = numel(Des_ori);
22
```
52 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>17</sup> Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design 53

```
23 % Generate the output parameters
24 for k=1:R
25 F(4*(k-1)+1)=F_ori(2*k-1); F(4*(k-1)+2)=F_ori(2*k);
26 Des(2*(k-1)+1)=Des_ori(k); Dev(2*(k-1)+1)=Dev_ori(k);
27 end
28 for k=1:R-1
29 F(4*(k-1)+3)=F_ori(2*k)+alpha;F(4*(k-1)+4)=F_ori(2*k+1)-alpha;
30 aid1=max(Des_ori(k)+Dev_ori(k),Des_ori(k+1)+Dev_ori(k+1));
31 aid2=min(Des_ori(k)-Dev_ori(k),Des_ori(k+1)-Dev_ori(k+1));
32 if itype==1
33 Des(2*k) = (aid1+aid2)/2;
34 Dev(2*k) = aid1-Des(2*k);
35 elseif itype==2
36 Des(2*k) = 0;
37 Dev(2*k) = aid1;
38 else
39 error('Type should be either 1 or 2.');
40 end
41 end
42 temp = Des(ones(2,1),:);
43 F = F';
44 Des = temp(:);
45 Wt = (1./Dev)';
```
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

Ω *<sup>κ</sup>* = [0, 0.2, 0.2005, 0.2495, 0.25, 0.6, 0.6005, 0.6995, 0.7, 1] .

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] *WA* = [1000, 1.9782, 100, 1.9608, 100],

and

16 Will-be-set-by-IN-TECH

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

corresponding lower and upper limits, whereas the admissible deviation is the difference

1, 2, . . . , *R* − 1 the desired values, denoted by *D<sup>κ</sup>*, and the admissible deviations, denoted by

<sup>0</sup> for *Type B*. (17h)

max{*D<sup>κ</sup>* + *δκ*, *Dκ*+<sup>1</sup> + *δκ*+1} − *D<sup>κ</sup>* for *Type A*

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 2*R* 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 (17e), 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

max{*D<sup>κ</sup>* <sup>+</sup> *δκ*, *<sup>D</sup>κ*+<sup>1</sup> <sup>+</sup> *δκ*+1} for *Type B*. (17i)

between the upper limit and the above-mentioned desired value. Hence, in Ω*<sup>T</sup>*

[ max{*D<sup>κ</sup>* + *δκ*, *Dκ*+<sup>1</sup> + *δκ*+1} + min{*D<sup>κ</sup>* − *δκ*, *Dκ*+<sup>1</sup> − *δκ*+1}]

1 function [F,Des,Wt]=convert2constrt(F\_ori, Des\_ori,Dev\_ori,alpha,itype) 2 % This function converts the original deisgn specifications into Type A 3 % or Type B transition band constraints compatible specifications.

6 % - F\_ori contains the edges of the R bands, where R = length(Des\_s).

8 % - Dev\_ori contains the admissiable deviations in the R bands.

10 % - type=1 (2) generates Type A (B) transition band constraints.

14 % - Des contains the desired values on all the edges of 2R-1 bands.

7 % - Des\_ori contains the desired values in the R bands.

19 error('alpha should be a small positive number.');

*<sup>κ</sup>* for both types is the average of the

2 for *Type A*

*<sup>κ</sup>* for *κ* =

stated upper limit in the nearest passband region. The desired value in the *κ*th transition band Ω*<sup>T</sup>*

> *δ <sup>κ</sup>* =

as Implementation I or II in this contribution.

9 % - alpha is a small positive constant.

17 % Check if the input data is correct

13 % - F contains the edges of the 2R-1 bands.

15 % - Wt contains the weights in the 2R-1 bands.

*δ*

and

4 %

11 %

16

22

20 end

5 % Input parameters:

12 % Output parameters

18 if (alpha < 0)

21 R = numel(Des\_ori);

*<sup>κ</sup>*, are as follows:

*D<sup>κ</sup>* =

```
DB = [0, 0, 1, 0, 0]
WB = [1000, 0.9901, 100, 0.9901, 100].
```
The relevant MATLAB commands are

```
1 >> % Data for Type A design
2 >> [F1,A1,W1] = convert2constrt([0 0.2 0.25 0.6 0.7 1],[0 1 0],...
3 [.001 .01 .01],0.0005,1);
4 >> N1 = 103; firr_coeff1 = firremez_imp1(N1,F1,A1,W1);
5 >> % Data for Type B design
6 >> [F2,A2,W2] = convert2constrt([0 0.2 0.25 0.6 0.7 1],[0 1 0],...
7 [.001 .01 .01],0.0005,2); firr_coeff2 = firremez_imp1(N1,F2,A2,W2);
```
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.

**Figure 4.** Magnitude response of the two bandpass filters of Example 3.

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

$$
\omega\_{p1} = 0.15\pi, \omega\_{s1} = 0.3\pi, \omega\_{s2} = 0.6\pi, \omega\_{p2} = 0.65\pi, \delta\_{p1} = \delta\_{p2} = 0.01, \text{ and } \delta\_{s} = 0.001.
$$

−100

implementation of the algorithm.

Magnitude (dB)

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

**Figure 6.** Magnitude response of the two bandstop filters of Example 4.

**Figure 5.** Magnitude response of the bandstop filter of Example 4.

Angular frequency

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

Angular frequency

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

1+δ p and −δ s

linear passband details

angular frequency 0 0.05π 0.1π 0.15π

0.99 1 1.01

> 1+δ p and −(1−δ p )

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

Transition band peak (linear scale)

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

0.2π 0.3π

0

50

−80

−60

−40

−20

Magnitude (dB)

0

20

40

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

```
1 >> [N,F,A,W] = firpmord([0.15 0.3 0.6 0.65],[1 0 1],[0.01 0.001 0.01]);
2 >> firr_coeff = firremez_imp1(N-4,F,A,W);
```
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

```
1 >> % Data for Type A design
2 >> [F1,A1,W1] = convert2constrt([0 0.15 0.3 0.6 0.65 1],[1 0 1],...
3 [0.01 0.001 0.01],0.0005,1);
4 >> N1 = 102; firr_coeff1 = firremez_imp1(N1,F1,A1,W1);
5 >> % Data for Type B design
6 >> [F2,A2,W2] = convert2constrt([0 0.15 0.3 0.6 0.65 1],[1 0 1],...
7 [0.01 0.001 0.01],0.0005,2); firr_coeff2 = firremez_imp1(N1,F2,A2,W2);
```
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 −*δ<sup>s</sup>* = −0.001. Furthermore, the response in the first passband remains equiripple.

**Figure 5.** Magnitude response of the bandstop filter of Example 4.

18 Will-be-set-by-IN-TECH

1+δ p and −δ s

1+δ p and −(1−δ p )

0.99 1 1.01

**Figure 4.** Magnitude response of the two bandpass filters of Example 3.

Angular frequency

*ωp*<sup>1</sup> = 0.15*π*, *ωs*<sup>1</sup> = 0.3*π*, *ωs*<sup>2</sup> = 0.6*π*, *ωp*<sup>2</sup> = 0.65*π*, *δp*<sup>1</sup> = *δp*<sup>2</sup> = 0.01, and *δ<sup>s</sup>* = 0.001.

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

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 −*δ<sup>s</sup>* =

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

1 >> [N,F,A,W] = firpmord([0.15 0.3 0.6 0.65],[1 0 1],[0.01 0.001 0.01]);

2 >> [F1,A1,W1] = convert2constrt([0 0.15 0.3 0.6 0.65 1],[1 0 1],...

6 >> [F2,A2,W2] = convert2constrt([0 0.15 0.3 0.6 0.65 1],[1 0 1],... 7 [0.01 0.001 0.01],0.0005,2); firr\_coeff2 = firremez\_imp1(N1,F2,A2,W2);

−0.001. Furthermore, the response in the first passband remains equiripple.

4 >> N1 = 102; firr\_coeff1 = firremez\_imp1(N1,F1,A1,W1);

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

0.3π 0.4π 0.5π 0.6π

linear passband details

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

2 >> firr\_coeff = firremez\_imp1(N-4,F,A,W);

1 >> % Data for Type A design

3 [0.01 0.001 0.01],0.0005,1);

5 >> % Data for Type B design

commands are

Magnitude (dB)

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.

**Figure 6.** Magnitude response of the two bandstop filters of Example 4.



*Condition 3*: The sign of **wei\_err** as a function of *k* alternates at the consecutive enteries of

*Step 1:* Find all the values of *k* in the range 1 ≤ *k* ≤ *ngrid*, where a local extremum of

*Step 2:* Extract those entries from **aid1**, where the absolute value of **wei\_err** is larger than or

*Step 3:* Split the range 1 ≤ *κ* ≤ *n***aid2**, where *n***aid2** is the length of **aid2**, into the sub-ranges

enteries the *nz* enteries to be included in the "real" vector **real** such that the largest absolute values of **wei\_err**(*k*), when regarded as a function of *k*, are retained subject to the condition that the maxima and minima alternate at the consecutive extracted enteries. If the length of

**real** . If *<sup>n</sup>*(**start**)

**real** if the absolute value of **wei\_err**(*k*) at *k* =

(*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*) for *<sup>m</sup>* <sup>=</sup> 1, 2, . . . , *<sup>n</sup>***aid2** in such a way that the signs of

**real** of length *nz***start** such that its *m*th entry is the value among

**real** generated at **Segment 2**, this main step extracts from its

**real** and no extraction is required. Otherwise, the extraction is

(**start**)

(**start**) **real** and *n*˜

(**start**)

**real** optional pairs such that the first pair consists of the

(**start**)

**real** is the length of

(*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*), at which <sup>|</sup>**wei\_err**(**aid2**(*κ*))<sup>|</sup> achieves its maximum

(*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*) are the same.

**real** − *nz* is an odd integer, then the first

(**start**)

**real** should be simultaneously

**real** − 1 consecutive entry pairs.

(**init**)

**real** (*<sup>k</sup>* <sup>=</sup> *<sup>n</sup>*(**start**)

(**start**)

**real** ). Go to the next

**real** , respectively. If

**real** , is zero, then set

**real** (1) is less

**wei\_err**(*k*) occurs. Store these values of *k* in the ascending order into the vector **aid1**.

**real** serves as a start-up vector for generating the "real" vector **real** at **Segment**

 (**start**) **real** .

*κ*(*low*)

 (**start**)

step.

*n*˜ (**start**)

*STEP C*: Since *n*˜

*Step 1:* Set

*Step 2:* If *nz*(**init**)

**real** =

This vector

(**start**)

*Step 4:* Generate the vector

value. Go to *Step 1* under **Segment 3**.

performed using the following steps:

(**start**)

discarded. There are altogether *n*˜

first and last entries of ˜

is the smallest. Go to *STEP C*.

(**init**) **real** =

(**init**)

**Segment 3:** Based on the vector

**real** is *nz*, then **real** ≡

*STEP A*: Denote the length of

(last) entry is discarded from

**aid2**(*κ*) for *κ*(*low*)

**3**. This segment is carried out using the following four steps:

equal to |*dev*|, and store these entries into the vector **aid2**.

**wei\_err**(**aid2**(*κ*)) in the *m*th sub-range as given by *κ*(*low*)

(**start**)

(**start**)

than or equal (larger) than the corresponding value at *k* =

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

**real** − *nz* = 0, then stop. Otherwise go to the next step.

This segment is carried out using the following seven steps:

**real** . If *nz*(**init**)

(**start**)

(**start**)

(**start**)

(**start**)

**real** by *<sup>n</sup>*(**start**)

**real** <sup>−</sup> *nz* is an even integer, two entries of ˜

**real** and the remaining ones are *n*˜

**real** <sup>−</sup> *nz*, where *nz*(**init**)

**real** and go to *Step 1* under **Segment 4**. Otherwise, go to the next step.

**real** − *nz* is an even integer, then go to *Step 4*. Otherwise, go to the next step.

The pair to be discarded is the pair, where the maximum of two absolute values of **wei\_err**(*k*)

(**start**)

(**start**)

\*Time in seconds

�Transition bands are excluded.

<sup>⊥</sup>*<sup>A</sup> Type A* transition band constraints are used.

<sup>⊥</sup>*<sup>B</sup> Type B* transition band constraints are used.

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

*Condition 1*: At each entry of (**start**) **real** , **wei\_err**(*k*), when regarded as a function of *k*, achieves a local extremum whose absolute value is larger than or equal to |*dev*|, where |*dev*| is determined according to (6).

*Condition 2*: In case of several consecutive local extrema of **wei\_err**(*k*) with the same sign for *<sup>k</sup>*(*low*) <sup>≤</sup> *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*(*upp*) , only one entry is included in (**start**) **real** and its value is the value of *k* for *<sup>k</sup>*(*low*) <sup>≤</sup> *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*(*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 of (**start**) **real** .

20 Will-be-set-by-IN-TECH

Example Original Implementation Implementation I

 0.152 0.125 0.032 0.032 0.235 0.136 0.064 0.064 � 0.159 0.120 0.032 0.032 ⊥*<sup>A</sup>* 0.262 0.191 0.056 0.056 ⊥*<sup>B</sup>* 0.307 0.184 0.061 0.061 � 0.198 0.138 0.047 0.047 ⊥*<sup>A</sup>* 0.272 0.169 0.051 0.051 ⊥*<sup>B</sup>* 0.318 0.186 0.055 0.055

\*Time in seconds

**4. Implementation II**

**4.1. Theoretical formulation**

main step generates the vector

*Condition 1*: At each entry of

according to (6).

*<sup>k</sup>*(*low*) <sup>≤</sup> *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*(*upp*)

*<sup>k</sup>*(*low*) <sup>≤</sup> *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*(*upp*)

with the accompanying steps are as follows.

order subject to the following three conditions:

�Transition bands are excluded.

<sup>⊥</sup>*<sup>A</sup> Type A* transition band constraints are used. <sup>⊥</sup>*<sup>B</sup> Type B* transition band constraints are used.

same optimum solution, especially in multiband cases.

**Table 1.** Performance Comparison of Original Implementation and Implementation I

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

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

**Segment 2:** Based on the values of **wei\_err**(*k*) for 1 ≤ *k* ≤ *ngrid* generated at **Segment 1**, this

local extremum whose absolute value is larger than or equal to |*dev*|, where |*dev*| is determined

*Condition 2*: In case of several consecutive local extrema of **wei\_err**(*k*) with the same sign for

, where |**wei\_err**(*k*)| achieves its maximum value.

**real** to include as many entries as possible in the ascending

**real** , **wei\_err**(*k*), when regarded as a function of *k*, achieves a

**real** and its value is the value of *k* for

(**start**)

(**start**)

, only one entry is included in

(**start**)

Total Time\* Self Time Total Time Self Time

This vector (**start**) **real** 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 range 1 ≤ *k* ≤ *ngrid*, where a local extremum of **wei\_err**(*k*) occurs. Store these values of *k* in the ascending order into the vector **aid1**.

*Step 2:* Extract those entries from **aid1**, where the absolute value of **wei\_err** is larger than or equal to |*dev*|, and store these entries into the vector **aid2**.

*Step 3:* Split the range 1 ≤ *κ* ≤ *n***aid2**, where *n***aid2** is the length of **aid2**, into the sub-ranges *κ*(*low*) (*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*) for *<sup>m</sup>* <sup>=</sup> 1, 2, . . . , *<sup>n</sup>***aid2** in such a way that the signs of **wei\_err**(**aid2**(*κ*)) in the *m*th sub-range as given by *κ*(*low*) (*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*) are the same.

*Step 4:* Generate the vector (**start**) **real** of length *nz***start** such that its *m*th entry is the value among **aid2**(*κ*) for *κ*(*low*) (*m*) <sup>≤</sup> *<sup>κ</sup>* <sup>≤</sup> *<sup>κ</sup>*(*upp*)(*m*), at which <sup>|</sup>**wei\_err**(**aid2**(*κ*))<sup>|</sup> achieves its maximum value. Go to *Step 1* under **Segment 3**.

**Segment 3:** Based on the vector (**start**) **real** 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 of **wei\_err**(*k*), when regarded as a function of *k*, are retained subject to the condition that the maxima and minima alternate at the consecutive extracted enteries. If the length of (**start**) **real** is *nz*, then **real** ≡ (**start**) **real** and no extraction is required. Otherwise, the extraction is performed using the following steps:

*STEP A*: Denote the length of (**start**) **real** by *<sup>n</sup>*(**start**) **real** . If *<sup>n</sup>*(**start**) **real** − *nz* is an odd integer, then the first (last) entry is discarded from (**start**) **real** if the absolute value of **wei\_err**(*k*) at *k* = (**start**) **real** (1) is less than or equal (larger) than the corresponding value at *k* = (**start**) **real** (*<sup>k</sup>* <sup>=</sup> *<sup>n</sup>*(**start**) **real** ). Go to the next step.

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

*STEP C*: Since *n*˜ (**start**) **real** <sup>−</sup> *nz* is an even integer, two entries of ˜ (**start**) **real** should be simultaneously discarded. There are altogether *n*˜ (**start**) **real** optional pairs such that the first pair consists of the first and last entries of ˜ (**start**) **real** and the remaining ones are *n*˜ (**start**) **real** − 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:* Set (**init**) **real** = (**start**) **real** . If *nz*(**init**) **real** <sup>−</sup> *nz*, where *nz*(**init**) **real** is the length of (**init**) **real** , is zero, then set **real** = (**init**) **real** and go to *Step 1* under **Segment 4**. Otherwise, go to the next step.

*Step 2:* If *nz*(**init**) **real** − *nz* is an even integer, then go to *Step 4*. Otherwise, go to the next step. *Step 3:* If |**wei\_err**( (**init**) **real** (1))|≤|**wei\_err**( (**init**) **real** (*nz*(**init**) **real** ))|, then discard the first entry from (**init**) **real** . Otherwise, discard the last entry from (**init**) **real** . Go to the next step.

*Step 4:* If *nz*(**init**) **real** <sup>−</sup> *nz*, where *nz*(**init**) **real** is the length of the remaining vector (**init**) **real** , is zero, then set **real** = (**init**) **real** 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:

$$\begin{split} \mathtt{wei\\_comp}(m) &= \max\left( |\mathtt{wei\\_err}(\ell^{(\mathtt{init})}\_{\mathtt{real}}(m))|, |\mathtt{wei\\_err}(\ell^{(\mathtt{init})}\_{\mathtt{real}}(m+1))| \right) \\ &\qquad \text{for } m = 1, 2, \ldots, n\tau^{(\mathtt{init})}\_{\mathtt{real}} - 1. \end{split} \tag{18}$$

 add(2:2:nz) = -add(2:2:nz); dnum = ad\*des(l\_trial)'; dden = add\*(ad./wt(l\_trial))';

 % Step 2: Lagrange ordinate vector y y = des(l\_trial) + dev\*add./wt(l\_trial); % Step 3: Overall abscissa vector x\_all x\_all = cos(2\*pi\*grid(1:ngrid));

err\_den = err\_num; % and err\_den.

err\_num = err\_num + y(jj)\*aid;

% STEP II DETERMINE THE VECTOR l\_real\_start

err\_cy = err\_num./err\_den;

wei\_err(l\_trial)=dev\_vect;

% Step 2: Determine l\_aid2.

abs(wei\_err(l\_aid2))));

else % of l\_real\_init.

 % STEP III DETERMINE THE VECTOR l\_real l\_real\_init = l\_real\_start; % Step 1

% Step 1: Find l\_aid1.

end

 end end

else

 end end

l\_real = l\_real\_init;

% STEP IV: TEST CONVERGENCE

dev = -dnum/dden; % Step 1: Current value of deviation.

% by alternatingly using the values of dev and -dev.

l\_aid1 = find(diff(sign(diff([0 wei\_err 0]))));

 cumsum([1,(wei\_err(l\_aid2(2:end))>=0) ... ~=(wei\_err(l\_aid2(1:end-1))>=0)]),...

while numel(l\_real\_init) > nz % Step 4

wei\_comp=max(wei\_real(1:end-1),...

 l\_aid2 = l\_aid1(abs(wei\_err(l\_aid1)) >= abs(dev)); [~,ind] = max(sparse(1:length(l\_aid2),... % Step 3

l\_real\_start = l\_aid2(ind); % Step 4: Determine l\_real\_start.

 if rem(numel(l\_real\_init) - nz,2) == 1 % Step 2: odd difference. if abs(wei\_err(l\_real\_init(1))) <= abs(wei\_err(l\_real\_init(end))) l\_real\_init(1) = []; % Step 3: discard the first entry...

l\_real\_init(end) = []; % otherwise discard the last entry.

 wei\_real(2:end)); % End of Step 5 if max(abs(wei\_err(l\_real\_init(1))),... % Start of Step 6

l\_real\_init = l\_real\_init(2:end-1); % End of Step 6

 [~,ind\_omit]=min(wei\_comp); % Start: Step 7 l\_real\_init(ind\_omit:ind\_omit+1) = []; % End: Step 7

 if (l\_real == l\_trial) % Step 1: The real and trial vectors... break; % coincide. Hence, stop. Remez loop ended successfully.

wei\_real=abs(wei\_err(l\_real\_init)); % Start of Step 5

abs(wei\_err(l\_real\_init(end))))<= min(wei\_comp)

err\_num = zeros(1,ngrid); % Step 4: Initializations of err\_num...

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

 wei\_err = (err\_cy - des).\*wt; % Step 7: Generate the vector wei\_err. dev\_vect = ones(size(l\_trial)); % Steps 8-10: Fill in the undefined dev\_vect(2:2:length(l\_trial))= -dev\_vect(2:2:length(l\_trial)); dev\_vect = dev\_vect \* dev; % entries of wei\_err at l\_trial(1:nz)...

 for jj = 1:nz % Steps 5 and 6: Intermediate evaluations for... aid = ad(jj)./(x\_all - x(jj)); % obtaining the weighted error... err\_den = err\_den + aid; % wei\_err(k) at all the grid points.

*Step 6:* If max |**wei\_err**( (**init**) **real** (1))|, |**wei\_err**( (**init**) **real** (*nz*(**init**) **real** ))| is less than or equal to the largest entry of **wei\_comp**, then remove the first and last entries from (**init**) **real** 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 *m*th entry, then remove the *m*th and the (*m* + 1)th entries from (**init**) **real** . 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.

```
1 function [x,y,ad,dev] = remez_imp2(nz,iext,ngrid,grid,des,wt)
2 % remez_imp1 implements the Segments 1 - 4 described in the preceding
3 % section, the function needs to be inserted within the MATLAB function
4 % firpm. The input argument values come directly from the function firpm
5 % and the output arguments are required to perform the Inverse Fourier
6 % transform in order to calculate the filter coefficients. In case of
7 % any issues send an e-mail to muhammad"dot"ahsan"at"tut "dot" fi.
8 % Last updated 04.15.2012 4:54 AM (UTC/GMT+2)
9
10 % INITIALIZATIONS PHASE
11 niter = 1; % Initialize the iteration counter.
12 itrmax = 250; % Maximum number of iterations.
13 l_trial = iext(1:nz)'; % Startup value of l_trial.
14
15 % ITERATION PHASE
16 % REMEZ LOOP FOR LOCATING DESIRED nz INDICES AMONG THE GRID POINTS
17 while (niter < itrmax)
18
19 % STEP I: BASED ON THE PRESENT 'TRIAL' VECTOR l_trial, GENERATE THE
20 % WEIGHTED ERROR FUNCTION wei_err(k) AT ALL THE GRID POINTS
21 x = cos(2*pi*grid(l_trial)); % Step 1: Lagrange abscissa vector x.
22 A = x'*ones(1,nz)-ones(nz,1)*x;
23 A(eye(nz)==1) = 1;
24 ad = prod(A);
25 ad = ad * (-2)^(nz-1); % Step 1: Lagrange coefficient vector ad...
26 ad = 1./ad; % found efficiently without using the function remezdd.
27 add = ones(size(ad));
```
 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>23</sup> Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design 59

Will-be-set-by-IN-TECH

(**init**) **real** (*nz*(**init**)

(**init**)

**real** and go to *Step 1* under **Segment 4**. Otherwise, go to the next step.

(**init**)

*Step 7:* Find the entry of **wei\_comp** with the smallest value. If this is the *m*th entry, then

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

(**init**) **real** (*nz*(**init**)

(**init**)

for *<sup>m</sup>* = 1, 2, . . . , *nz*(**init**)

**real** is the length of the remaining vector

**real** (*m*))|, |**wei\_err**(

**real** − 1.

**real** ))| 

**real** . Go to *Step 4*.

**real** . Go to the next step.

**real** ))|, then discard the first entry from

(**init**)

**real** (*m* + 1))|

(**init**)

(**init**)

**real** , is zero, then

(18)

is less than or equal to the

**real** and go to *Step 4*.

*Step 3:* If |**wei\_err**(

(**init**)

Otherwise, go to the next step.

**4.2. MATLAB code snippet**

% INITIALIZATIONS PHASE

% ITERATION PHASE

while (niter < itrmax)

 A(eye(nz)==1) = 1; ad = prod(A);

add = ones(size(ad));

A = x'\*ones(1,nz)-ones(nz,1)\*x;

*Step 4:* If *nz*(**init**)

set **real** =

*Step 6:* If max

Implementation I.

 (**init**) (**init**)

**real** . Otherwise, discard the last entry from

**real** <sup>−</sup> *nz*, where *nz*(**init**)

**wei\_comp**(*m*) =max


remove the *m*th and the (*m* + 1)th entries from

% Last updated 04.15.2012 4:54 AM (UTC/GMT+2)

 niter = 1; % Initialize the iteration counter. itrmax = 250; % Maximum number of iterations. l\_trial = iext(1:nz)'; % Startup value of l\_trial.

% REMEZ LOOP FOR LOCATING DESIRED nz INDICES AMONG THE GRID POINTS

 % STEP I: BASED ON THE PRESENT 'TRIAL' VECTOR l\_trial, GENERATE THE % WEIGHTED ERROR FUNCTION wei\_err(k) AT ALL THE GRID POINTS

x = cos(2\*pi\*grid(l\_trial)); % Step 1: Lagrange abscissa vector x.

 ad = ad \* (-2)^(nz-1); % Step 1: Lagrange coefficient vector ad... ad = 1./ad; % found efficiently without using the function remezdd.

**real** (1))|≤|**wei\_err**(

*Step 5:* Determine the entries of the vector **wei\_comp** as follows:

(**init**)


**real** (1))|, |**wei\_err**(

largest entry of **wei\_comp**, then remove the first and last entries from

function [x,y,ad,dev] = remez\_imp2(nz,iext,ngrid,grid,des,wt)

 % remez\_imp1 implements the Segments 1 - 4 described in the preceding % section, the function needs to be inserted within the MATLAB function % firpm. The input argument values come directly from the function firpm % and the output arguments are required to perform the Inverse Fourier % transform in order to calculate the filter coefficients. In case of % any issues send an e-mail to muhammad"dot"ahsan"at"tut "dot" fi.

```
28 add(2:2:nz) = -add(2:2:nz);
29 dnum = ad*des(l_trial)';
30 dden = add*(ad./wt(l_trial))';
31 dev = -dnum/dden; % Step 1: Current value of deviation.
32 % Step 2: Lagrange ordinate vector y
33 y = des(l_trial) + dev*add./wt(l_trial);
34 % Step 3: Overall abscissa vector x_all
35 x_all = cos(2*pi*grid(1:ngrid));
36 err_num = zeros(1,ngrid); % Step 4: Initializations of err_num...
37 err_den = err_num; % and err_den.
38 for jj = 1:nz % Steps 5 and 6: Intermediate evaluations for...
39 aid = ad(jj)./(x_all - x(jj)); % obtaining the weighted error...
40 err_den = err_den + aid; % wei_err(k) at all the grid points.
41 err_num = err_num + y(jj)*aid;
42 end
43 err_cy = err_num./err_den;
44 wei_err = (err_cy - des).*wt; % Step 7: Generate the vector wei_err.
45 dev_vect = ones(size(l_trial)); % Steps 8-10: Fill in the undefined
46 dev_vect(2:2:length(l_trial))= -dev_vect(2:2:length(l_trial));
47 dev_vect = dev_vect * dev; % entries of wei_err at l_trial(1:nz)...
48 % by alternatingly using the values of dev and -dev.
49 wei_err(l_trial)=dev_vect;
50
51 % STEP II DETERMINE THE VECTOR l_real_start
52 % Step 1: Find l_aid1.
53 l_aid1 = find(diff(sign(diff([0 wei_err 0]))));
54 % Step 2: Determine l_aid2.
55 l_aid2 = l_aid1(abs(wei_err(l_aid1)) >= abs(dev));
56 [~,ind] = max(sparse(1:length(l_aid2),... % Step 3
57 cumsum([1,(wei_err(l_aid2(2:end))>=0) ...
58 ~=(wei_err(l_aid2(1:end-1))>=0)]),...
59 abs(wei_err(l_aid2))));
60 l_real_start = l_aid2(ind); % Step 4: Determine l_real_start.
61
62 % STEP III DETERMINE THE VECTOR l_real
63 l_real_init = l_real_start; % Step 1
64 if rem(numel(l_real_init) - nz,2) == 1 % Step 2: odd difference.
65 if abs(wei_err(l_real_init(1))) <= abs(wei_err(l_real_init(end)))
66 l_real_init(1) = []; % Step 3: discard the first entry...
67 else % of l_real_init.
68 l_real_init(end) = []; % otherwise discard the last entry.
69 end
70 end
71 while numel(l_real_init) > nz % Step 4
72 wei_real=abs(wei_err(l_real_init)); % Start of Step 5
73 wei_comp=max(wei_real(1:end-1),...
74 wei_real(2:end)); % End of Step 5
75 if max(abs(wei_err(l_real_init(1))),... % Start of Step 6
76 abs(wei_err(l_real_init(end))))<= min(wei_comp)
77 l_real_init = l_real_init(2:end-1); % End of Step 6
78 else
79 [~,ind_omit]=min(wei_comp); % Start: Step 7
80 l_real_init(ind_omit:ind_omit+1) = []; % End: Step 7
81 end
82 end
83 l_real = l_real_init;
84
85 % STEP IV: TEST CONVERGENCE
86 if (l_real == l_trial) % Step 1: The real and trial vectors...
87 break; % coincide. Hence, stop. Remez loop ended successfully.
```
24 Will-be-set-by-IN-TECH 60 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>25</sup>

```
88 else
89 l_trial = l_real; % Step 2: Otherwise, replace the values of...
90 niter = niter + 1; % l_trial with the values of l_real and...
91 end % continue.
92 end % END OF THE OVERALL REMEZ LOOP
```
## **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:

−100

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

**Figure 8.** Magnitude response of the five-band filter of Example 6.

9 by using a dashed red line and a dot-dashed black line, respectively.

Magnitude (dB)

**Figure 7.** Magnitude response of the five-band filter of Example 5.

Angular frequency

Passband details (linear)

0.4π 0.5π 0.6π 0.7π

Angular frequency

The minimum filter order required to meet these specifications is 100. 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 is 101. The responses of the resulting filters meeting these additional constraints are depicted in Fig.

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

0.99 0.995 1 1.005 1.01

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

First passband (linear scale)

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

0.25π 0.3π 0.35π 0.4π 0.45π

−80

−60

−40

−20

Magnitude (dB)

0

0.99 1 1.01

20

40

$$
\omega\_{s1} = 0.17\pi, \omega\_{p1} = 0.23\pi, \omega\_{p2} = 0.47\pi, \omega\_{s2} = 0.53\pi, \omega\_{s3} = 0.67\pi, \omega\_{p3} = 0.73\pi,
$$

$$
\omega\_{p4} = 0.82\pi, \omega\_{s4} = 0.88\pi, \delta\_{s1} = \delta\_{s2} = \delta\_{s3} = 0.001, \text{ and } \delta\_{p1} = \delta\_{p2} = 0.01.
$$

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

1 >> [n,f,a,w] = firpmord([.17 .23 .47 .53 .67 .73 .82 .88],... 2 [0 1 0 1 0], [.001 .01 .001 .01 .001]); 3 >> firr\_coeff = firremez\_imp2(n+4,f,a,w);

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:

$$
\omega\_{p1} = 0.1\pi, \omega\_{s1} = 0.15\pi, \omega\_{s2} = 0.3\pi, \omega\_{p2} = 0.35\pi, \omega\_{p3} = 0.75\pi, \omega\_{s3} = 0.8\pi, \omega\_{s4} = 0.85\pi,
$$

$$
\omega\_{p4} = 0.9\pi, \delta\_{p1} = \delta\_{p2} = \delta\_{p3} = 0.01, \text{ and } \delta\_{s1} = \delta\_{s2} = 0.001.
$$

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

1 >> [n,f,a,w] = firpmord([.1 .15 .3 .35 .75 .8 .85 .9],[1 0 1 0 1],... 2 [.01 .001 .01 .001 .01]); firr\_coeff = firremez\_imp2(n,f,a,w);

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:

$$
\omega\_{s1} = 0.15\pi, \omega\_{p1} = 0.2\pi, \omega\_{p2} = 0.45\pi, \omega\_{s2} = 0.55\pi, \omega\_{s3} = 0.7\pi, \omega\_{p3} = 0.8\pi, \omega\_{p4} = 0.85\pi,
$$

$$
\omega\_{s4} = 0.93\pi, \delta\_{s1} = \delta\_{s2} = \delta\_{s3} = 0.001, \text{ and } \delta\_{p1} = \delta\_{p2} = 0.01.
$$

**Figure 7.** Magnitude response of the five-band filter of Example 5.

24 Will-be-set-by-IN-TECH

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

**Example 5:** It is desired to design a five-band filter with two passbands and three stopbands

*ωs*<sup>1</sup> = 0.17*π*, *ωp*<sup>1</sup> = 0.23*π*, *ωp*<sup>2</sup> = 0.47*π*, *ωs*<sup>2</sup> = 0.53*π*, *ωs*<sup>3</sup> = 0.67*π*, *ωp*<sup>3</sup> = 0.73*π*, *ωp*<sup>4</sup> = 0.82*π*, *ωs*<sup>4</sup> = 0.88*π*, *δs*<sup>1</sup> = *δs*<sup>2</sup> = *δs*<sup>3</sup> = 0.001, and *δp*<sup>1</sup> = *δp*<sup>2</sup> = 0.01.

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

*ωp*<sup>1</sup> = 0.1*π*, *ωs*<sup>1</sup> = 0.15*π*, *ωs*<sup>2</sup> = 0.3*π*, *ωp*<sup>2</sup> = 0.35*π*, *ωp*<sup>3</sup> = 0.75*π*, *ωs*<sup>3</sup> = 0.8*π*, *ωs*<sup>4</sup> = 0.85*π*, *ωp*<sup>4</sup> = 0.9*π*, *δp*<sup>1</sup> = *δp*<sup>2</sup> = *δp*<sup>3</sup> = 0.01, and *δs*<sup>1</sup> = *δs*<sup>2</sup> = 0.001.

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

*ωs*<sup>1</sup> = 0.15*π*, *ωp*<sup>1</sup> = 0.2*π*, *ωp*<sup>2</sup> = 0.45*π*, *ωs*<sup>2</sup> = 0.55*π*, *ωs*<sup>3</sup> = 0.7*π*, *ωp*<sup>3</sup> = 0.8*π*, *ωp*<sup>4</sup> = 0.85*π*, *ωs*<sup>4</sup> = 0.93*π*, *δs*<sup>1</sup> = *δs*<sup>2</sup> = *δs*<sup>3</sup> = 0.001, and *δp*<sup>1</sup> = *δp*<sup>2</sup> = 0.01.

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

1 >> [n,f,a,w] = firpmord([.1 .15 .3 .35 .75 .8 .85 .9],[1 0 1 0 1],... 2 [.01 .001 .01 .001 .01]); firr\_coeff = firremez\_imp2(n,f,a,w);

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

1 >> [n,f,a,w] = firpmord([.17 .23 .47 .53 .67 .73 .82 .88],...

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

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

89 l\_trial = l\_real; % Step 2: Otherwise, replace the values of... 90 niter = niter + 1; % l\_trial with the values of l\_real and...

88 else

91 end % continue.

92 end % END OF THE OVERALL REMEZ LOOP

constraints described in Subsection 3.3.

meeting the following specifications:

2 [0 1 0 1 0], [.001 .01 .001 .01 .001]); 3 >> firr\_coeff = firremez\_imp2(n+4,f,a,w);

with following specifications:

with following specifications:

**4.3. Performance comparison**

**Figure 8.** Magnitude response of the five-band filter of Example 6.

The minimum filter order required to meet these specifications is 100. 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 is 101. 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.

26 Will-be-set-by-IN-TECH 62 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>27</sup>

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

1 >> % Filter design without transition band constraints

5 >> % Filter design with Type A transition band constraints

9 >> % Filter design with Type B transition band constraints

3 [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01]);

8 >> N1 = 104; htbi1 = firremez\_imp2(N1,F1,A1,W1);

12 >> N2 = 102; htbi2 = firremez\_imp2(N2,F2,A2,W2);

−100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 10

\*Time in seconds

�Transition bands are excluded.

<sup>⊥</sup>*AType A* transition band constraints are used. <sup>⊥</sup>*BType B* transition band constraints are used.

**Figure 10.** Magnitude response of three five-band filters of Example 8.

Magnitude (dB)

4 >> h = firremez\_imp2(n-4,f,a,w);

2 >> [n,f,a,w] = firpmord([0.17 0.27 0.47 0.52 0.69 0.79 0.87 0.92],...

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

6 >> [F1,A1,W1] = convert2constrt([0 0.17 0.27 0.47 0.52 0.69 0.79,... 7 0.87 0.92 1], [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01],0.0005,1);

10 >> [F2,A2,W2] = convert2constrt([0 0.17 0.27 0.47 0.52 0.69 0.79,... 11 0.87 0.92 1], [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01],0.0005,2);

Angular frequency

Example Original Implementation Implementation II Reduction Percentage(%)

 34 0.223 7 0.024 79 89 35 0.269 16 0.040 54 85 � 15 0.160 15 0.033 − 79 ⊥*<sup>A</sup>* 93 0.617 23 0.055 75 91 ⊥*<sup>B</sup>* 38 0.283 20 0.051 47 82 � 11 0.146 11 0.029 − 80 ⊥*<sup>A</sup>* 49 0.360 37 0.082 24 77 ⊥*<sup>B</sup>* 66 0.398 23 0.057 65 86

**Table 2.** Performance Comparison of Original Implementation and Implementation II.

Iterations Time\* Iterations Time Iterations Time

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

No constraint 1+δ p and −δ s

1+δ p and −(1−δ p )

```
1 >> % Filter design without transition band constraints
2 >> [n,f,a,w] = firpmord([.15 .2 .45 .55 .7 .8 .85 .93],...
3 [0 1 0 1 0],[.001 .01 .001 .01 .001]);
4 >> N1 = 101;h = firremez_imp2(n-2,f,a,w);
5 >> % Filter design with Type A transition band constraints
6 >> [F1,A1,W1] = convert2constrt([0 .15 .2 .45 .55 .7 .8 .85 .93 1],...
7 [0 1 0 1 0],[.001 .01 .001 .01 .001],0.0005,1);
8 >> htbi1 = firremez_imp2(N1,F1,A1,W1);
9 >> % Filter design with Type B transition band constraints
10 >> [F2,A2,W2] = convert2constrt([0 .15 .2 .45 .55 .7 .8 .85 .93 1],...
11 [0 1 0 1 0],[.001 .01 .001 .01 .001],0.0005,2);
12 >> htbi2 = firremez_imp2(N1,F2,A2,W2);
```
**Figure 9.** Magnitude responses of the three five-band filters of Example 7.

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

$$
\omega\_{p1} = 0.17\pi, \omega\_{s1} = 0.27\pi, \omega\_{s2} = 0.47\pi, \omega\_{p2} = 0.52\pi, \omega\_{p3} = 0.69\pi, \omega\_{s3} = 0.79\pi, \omega\_{s4} = 0.87\pi,
$$

$$
\omega\_{p4} = 0.92\pi, \delta\_{p1} = \delta\_{p2} = \delta\_{p3} = 0.01, \text{ and } \delta\_{s1} = \delta\_{s2} = 0.001.
$$

The minimum filter order to meet these specifications is 102. 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, 2 dB, 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 and 102, 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

62 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design <sup>27</sup> Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design 63

```
1 >> % Filter design without transition band constraints
2 >> [n,f,a,w] = firpmord([0.17 0.27 0.47 0.52 0.69 0.79 0.87 0.92],...
3 [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01]);
4 >> h = firremez_imp2(n-4,f,a,w);
5 >> % Filter design with Type A transition band constraints
6 >> [F1,A1,W1] = convert2constrt([0 0.17 0.27 0.47 0.52 0.69 0.79,...
7 0.87 0.92 1], [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01],0.0005,1);
8 >> N1 = 104; htbi1 = firremez_imp2(N1,F1,A1,W1);
9 >> % Filter design with Type B transition band constraints
10 >> [F2,A2,W2] = convert2constrt([0 0.17 0.27 0.47 0.52 0.69 0.79,...
11 0.87 0.92 1], [1,0,1,0,1],[0.01,0.001,0.01,0.001,0.01],0.0005,2);
12 >> N2 = 102; htbi2 = firremez_imp2(N2,F2,A2,W2);
```
**Figure 10.** Magnitude response of three five-band filters of Example 8.


\*Time in seconds

26 Will-be-set-by-IN-TECH

Angular frequency

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

*ωp*<sup>1</sup> = 0.17*π*, *ωs*<sup>1</sup> = 0.27*π*, *ωs*<sup>2</sup> = 0.47*π*, *ωp*<sup>2</sup> = 0.52*π*, *ωp*<sup>3</sup> = 0.69*π*, *ωs*<sup>3</sup> = 0.79*π*, *ωs*<sup>4</sup> = 0.87*π*, *ωp*<sup>4</sup> = 0.92*π*, *δp*<sup>1</sup> = *δp*<sup>2</sup> = *δp*<sup>3</sup> = 0.01, and *δs*<sup>1</sup> = *δs*<sup>2</sup> = 0.001.

The minimum filter order to meet these specifications is 102. 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, 2 dB, 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 and 102, respectively, whereas the responses of the resulting filters are depicted in Fig. 10 using a dashed red line

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

0 0.1π 0.2π 0.3π 0.4π 0.5π 0.6π 0.7π 0.8π 0.9π π

No constraints 1+δ p and −δ s

1+δ p and −(1−δ p )

**Figure 9.** Magnitude responses of the three five-band filters of Example 7.

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

1 >> % Filter design without transition band constraints 2 >> [n,f,a,w] = firpmord([.15 .2 .45 .55 .7 .8 .85 .93],...

7 [0 1 0 1 0],[.001 .01 .001 .01 .001],0.0005,1);

11 [0 1 0 1 0],[.001 .01 .001 .01 .001],0.0005,2);

5 >> % Filter design with Type A transition band constraints

9 >> % Filter design with Type B transition band constraints

6 >> [F1,A1,W1] = convert2constrt([0 .15 .2 .45 .55 .7 .8 .85 .93 1],...

10 >> [F2,A2,W2] = convert2constrt([0 .15 .2 .45 .55 .7 .8 .85 .93 1],...

3 [0 1 0 1 0],[.001 .01 .001 .01 .001]); 4 >> N1 = 101;h = firremez\_imp2(n-2,f,a,w);

8 >> htbi1 = firremez\_imp2(N1,F1,A1,W1);

12 >> htbi2 = firremez\_imp2(N1,F2,A2,W2);

−100

and a dot-dashed black line, respectively.

with following specifications:

−80

−60

−40

Magnitude (dB)

−20

0

20

�Transition bands are excluded.

<sup>⊥</sup>*AType A* transition band constraints are used.

<sup>⊥</sup>*BType B* transition band constraints are used.

**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.

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,

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

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

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

[1] Ahsan, M. (2008). Design 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,

[2] Ahsan, M. & Saramäki, T. (2009). Significant improvements in translating the Parks-McClellan algorithm from its FORTRAN code to its corresponding MATLAB

[3] Ahsan, M. & Saramäki, T. (2011). A MATLAB Based Optimum Multiband FIR Filters Design Program Following the Original Idea of the Remez Multiple Exchange Algorithm, In:*IEEE Symp. Circuits Syst.*, (May 2011), Rio de Janiero Brazil, pp 137-140. [4] Cheney, W. E. (1966). *Introduction to Approximation Theory*, AMS Chalsea Publishing. [5] McClellan, J. H. & Parks, T. W. (1972). A program for the design of linear phase finite impulse response digital filters, In: *IEEE Trans. Audio Electroacoust.*, Vol. AU-20, (Aug.

[6] McClellan, J. H. & Parks, T. W. (1973). A unified approach to the Design of Optimum FIR linear phase digital filters, In: *IEEE Trans. Circuit Theory*, Vol. 20, (Nov. 1973) pp.

[7] McClellan, J. H. & Parks, T. W. & Rabiner, L. R. (1973). A computer program for designing optimum FIR linear phase digital filters, In: *IEEE Trans. Audio Electroacoust.*,

[8] Novodvorskii, E. P. & Pinsker, I. S. (1951). The process of equating maxima, In: *Uspekhi*

[9] Rabiner, L. R., Kaiser, J. F. & Schafer, R. W. (1974). Some considerations in the design of multiband finite-impulse-response digital filters, In: *IEEE Trans. Acoust. Speech, Signal*

[10] Rabiner, R. L., Gold, G. (1975). *Theory and Application of DIGITAL SIGNAL PROCESSING*,

*Mat. Nauk (USSR)*, Vol. 6, (1951) pp. 174-181 (Engl. transl. by A. Schenitzer).

code, In:*IEEE Symp. Circuits Syst.*, (May 2009), Taipei, Taiwan, pp. 289-292.

whereas the reduction in the CPU execution time is around **80** percent.

time.

function efficiently.

**Author details**

**6. References**

Muhammad Ahsan and Tapio Saramäki

Finland, 107 pages.

1972) pp. 195-199.

Vol. 21, (Dec. 1973) pp. 506-526.

*Processing*, Vol. 22, (Dec. 1974) pp. 462-472.

697-701.

Prentice Hall.

*Tampere University of Technology, Tampere, Finland*

## **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.

## **Author details**

28 Will-be-set-by-IN-TECH

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

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

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

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

implementation is highly efficient in the design of multiband filters.

**5. Conclusion**

implementation.

points under consideration.

Muhammad Ahsan and Tapio Saramäki *Tampere University of Technology, Tampere, Finland*

## **6. References**

	- [11] Remez, E. (1934). Sur le calcul effectifdes polynomes d'approximation de Tchebychef, In: *Compt. Rend. Acad. Sci*, Vol. 199, (1934) pp. 337-340.
	- [12] Rice, R. J. (1964). *The Approximation of Functions. Volume 1: Linear Theory*, Addison-Wesley Pub (Sd).

**MATLAB/SIMULINK and Its Engineering** 

**Applications** 


**MATLAB/SIMULINK and Its Engineering Applications** 

30 Will-be-set-by-IN-TECH

[11] Remez, E. (1934). Sur le calcul effectifdes polynomes d'approximation de Tchebychef,

[12] Rice, R. J. (1964). *The Approximation of Functions. Volume 1: Linear Theory*,

[13] Rivlin, J. T. (2010). *An Introduction to the Approximation of Functions*, Dover Publications. [14] Saramäki, T. (1981). Design of digital filters requiring a small number of arithmetic operations, In:*Dept. of Electrical Engineering, Tampere University of Technology*, Dr. Tech.

[15] Saramäki, T. (1987). Efficient iterative algorithms for the design of optimum IIR filters with arbitrary specifications, In:*Proc. Int. Conf. Digital Signal Process.*, Florence,Italy, (Sep.

[16] Saramäki, T. (1992). An efficient Remez-type algorithm for the design of optimum IIR filters with arbitrary partially constrained specifications, In:*IEEE Symp. Circuits Syst.*,

[17] Saramäki, T. (1993). Finite impulse response filter design, In:*Handbook for Digital Signal Processing*, Mitra, S. K. & Kaiser, J. F. , (Eds.), Ch. 4, (1993) New York, NY:John Wiley and

[18] Saramäki, T. (1994). Generalizations of classical recursive digital filters and their design with the aid of a Remez-type Algorithm, In:*IEEE Symp. Circuits Syst.*, Vol. 2, (May 1994),

[19] Saramäki, T. & Renfors, M. (1995). A Remez-type algorithm for designing digital filters composed of all-pass sections based on phase approximations, In:*Proc. 38th Midwest*

[20] Temes, G. C. & Calahan, D. A. (1967). Computer-aided network optimization the state-of-the-art, In:*Proc. IEEE*, Vol. 55, (Nov. 1967), pp. 1832-1863. Nov. 1967. [21] The MathWorks Inc. (2009). Filter design toolbox user's guide, In:*MATLAB Product Help*,

*Symp. Circuits Syst.*, Vol. 1, (Aug. 1995), Rio de Janiero Brazil, pp. 571-575.

Version 4.6, (Sep. 2009), The MathWorks Inc., Natick, MA.

In: *Compt. Rend. Acad. Sci*, Vol. 199, (1934) pp. 337-340.

Dissertation, Publ. 12, Tampere, Finland, 1981,226 pages.

Vol. 5, (May 1992) San Diego CA, pp. 2577-2580.

Addison-Wesley Pub (Sd).

1987) pp. 32-36.

Sons, pp. 155-277.

London UK, pp. 549-552.

**Chapter 4** 

© 2012 Menghal and Laxmi, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Menghal and Laxmi, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

opportunity helps the students of all courses to realise the limitations of hardware.

**A Virtual Laboratory: Teaching and Research** 

A virtual laboratory for Automatic Control Engineering can provide easy access to university students with regard to engineering applications at anytime and from any computing environment. This interactive learning environment, consist of simulations, demonstrations and exercises, which can fulfill the role of a bridge from passive learning to active engagement and accordingly stimulate deeper thinking; grounding a problem basedlearning environment. The applications are also very important for relating theory to practice, so that the students develop engineering judgment and understand how process behavior can be captured using mathematical models. The undergraduate control engineering at engineering colleges is based on a strong "hands-on" laboratory experience. Regardless of how many fine lectures are given or how many homework problems assigned, the students do not see how control systems work in the real world until they get into the laboratory.They do not understand that they can modify the performance of a physical system to meet design specifications. Only after they complete the laboratory course do they understand the power that they have to become "control gods." After completing the lab experiments on virtual laboratory, the students first characterize the performance of a second order system (a dc servomechanism ES-130).Virtual Laboratory (VLab) has been developed to control Engineering by using MATLAB/SIMULINK. This chapter will also emphasize on the use of Mathematical Modeling and simulation of Feed back Servo trainer (33-100) and study their behavior by using the MATLAB/SIMULINK models and Graphical User Interface (GUI). A graphical user interface is developed which is user friendly and does not require the knowledge of MATLAB. This user can change the parameters of the systems as per his choice or required condition, this computational tool as a part of laboratory experiments will enhance laboratory experience by providing students with the opportunity to compare the practical results with those obtained by computer simulation. Such an

**Tool in Control Engineering Education** 

Prashant M. Menghal and A. Jaya Laxmi

http://dx.doi.org/10.5772/46468

**1. Introduction** 

Additional information is available at the end of the chapter

Prashant M. Menghal and A. Jaya Laxmi

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46468

## **1. Introduction**

A virtual laboratory for Automatic Control Engineering can provide easy access to university students with regard to engineering applications at anytime and from any computing environment. This interactive learning environment, consist of simulations, demonstrations and exercises, which can fulfill the role of a bridge from passive learning to active engagement and accordingly stimulate deeper thinking; grounding a problem basedlearning environment. The applications are also very important for relating theory to practice, so that the students develop engineering judgment and understand how process behavior can be captured using mathematical models. The undergraduate control engineering at engineering colleges is based on a strong "hands-on" laboratory experience. Regardless of how many fine lectures are given or how many homework problems assigned, the students do not see how control systems work in the real world until they get into the laboratory.They do not understand that they can modify the performance of a physical system to meet design specifications. Only after they complete the laboratory course do they understand the power that they have to become "control gods." After completing the lab experiments on virtual laboratory, the students first characterize the performance of a second order system (a dc servomechanism ES-130).Virtual Laboratory (VLab) has been developed to control Engineering by using MATLAB/SIMULINK. This chapter will also emphasize on the use of Mathematical Modeling and simulation of Feed back Servo trainer (33-100) and study their behavior by using the MATLAB/SIMULINK models and Graphical User Interface (GUI). A graphical user interface is developed which is user friendly and does not require the knowledge of MATLAB. This user can change the parameters of the systems as per his choice or required condition, this computational tool as a part of laboratory experiments will enhance laboratory experience by providing students with the opportunity to compare the practical results with those obtained by computer simulation. Such an opportunity helps the students of all courses to realise the limitations of hardware.

© 2012 Menghal and Laxmi, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Menghal and Laxmi, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

With the growing popularity and possibilities of the Internet, web-based teaching is becoming more and more popular in education. The new trend focuses on developing more effective and efficient teaching methods for large groups of students by using interactive web based material. Control systems curricula are often viewed by students as theoretical and highly mathematical. Students are often unable to relate theory to applications in the real world. The obvious solution to this problem, is to include the virtual lab experiment in the control curriculum. Simulation tools are frequently used as an educational aid in automatic control courses. Initially, analog computers with electronic circuits were used for simulating different types of physical processes. When digital computers were present several simulation packages based on numerical techniques were obatinable. Personal computers with low price are universally acknowledged. MATLAB/SIMULINK is a Windows based engineering and science toolbox, which offers valuable interactive demonstrations or the possibility to easily create different simulations related to the theory. It is an integrated technical computing environment that combines numeric computation, advanced graphics and visualization, through powerful Graphical User Interface (GUI). The undergraduate(UG) level engineering lecture class are co- requisites with the lab experiments selected to support and complement the lectures. Most of the universities cover the following lab experiments at UG level:

A Virtual Laboratory: Teaching and Research Tool in Control Engineering Education 71

It is usually recognized that, by studying and experimenting the simulation model a student can experiment with systems that are impossible. He can also work in a laboratory which are being potentially dangerous or huge dimension or very a expensive processes. By using a virtual laboratory, it is possible to obtain some knowledge about its model, operation or

Modeling a system and writing the simulation program contributes to better understanding of its physical principles and properties. By defining the equations of the different parts of the system and how they interact with each other, a student can obtain the clear understanding about the system structure and the way it operates. The feedback mechanical unit 33-100 is a electromechanical unit which comprises of a dc motor, analog tachogenerator, analog input and output potentiometers, absolute and incremental digital encoders and magnetic break as shown in fig.1. The main component of feedback mechanical unit is a DC motor. The mathematical model of DC motor in armature control

Φ = Kf if

mode has been carried out by writing differential equations which are as follows:

**Figure 2.** Feedback Analog Unit 33-002

stability of the system without using a actual hardware.

**2. Mathematical modelling and simulation** 

The air gap flux is proportional to field current,

Where Kf is a constant.


The above listed experiments are performed with the help of feedback mechanical unit 33- 100 and feedback analog unit 33-002 which are as shown in fig.1 and Fig.2.

**Figure 1.** Feed back Mechanical Unit 33-100

**Figure 2.** Feedback Analog Unit 33-002

the following lab experiments at UG level:

**Figure 1.** Feed back Mechanical Unit 33-100

i. Modeling mechanical, electrical and electromechanical systems.

ii. Transfer function, Block diagram and Manson's rules. iii. DC Servo Mechanism: Open loop and Closed loop systems.

vi. Frequency response analysis: Bode and Nyquist plots etc.

100 and feedback analog unit 33-002 which are as shown in fig.1 and Fig.2.

The above listed experiments are performed with the help of feedback mechanical unit 33-

iv. Steady state and transient time response analysis. v. Root locus and introduction to root locus design.

With the growing popularity and possibilities of the Internet, web-based teaching is becoming more and more popular in education. The new trend focuses on developing more effective and efficient teaching methods for large groups of students by using interactive web based material. Control systems curricula are often viewed by students as theoretical and highly mathematical. Students are often unable to relate theory to applications in the real world. The obvious solution to this problem, is to include the virtual lab experiment in the control curriculum. Simulation tools are frequently used as an educational aid in automatic control courses. Initially, analog computers with electronic circuits were used for simulating different types of physical processes. When digital computers were present several simulation packages based on numerical techniques were obatinable. Personal computers with low price are universally acknowledged. MATLAB/SIMULINK is a Windows based engineering and science toolbox, which offers valuable interactive demonstrations or the possibility to easily create different simulations related to the theory. It is an integrated technical computing environment that combines numeric computation, advanced graphics and visualization, through powerful Graphical User Interface (GUI). The undergraduate(UG) level engineering lecture class are co- requisites with the lab experiments selected to support and complement the lectures. Most of the universities cover

> It is usually recognized that, by studying and experimenting the simulation model a student can experiment with systems that are impossible. He can also work in a laboratory which are being potentially dangerous or huge dimension or very a expensive processes. By using a virtual laboratory, it is possible to obtain some knowledge about its model, operation or stability of the system without using a actual hardware.

## **2. Mathematical modelling and simulation**

Modeling a system and writing the simulation program contributes to better understanding of its physical principles and properties. By defining the equations of the different parts of the system and how they interact with each other, a student can obtain the clear understanding about the system structure and the way it operates. The feedback mechanical unit 33-100 is a electromechanical unit which comprises of a dc motor, analog tachogenerator, analog input and output potentiometers, absolute and incremental digital encoders and magnetic break as shown in fig.1. The main component of feedback mechanical unit is a DC motor. The mathematical model of DC motor in armature control mode has been carried out by writing differential equations which are as follows:

The air gap flux is proportional to field current,

Φ = Kf if

Where Kf is a constant.

**Figure 4.** Simulink Model of DC Motor Model

**Figure 5.** Front End User Interface of Virtual Laboratory

The aim of this experiment is to acquaint the user with Open Loop control system characteristics. The circuit diagram for performing this experiment on Analog Servo trainer

The simulation model created with the help of SIMULINK is shown in fig. 7. The response

**2.1. Open loop characteristics** 

– Analog unit 33-110 is as shown in fig. 6.

of the open loop system is shown in fig. 8.

**Figure 3.** Armature controlled DC Motor

The torque Tm developed by the motor is proportional to the product of the armature current and air gap flux,

$$\mathbf{T}\_{\mathrm{m}} = \mathbf{K}\_{1} \mathbf{K} \epsilon \mathbf{i} \epsilon \mathbf{i} \mathbf{a}$$

Where K1 is constant. In the armature controlled DC motor field current is kept constant, so the above equation can be written as

$$\mathbf{T}\_{\mathbf{M}} = \mathbf{K} \tau \mathbf{i}\_{\mathbf{a}}$$

Where KT is known as motor torque constant. The motor back emf being proportional to speed is given as

$$\mathbf{e} \mathbf{e} = \mathbf{K} \mathbb{B} \mathbf{d} \theta / \mathbf{d} \mathbf{t}$$

where Kb equals to back emf constant.

The differential equation of the armature circuit is

$$\mathbf{L\_{a}d/dt \ (i\_{b}) + R\_{a}i\_{a} + e\_{b} = e\_{a}}$$

The torque equation is

$$\mathbf{J}\mathbf{d}^2/\mathbf{d}\mathbf{t}^2\left(\Theta\right) + \mathbf{f}\_\circ \mathbf{d}\Theta/\mathbf{d}\mathbf{t} = \mathbf{T}\mathbf{M} = \mathbf{K}\tau\mathbf{\dot{u}}$$

Thus, from above equations, the transfer function of DC motor is as follows:

$$\mathbf{T(s) = K\pi/\left[s/(R\_A + s\ L\_A)(Js + f\_\circ) + K\tau K\nu\right]}\tag{1}$$

The simulink model of DC motor is derived from the equation 1 is shown in fig. 4. It is the universal model for performing the control system engineering practical in virtual laboratory.

The front end user interface is created to perform the control system practicals with the help of GUI platform of the MATLAB as shown in fig. 5.

**Figure 4.** Simulink Model of DC Motor Model

72 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

The torque Tm developed by the motor is proportional to the product of the armature current

Tm = K1 Kf if ia Where K1 is constant. In the armature controlled DC motor field current is kept constant, so

TM = KT ia Where KT is known as motor torque constant. The motor back emf being proportional to

eb = Kb dθ/dt

La d/dt (ia) + Ra ia + eb = ea

Jd2/dt2 (θ) + fo dθ/dt = TM = KT ia

The simulink model of DC motor is derived from the equation 1 is shown in fig. 4. It is the universal model for performing the control system engineering practical in virtual

The front end user interface is created to perform the control system practicals with the help

T(s) = KT**/** [s{(R a + s La )(Js + f o ) + KT Kb}] (1)

Thus, from above equations, the transfer function of DC motor is as follows:

**Figure 3.** Armature controlled DC Motor

the above equation can be written as

where Kb equals to back emf constant.

The differential equation of the armature circuit is

of GUI platform of the MATLAB as shown in fig. 5.

and air gap flux,

speed is given as

The torque equation is

laboratory.


**Figure 5.** Front End User Interface of Virtual Laboratory

### **2.1. Open loop characteristics**

The aim of this experiment is to acquaint the user with Open Loop control system characteristics. The circuit diagram for performing this experiment on Analog Servo trainer – Analog unit 33-110 is as shown in fig. 6.

The simulation model created with the help of SIMULINK is shown in fig. 7. The response of the open loop system is shown in fig. 8.

**Figure 8.** Output Response of Open Loop Control System

Laboratory performs the following experiments:

c. Steady State Characteristics of DC Motor.

e. Position Control with Controller (P, PI, PD, PID).

f. Speed Control of DC Motor with Controller (P, PI, PD and PID).

a. Speed Vs Input voltage. (b) Step Response. (c) Root locus. (d) Bode plots.

d. Transient Response of DC Motor.

b. Nyquist plots. (f) Pole Zero Map.

a. Open Loop Control System. b. Closed Loop Control System

practicals.

by modeling and simulating of each experiment with the help of MATLAB/SIMULINK. A complete Virtual Laboratory has been developed and successfully implemented. The Virtual

An interactive user interface has been developed using GUI feature of MATLAB to ease and help the user in better understanding and performing of above mentioned practicals. The following are the graphs generated to show the system response in the above mentioned

**Figure 6.** Open Loop system

**Figure 7.** Simulink Model of Open Loop Control System

**Figure 8.** Output Response of Open Loop Control System

by modeling and simulating of each experiment with the help of MATLAB/SIMULINK. A complete Virtual Laboratory has been developed and successfully implemented. The Virtual Laboratory performs the following experiments:

a. Open Loop Control System.

74 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**O/P**

**DC Input** 

**Figure 6.** Open Loop system

**Figure 7.** Simulink Model of Open Loop Control System


An interactive user interface has been developed using GUI feature of MATLAB to ease and help the user in better understanding and performing of above mentioned practicals. The following are the graphs generated to show the system response in the above mentioned practicals.


A Virtual Laboratory: Teaching and Research Tool in Control Engineering Education 77

Every practical system takes finite time to reach its steady state and during this period , it oscillates or increases exponentially.The behaviour of the system gets decided by the type and location of closed loop poles in s–plane. The closed loop poles are dependent on selection of the parmeters of the system. Every system has a tendency to oppose the oscillatory behaviour of the system which is called as **damping.** Now this tendency controls the type of closed loop poles and hence the nature of the response. This damping is measured by a factor called damping ratio of the system. Damping ratio indicates how much dominant the opposition from the system is to the oscillations in the output. In some systems it will be low where system will oscillate but slowly i.e. with damped frequency. If damping ratio is high, system output will not oscillate at all and not only it will be exponential, but also so slow that it will take a very long time to reach a steady state. That is why all practical systems are designed for the damping ratio less than 1 i.e. underdamped.

A Proportional, Integral, Derivative controller (PID controller) is a controller which attempts to correct the error between a measured process variable and a desired set point by calculating and then outputing a corrective action that can adjust the process accordingly. The PID controller calculation (algorithm) involves three separate parameters; the Proportional, the Integral and the Derivative values. The Proportional value determines the reaction to the current error, the Integral determines the reaction based on the sum of recent errors and the Derivative determines the reaction to the rate at which the error has been changing. The weighted sum of these three actions is used to adjust the process via a control element. By "tuning" the three constants in the PID controller algorithm, the PID can provide control action designed for specific process requirements. The response of the controller can be described in terms of the responsiveness of the controller to an error, the degree to which the controller overshoots the set point and the degree of system oscillation. The use of the PID algorithm for control does not guarantee optimal control of the system or system stability. Some applications

**3. Need of controller (P, PD, PI & PID)** 

**Figure 10.** PID Controller

## **2.2. Simulation results**

The simulation results which are obtained by Virtual Laboratory are as follows:

a. Closed Loop Control System

**Figure 9.** (a) Closed loop with Negative Feedback; (b) Closed loop with Positive Feedback

## **3. Need of controller (P, PD, PI & PID)**

76 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

The simulation results which are obtained by Virtual Laboratory are as follows:

(a)

**Figure 9.** (a) Closed loop with Negative Feedback; (b) Closed loop with Positive Feedback

(b)

**2.2. Simulation results** 

a. Closed Loop Control System

Every practical system takes finite time to reach its steady state and during this period , it oscillates or increases exponentially.The behaviour of the system gets decided by the type and location of closed loop poles in s–plane. The closed loop poles are dependent on selection of the parmeters of the system. Every system has a tendency to oppose the oscillatory behaviour of the system which is called as **damping.** Now this tendency controls the type of closed loop poles and hence the nature of the response. This damping is measured by a factor called damping ratio of the system. Damping ratio indicates how much dominant the opposition from the system is to the oscillations in the output. In some systems it will be low where system will oscillate but slowly i.e. with damped frequency. If damping ratio is high, system output will not oscillate at all and not only it will be exponential, but also so slow that it will take a very long time to reach a steady state. That is why all practical systems are designed for the damping ratio less than 1 i.e. underdamped.

**Figure 10.** PID Controller

A Proportional, Integral, Derivative controller (PID controller) is a controller which attempts to correct the error between a measured process variable and a desired set point by calculating and then outputing a corrective action that can adjust the process accordingly. The PID controller calculation (algorithm) involves three separate parameters; the Proportional, the Integral and the Derivative values. The Proportional value determines the reaction to the current error, the Integral determines the reaction based on the sum of recent errors and the Derivative determines the reaction to the rate at which the error has been changing. The weighted sum of these three actions is used to adjust the process via a control element. By "tuning" the three constants in the PID controller algorithm, the PID can provide control action designed for specific process requirements. The response of the controller can be described in terms of the responsiveness of the controller to an error, the degree to which the controller overshoots the set point and the degree of system oscillation. The use of the PID algorithm for control does not guarantee optimal control of the system or system stability. Some applications may necessitate using only one or two modes to provide the appropriate system control. This is achieved by setting the gain of undesired control outputs to zero. A PID controller will be called a PI, PD, P or I controller in the absence of the respective control actions.

A proportional controller (Kp) will have the effect of reducing the rise time and will reduce, but never eliminate, the steady-state error. An integral control (Ki) will have the effect of eliminating the steady-state error, but it may make the transient response worse. A derivative control (Kd) will have the effect of increasing the stability of the system, reducing the overshoot, and improving the transient response. Effects of each of controllers Kp, Kd, and Ki on a closed-loop system are summarized in the table shown below:


**Table 1.** Effects of Controller of Time Response.

The transfer function for PID controller is

$$\text{Transfer Function} = \text{K}\_{\text{P}} + \text{K}\_{\text{d}} \text{S} + \text{K}\_{\text{d}} \text{\textasciic} \tag{2}$$

A Virtual Laboratory: Teaching and Research Tool in Control Engineering Education 79

**Figure 11.** Steady State Characteristics of DC Motor

**Figure 12.** Transient response of DC Motor

### **3.1. Simulation results & analysis**

The simulation results which are obtained by virtual laboratory for DC motors with controller are as follows.

The simulation results of transient behavior of DC motor are shown from Fig. 12 to Fig. 33. The following inferences are drawn from the simulation results of DC Servomechanism:


**Figure 11.** Steady State Characteristics of DC Motor

78 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

called a PI, PD, P or I controller in the absence of the respective control actions.

and Ki on a closed-loop system are summarized in the table shown below:

**Table 1.** Effects of Controller of Time Response.

The transfer function for PID controller is

**3.1. Simulation results & analysis** 

steady state value faster as shown in Fig. 12

shown in Fig. 14, Fig. 15 and Fig. 16.

gain of Derivative and Integral Controller.

controller are as follows.

may necessitate using only one or two modes to provide the appropriate system control. This is achieved by setting the gain of undesired control outputs to zero. A PID controller will be

A proportional controller (Kp) will have the effect of reducing the rise time and will reduce, but never eliminate, the steady-state error. An integral control (Ki) will have the effect of eliminating the steady-state error, but it may make the transient response worse. A derivative control (Kd) will have the effect of increasing the stability of the system, reducing the overshoot, and improving the transient response. Effects of each of controllers Kp, Kd,

**CONTROLLER RISE TIME OVERSHOOT SETTLING TIME STEADY STATE** 

Kp Decrease Increase Small Change Decrease Ki Decrease Increase Increase Eliminate Kd Small Change Decrease Decrease Small Change

The simulation results which are obtained by virtual laboratory for DC motors with

The simulation results of transient behavior of DC motor are shown from Fig. 12 to Fig. 33. The following inferences are drawn from the simulation results of DC Servomechanism:

i. From the output response of the DC motor, it is observed that for every increase in input voltage the output response (i.e. Current and Speed) also increases and reaches its

ii. When proportional controller is used the response is oscillatory and underdamped. It takes more settling time to reach the steady state value as shown in Fig. 13. By using PD controller, the peak overshoot is reduced and thus transient response is improved as shown in Fig.14. System without controller is type I system. The integral controller of the system increases from I to II. As type of system is increased, steady state errors are reduced as shown in Fig. 15. Transient response as well as steady state response can be improved by adjusting value of derivative gain and integral gain as shown in Fig. 16. iii. By using PID controller the transient as well as steady state response improves as

iv. The comparison of the transient response of DC motor is shown in Fig. 33 with P, PI, PD and PID controller. As the Integral gain is increased the system response becomes more sluggish, however the tuning of the controller can be judged by deciding on the proper

Transfer Function = Kp + KdS+ Ki/ (2)

**ERROR** 

**Figure 12.** Transient response of DC Motor

**Figure 15.** Transient Response of the DC motor with PI Controller

**Figure 16.** Transient Response of the DC motor with PID Controller

**Figure 13.** Transient Response of the DC motor without Controller.

**Figure 14.** Transient Response of the DC motor with PD Controller

**Figure 15.** Transient Response of the DC motor with PI Controller

80 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 13.** Transient Response of the DC motor without Controller.

**Figure 14.** Transient Response of the DC motor with PD Controller

**Figure 16.** Transient Response of the DC motor with PID Controller

82 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 20.** Root Locus of the DC motor with PID Controller

**Figure 21.** Bode plot of the DC motor without Controller

**Figure 22.** Bode plot of the DC motor with PD Controller

**Figure 17.** Root Locus of the DC motor without Controller

**Figure 18.** Root Locus of the DC motor with PD Controller

**Figure 19.** Root Locus of the DC motor with PI Controller

**Figure 20.** Root Locus of the DC motor with PID Controller

**Figure 17.** Root Locus of the DC motor without Controller

**Figure 18.** Root Locus of the DC motor with PD Controller

**Figure 19.** Root Locus of the DC motor with PI Controller

**Figure 21.** Bode plot of the DC motor without Controller

**Figure 22.** Bode plot of the DC motor with PD Controller

**Figure 26.** Nyquist plot of the DC motor with PD Controller

**Figure 27.** Nyquist plot of the DC motor with PI Controller

**Figure 28.** Nyquist plot of the DC motor with PID Controller

**Figure 23.** Bode plot of the DC motor with PI Controller

**Figure 24.** Bode plot of the DC motor with PID Controller

**Figure 25.** Nyquist plot of the DC motor without Controller

**Figure 26.** Nyquist plot of the DC motor with PD Controller

**Figure 23.** Bode plot of the DC motor with PI Controller

**Figure 24.** Bode plot of the DC motor with PID Controller

**Figure 25.** Nyquist plot of the DC motor without Controller

**Figure 27.** Nyquist plot of the DC motor with PI Controller

**Figure 28.** Nyquist plot of the DC motor with PID Controller

Nyquist Diagram

P PI PD PID

Frequency (rad/sec)

100 102 104

Bode Diagram

Real Axis


**Figure 32.** Pole Zero Map of the DC motor with PID Controller

Step Response

Pole-Zero Map

Time (sec)

0 2 4 6 8 10 12

Imaginary Axis



50 100

0

0.5

1

1.5

Amplitude

Real Axis


**Figure 33.** Comparison of the Transient Response of DC Motor with P, PI, PD, PID Controller

x 105

**4. Integration of Virtual laboratory into Control Engineering courses** 

The Virtual laboratory for control system has been developed using MATLAB/SIMULINK and integrated into the Control Engineering Curriculum of UG courses. It consists of the complete details of mathematical modeling using the modules available in SIMULINK for open loop, closed loop, PID control, Transient, frequency response of DC Motor, Position and speed control of DC motor using controller experiments. Then the programming of each

Imaginary Axis



1 2



0 100

Phase (deg)

Magnitude (dB)

**Figure 29.** Pole Zero Map of the DC motor without Controller

**Figure 30.** Pole Zero Map of the DC motor with PD Controller

**Figure 31.** Pole Zero Map of the DC motor with PI Controller

**Figure 32.** Pole Zero Map of the DC motor with PID Controller

**Figure 29.** Pole Zero Map of the DC motor without Controller

**Figure 30.** Pole Zero Map of the DC motor with PD Controller

**Figure 31.** Pole Zero Map of the DC motor with PI Controller

**Figure 33.** Comparison of the Transient Response of DC Motor with P, PI, PD, PID Controller

#### **4. Integration of Virtual laboratory into Control Engineering courses**

The Virtual laboratory for control system has been developed using MATLAB/SIMULINK and integrated into the Control Engineering Curriculum of UG courses. It consists of the complete details of mathematical modeling using the modules available in SIMULINK for open loop, closed loop, PID control, Transient, frequency response of DC Motor, Position and speed control of DC motor using controller experiments. Then the programming of each

selected module is carried out. For the ease of an user, a graphical user interface is modeled for all the above said experiments on a single platform as shown in Fig. 5. The platform includes

A Virtual Laboratory: Teaching and Research Tool in Control Engineering Education 89

graphics, to illustrate transient and steady-state performance and stability analysis of control system under various parameter controls. The user can change the parameters of the systems as per his choice or required condition. Such an opportunity helps students of all courses realise the limitations of hardware. A Virtual laboratory for Automatic Control (AC) allows students an easy access to different applications, simulations related to the theory they studied. These interactive demos present in a tutorial manner the influence of the different parameters of the mathematical model to the system behavior. These simulations provide a more intuitive and more practical approach for the abstract theory of Automatic Control. The advantage of the approach presented here is the use of the available simulation tools. The user can focus on the learning and understanding of problems and concepts, as he/she doesn't have to master the MATLAB programming environment. The scope of this Vlab is extending to the remote control laboratory using MATLAB /SIMULINK applications to control system engineering. The MATLAB/SIMULINK models and GUI representation of DC Motor using controller will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, with the computer-generated graphics, to illustrate transient and steady-state performance and stability analysis of DC motor under various parameter controls. The user can change the parameters of the system as per his/her choice or required condition. Thus, this computational tool, as a part of laboratory experiments will enhance practical experience by providing students with the opportunity to compare the results of laboratory experiments with those obtained by computer simulation. Such an

opportunity helps students of all courses realize the limitations of hardware.

*Radar & Control System Dept, Faculty of Electronics, Military College of Electronics and* 

*Electrical & Electronics Engg. Dept. Jawaharlal Nehru Technological University, Hyderabad College* 

Eric Hoffer, in **"Reflections On the Human Condition"** quotes "the hardest arithmetic to master is to enable us to count our blessings." In pursuit of accomplishing a goal, there is incessant need for constant stimulation and inspiration to persevere and attain. There are also times when obscurity threatens to conceal the desire to succeed with the drape of uncertainnities and hindrances and it is in those hours of trepidation that the Guru rekindles the spark within us with flames of guidance and mentoring. This is a humble effort on the part of me to undertake the enormous responsibility of expressing in words the emotions and gratitude felt towards all our gurus, without whose ardor and continuous assurances, this voyage of intense erudition would not have been possible. I would like to thank Head of Department (Radar & Control System), V. K. Pokhriyal, Dean Faculty of Electronics V S

*Mechanical Engineering, Secunderabad, Andhra Pradesh, India* 

*of Engineering, Kukatpally, Hyderabad, Andhra Pradesh, India* 

**Author details** 

A Jaya Laxmi

Prashant M. Menghal

**Acknowledgement** 


## **4.1. The educational use of the model**

The mathematical model of DC Motor is very useful for carrying out transient and steady state analysis and understanding the basic concepts of control systems with and without the controller. It is a fundamental system useful to any course of instruction including basic subjects in Engineering at undergraduate courses. The MATLAB/SIMULINK model and GUI representation of DC motor using controller will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, with the computer-generated graphics, to illustrate transient and steady-state performance and stability analysis of DC motor under various parameter controls. The user can change the parameters of the system as per his/her choice or required condition. Thus this computational tool as a part of laboratory experiments will enhance laboratory experience by providing students with the opportunity to compare the results of laboratory experiments with those obtained by computer simulation. Such an opportunity helps students of all courses realise the limitations of hardware.

## **5. Conclusion**

The MATLAB/SIMULINK models and GUI representation of closed loop, open loop, PID Control, etc. of Feedback Control System 33-001 will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, through the computer-generated graphics, to illustrate transient and steady-state performance and stability analysis of control system under various parameter controls. The user can change the parameters of the systems as per his choice or required condition. Such an opportunity helps students of all courses realise the limitations of hardware. A Virtual laboratory for Automatic Control (AC) allows students an easy access to different applications, simulations related to the theory they studied. These interactive demos present in a tutorial manner the influence of the different parameters of the mathematical model to the system behavior. These simulations provide a more intuitive and more practical approach for the abstract theory of Automatic Control. The advantage of the approach presented here is the use of the available simulation tools. The user can focus on the learning and understanding of problems and concepts, as he/she doesn't have to master the MATLAB programming environment. The scope of this Vlab is extending to the remote control laboratory using MATLAB /SIMULINK applications to control system engineering. The MATLAB/SIMULINK models and GUI representation of DC Motor using controller will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, with the computer-generated graphics, to illustrate transient and steady-state performance and stability analysis of DC motor under various parameter controls. The user can change the parameters of the system as per his/her choice or required condition. Thus, this computational tool, as a part of laboratory experiments will enhance practical experience by providing students with the opportunity to compare the results of laboratory experiments with those obtained by computer simulation. Such an opportunity helps students of all courses realize the limitations of hardware.

## **Author details**

88 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

includes

a. Buttons to perform all the experiments. b. Controls to change the system parameters.

separate window.

limitations of hardware.

**5. Conclusion** 

c. View the results of the experiments in the graphical form. d. Compare system response such time and frequency domain .

the experiments performed can be viewed.

g. SIMULINK diagram of the current experiment.

**4.1. The educational use of the model** 

h. An EXIT button to come out of the GUI page of project.

j. A button to see the SIMULINK diagram of each experiments.

i. A refresh button to clear all the graphs and diagram displayed earlier.

subjects in electronics and mechanical at undergraduate level.

k. Menu bar to see different values in the graph, to zoom in/out and other controls.

selected module is carried out. For the ease of an user, a graphical user interface is modeled for all the above said experiments on a single platform as shown in Fig. 5. The platform

e. A menu bar on activation of which brief description along with the transfer function of

f. A menu bar to simulate directly the SIMULINK circuits diagrams and see the results in

l. A course ware for each experiment is prepared and successfully integrated into the control engineering laboratory. Vlab is very useful for performing experiments and understanding the basic concepts of feedback systems in the control system laboratory. It is a fundamental system which is useful to any course of instruction including basic

The mathematical model of DC Motor is very useful for carrying out transient and steady state analysis and understanding the basic concepts of control systems with and without the controller. It is a fundamental system useful to any course of instruction including basic subjects in Engineering at undergraduate courses. The MATLAB/SIMULINK model and GUI representation of DC motor using controller will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, with the computer-generated graphics, to illustrate transient and steady-state performance and stability analysis of DC motor under various parameter controls. The user can change the parameters of the system as per his/her choice or required condition. Thus this computational tool as a part of laboratory experiments will enhance laboratory experience by providing students with the opportunity to compare the results of laboratory experiments with those obtained by computer simulation. Such an opportunity helps students of all courses realise the

The MATLAB/SIMULINK models and GUI representation of closed loop, open loop, PID Control, etc. of Feedback Control System 33-001 will definitely work as a teaching tool and support the classroom teaching by enabling the faculty, through the computer-generated Prashant M. Menghal *Radar & Control System Dept, Faculty of Electronics, Military College of Electronics and Mechanical Engineering, Secunderabad, Andhra Pradesh, India* 

A Jaya Laxmi *Electrical & Electronics Engg. Dept. Jawaharlal Nehru Technological University, Hyderabad College of Engineering, Kukatpally, Hyderabad, Andhra Pradesh, India* 

## **Acknowledgement**

Eric Hoffer, in **"Reflections On the Human Condition"** quotes "the hardest arithmetic to master is to enable us to count our blessings." In pursuit of accomplishing a goal, there is incessant need for constant stimulation and inspiration to persevere and attain. There are also times when obscurity threatens to conceal the desire to succeed with the drape of uncertainnities and hindrances and it is in those hours of trepidation that the Guru rekindles the spark within us with flames of guidance and mentoring. This is a humble effort on the part of me to undertake the enormous responsibility of expressing in words the emotions and gratitude felt towards all our gurus, without whose ardor and continuous assurances, this voyage of intense erudition would not have been possible. I would like to thank Head of Department (Radar & Control System), V. K. Pokhriyal, Dean Faculty of Electronics V S

Randhwa,Head of Institution SM Mehta SM,VSM\*\*,who unlocked for me the opportunities and the resources to explore potentials beyond my envision through his inexorable confidence in my capabilities. I am obliged to the unrelenting espousal and conviction of my Ph.D. supervisor, Dr A. Jaya Laxmi, who was a catalyst in leading towards the completion of this chapter.

**Chapter 5** 

© 2012 Pekař et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Pekař et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**On Finite-Dimensional Transformations of** 

Linear time-invariant time delay systems (LTI TDS) have usually been assumed to contain delay elements in input-output relations only. All the system dynamics has been hence modeled by point accumulations in the form of a set of ordinary differential equations. The Laplace transform then results in a transfer function expressed by a serial combination of a delayless term and a delay element. However, this conception is somewhat restrictive in effort to fit the real plant dynamics because inner feedbacks are often of the time-distributed

*Anisochronic* (or *hereditary*) TDS models, on the other hand, offer a more universal dynamics description applying both integrators and delay elements either in lumped or distributed form so that delays appear on the left side of a differential equation which is no longer *ordinary* (ODE) but rather *functional* (FDE) - this brings the concept of *internal* (or *state*) delays. In contrast to undelayed systems, the main difference in dynamics is that their spectra are infinite in general. In the further text, an abbreviation TDS means LTI TDS

Already in (Volterra, 1928) differential equations incorporating the past states when studying predator-pray models were formulated. The theory of these models has been then developed by many outstanding authors, see e.g. (Bellman & Cooke, 1963), (Krasovskii, 1963), (Kolmanovskii & Nosov, 1986), (Zítek, 1983), (Górecki et al., 1989), and especially (Hale & Verduyn Lunel, 1993) and (Nicolescu, 2001), to name a few. Aftereffect phenomenon is included in many processes, e.g. in chemical processes (Zítek & Hlava, 2001), heat exchange networks (Zítek, 1997), in models of mass flow in sugar factory (Findeisen et al., 1970), in metallurgic processes (Morávka & Michálek, 2008), etc. Plenty of

**Anisochronic Controllers Designed by** 

**Algebraic Means: A User Interface** 

Libor Pekař, Eva Kurečková and Roman Prokop

containing state delays with or without input-output delays.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46472

**1. Introduction** 

or delayed nature.

Once again we would like to express our heartfelt gratitude to each and every person who was pivotal in the successful architecting and completion of this chapter and without whom this chapter would not have been a reality.

## **6. References**


www.mathworks.com
