# Design of Digital IIR Filter with Conflicting Objectives Using Hybrid Gravitational Search Algorithm

Authors: D. S. Sidhu; J. S. Dhillon; Dalvir Kaur
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

In the recent years, the digital IIR filter design as a single objective optimization problem using evolutionary algorithms has gained much attention. In this paper, the digital IIR filter design is treated as a multiobjective problem by minimizing the magnitude response error, linear phase response error and optimal order simultaneously along with meeting the stability criterion. Hybrid gravitational search algorithm (HGSA) has been applied to design the digital IIR filter. GSA technique is hybridized with binary successive approximation (BSA) based evolutionary search method for exploring the search space locally. The relative performance of GSA and hybrid GSA has been evaluated by applying these techniques to standard mathematical test functions. The above proposed hybrid search techniques have been applied effectively to solve the multiparameter and multiobjective optimization problem of low-pass (LP), high-pass (HP), band-pass (BP), and band-stop (BS) digital IIR filter design. The obtained results reveal that the proposed technique performs better than other algorithms applied by other researchers for the design of digital IIR filter with conflicting objectives.


---

## Body

## 1. Introduction

The electrical signal is the best method to transfer information to longer distances. As the signal becomes complex, we require to process the signal to extricate the information available in the signal. Generally, signals are classified into two categories, that is, continuous-time and discrete-time signals. Digital signal processing (DSP) is a numerical manipulation of discrete-time signal data. With the advent of digital circuit technology, cheaper, simple, and faster digital computers with huge memory storage capability are able to perform complex real time digital signal processing. In the recent period, digital signal processing (DSP) is the area of interest for engineers and scientists due to their accuracy, reliability, flexibility, low cost, and small physical size. Filtering is an essential part of digital signal processing for different applications such as telecommunications, speech processing, image processing, consumer electronics, biomedical systems, industrial applications, and military electronics. As compared to the analog filters, digital filters are flexible, reliable, and versatile. The sampled values of input signal and the transfer function of digital filter can be easily stored in the memory [1, 2].Digital filters are classified as finite impulse response (FIR) filters and infinite impulse response (IIR) filters. IIR filter provides better response than FIR filters with the same order or same number of coefficients. The stability of the digital IIR filter can be obtained by limiting the parameter space in a suitable range. Multimodal error surface of IIR filter may force gradient based algorithm to get stuck in local minima [3]. In order to avoid the above problem and to find the global minima, global optimization techniques like genetic algorithm (GA), particle swarm optimization (PSO), differential evolution (DE), predator prey optimization (PPO), artificial bee colony (ABC), and so forth have been used by the researchers [3–7] to design digital IIR filters.The synergy between exploration and exploitation is a thrust area in the optimization process. The generalization of heuristics allows the wide applicability but it limits the efficiency of heuristics and simply delivers mediocre performance. The hybridization of general heuristics with problem-specific heuristics provides better results for more complicated problems. The technique which uses general heuristics as global search and a specific heuristic as local search is generally called hybrid algorithm. With a proper tuning, hybrid algorithms with a good explorative ability as a population-based global search algorithm also provide a good exploitive performance as a local search algorithm [8]. Many researchers have applied various hybrid optimization techniques for design of digital IIR and FIR filters. Singh et al. [9] had applied DE hybridized with pattern search technique for the design of digital IIR filter. DE has been used for global search technique and pattern search method as local search technique. Kaur and Dhillon [10] had proposed integrated cat swarm optimization (CSO) and DE for digital IIR filter design. So Kaur and Dhillon discussed hybridization of two global search techniques. Hua et al. [11] have hybridized DE and PSO for designing 2D FIR filter. Vasundhara et al. [12] had proposed hybridization of adaptive DE and PSO for FIR filter design. Both researchers have also considered the hybridization of two different global search techniques. Bindiya and Elias [13] had applied different modified metaheuristic algorithms for the design of multiplier-less nonuniform channel filters and compared their results. They concluded that gravitational search algorithm (GSA) is the fastest and also gives acceptable values of the frequency-response characteristics and complexity. Saha et al. [14] had applied GSA for design and simulation of FIR band-pass and band-stop filters and compared the results with results obtained by DE, PSO, RGA, and BBO. It has been concluded that GSA converges faster than other techniques. Saha et al. [15] had also applied GSA with wavelet mutation for design of IIR filter and concluded that GSAWM provides better results as compared to other techniques. The above researchers considered the design problem as a single objective optimization problem and applied the techniques for minimizing the magnitude response error of digital filter only with supplementary conditions.Yu and Xinjie [16] had applied cooperative coevolutionary genetic algorithm to design the digital IIR filter. Magnitude response and phase response are simultaneously optimized along with the minimum filter order. Wang et al. [17] applied local search operator enhanced multiobjective evolutionary algorithm (LS-MOEA) for design of digital IIR filter with multiple objectives. Wang et al. [17] has also applied nondominated sorted based genetic algorithm II (NSGA-II) for the above problem. Kaur et al. [18] applied real coded genetic algorithm (RCGA) for design of digital IIR filter with multiple objectives. The literature survey reveals that metaheuristic methods generally suffer from either the exploration or exploitation capabilities leading to lack of diversity and stagnation to local minima, respectively. Involvement of many parameters leads to parameter tuning problems and becomes time consuming. Mostly the aforementioned methods lack to maintain balance between exploration and exploitation.Any proposed optimization method may be a good choice to design digital filters that may have exclusive advantages like robust, global search capability, little information requirement, ease of implementation, parallelism, no requirement of continuous, and differentiable objective function. The intent of this paper is to propose the global search technique, gravitational search algorithm (GSA) to design the digital IIR filter. The global search technique is then hybridized with evolutionary search method for exploring the search space locally. With the increase in number of filter coefficients, the number of comparisons increases exponentially in evolutionary search technique during exploration of search space. To reduce such comparisons, a unique binary successive approximation (BSA) strategy has been implemented while exploring the search space locally. The opposition-based learning strategy is applied to have a start with good initial population. In this paper the design of digital IIR filter is treated as a multiobjective problem by simultaneously minimizing the magnitude response error, linear phase response error, and optimal order along with meeting the stability criterion. The proposed hybrid technique is applied to solve the multiparameter and multiobjective optimization problem of low-pass, high-pass, band-pass, and band-stop digital IIR filter design. Owing to the imprecise nature of decision maker judgement, it is presumed that decision maker may have fuzzy or imprecise goals for each objective [19]. Hybrid GSA is used to simulate trade-off between conflicting objectives that is magnitude response error and phase response. Once the trade-off has been obtained, fuzzy set theory helps decision maker to decide the optimal design of IIR digital filter over the noninferior trade-off curve.This paper is divided into five sections. In Section2, digital IIR filter problem is described. The gravitational search algorithm (GSA) and binary successive approximation (BSA) based evolutionary search method for optimal IIR filter design are explained in Section 3. The GSA and hybrid GSA have been applied to standard mathematical test functions in Section 4 to evaluate the relative performance. In Section 5, digital IIR filter has been designed by applying the proposed hybrid heuristic search technique and obtained results are evaluated and compared with results of other authors. The final conclusion and discussions are outlined in Section 6.


## 2. IIR Filter Design Problem

Digital infinite impulse response (IIR) filter design requires the optimization of set of filter coefficients which meet various specification of the filter such as pass-band width and gain, stop-band width and attenuation, edge frequencies of pass band and stop band, peak ripples in pass band and stop band, and linear phase response in pass band and transition band.The cascaded digital IIR filter has several first order and second order sections together and is easy to implement. Puttingz
=
e
j
ω, the structure of cascading type IIR filter can be expressed as [20](1)
H
ω
,
x
=
A
∏
i
=
1
M
1
+
a
i
e
-
j
ω
1
+
b
i
e
-
j
ω
∏
k
=
1
N
1
+
p
k
e
-
j
ω
+
q
k
e
-
2
j
ω
1
+
r
k
e
-
j
ω
+
s
k
e
-
2
j
ω
,where(2)
x
=
A
,
a
1
,
b
1
,
…
,
a
M
,
b
M
,
p
1
,
q
1
,
r
1
,
s
1
,
…
,
p
N
,
q
N
,
r
N
,
s
N
T
.Vectorx denotes the set of filter coefficients of first order and second order sections of dimensions D
×
1 with D
=
2
M
+
4
N
+
1.The coefficients of the transfer functionH
ω
,
x are approximated during the design of filter. The transfer function H
ω
,
x is compared this with the ideal transfer function H
I
ω. The magnitude error function E
1
x is obtained from these two values as given below:(3a)
Minimize
E
1
x
=
∑
i
k
H
d
ω
i
-
H
ω
i
,
x
,where (3b)
H
d
ω
i
=
1
,
for
ω
i
∈
pass band
0
,
for
ω
i
∈
stop band.The ripple magnitude of pass-bandδ
1
(
x
) and of stop-band δ
2
(
x
) are defined as below(4a)
δ
1
x
=
max
ω
i
⁡
H
ω
i
,
x
-
min
ω
i
⁡
H
ω
i
,
x
for
ω
i
∈
pass band
,
(4b)
δ
2
x
=
max
ω
i
⁡
H
ω
i
,
x
for
ω
i
∈
stop band.During the design of filter, the linear phase response is also optimized for both pass band and transition band. The expression for phase is defined as below:(5)
∅
=
tan
-
1
⁡
imaginary part of numerator
real part of numerator
-
tan
-
1
⁡
imaginary part of denominator
real part of denominator
.The first order difference in phase is determined as(6)
Δ
phase
=
Δ
∅
1
,
Δ
∅
2
,
…
,
Δ
∅
n
-
1
,where(7)
Δ
∅
i
=
Δ
∅
i
+
1
-
Δ
∅
i
i
=
1,2
,
…
,
n
-
1
.N is the total number of sampling points in pass band and transition band. The phase response will be linear if all the elements of Δ
phase have the equal values [18]. So the second objective function is to minimize variance of phase differences and is expressed as below:(8)
Minimize
E
2
=
var
⁡
Δ
∅
n
∅
n
∈
pass band
∪
transition band.Consider all these objectives, subject to the following stability constraints:(9a)
1
+
b
i
≥
0
i
=
1,2
,
…
,
M
,
(9b)
1
-
b
i
≥
0
i
=
1,2
,
…
,
M
,
(9c)
1
-
s
k
≥
0
k
=
1,2
,
…
,
N
,
(9d)
1
+
r
k
+
s
k
≥
0
k
=
1,2
,
…
,
N
,
(9e)
1
-
r
k
+
s
k
≥
0
k
=
1,2
,
…
,
N
.For the design of digital IIR filter with the above multiobjective optimization criterion, the fuzzy membership functions are assigned to magnitude response error function and phase response error function along with stability constraints. The fuzzy sets are stated by their membership functions. Such membership functions describe the degree of membership in certain fuzzy sets using values from 0 to 1. The membership function value 0 indicates incompatibility with the sets whereas 1 means full compatibility [19]. The member functions are used to convert minimization to maximization optimization problem and constraints are converted into objective function to be maximized.


### 2.1. Magnitude Response Error

The main objective to design digital IIR filter is to minimize the magnitude response error in a predefined pass-band and stop-band and within prescribed permissible ripples. In the prescribed pass-band the signal is allowed to pass and in the prescribed stop-band the signal is restricted.The membership function of magnitude response in pass-band is described in Figure1.Figure 1
Membership function of magnitude response during pass-band.Mathematical representation of membership function of magnitude response during pass-band is given as below:(10)
u
n
H
ω
p
=
0
H
ω
p
≤
H
ω
p
m
i
n
or
H
ω
p
≥
H
ω
p
m
a
x
H
ω
p
-
H
ω
p
m
i
n
1
-
δ
p
1
-
H
ω
p
m
i
n
H
ω
p
m
i
n
<
H
ω
p
<
1
-
δ
p
1
H
ω
p
m
a
x
-
H
ω
p
H
ω
p
m
a
x
-
1
+
δ
p
2
1
+
δ
p
2
<
H
ω
p
<
H
ω
p
m
a
x
1
1
-
δ
p
1
≤
H
ω
p
≤
1
+
δ
p
2
,where ω
p is the frequency in pass-band; H
(
ω
p
) is the magnitude response in pass-band; H
(
ω
p
)
m
i
n is the minimum magnitude in pass band; H
(
ω
p
)
m
a
x is maximum magnitude in the pass-band; δ
p
1 and δ
p
2 are the maximum tolerable deviations in magnitude in pass-band; δ
p
2 is kept almost nearer to zero; u
n
H
(
ω
p
) is the membership function of magnitude response in pass-band; n
∈
p
a
s
s
b
a
n
d.The membership function of magnitude response in stop-band is described in Figure2.Figure 2
Membership function of magnitude response during stop band.Mathematical expression of membership function of magnitude response during stop-band is given as below:(11)
μ
n
H
ω
s
=
0
H
ω
s
≥
H
ω
s
m
a
x
H
ω
s
m
a
x
-
H
ω
s
H
ω
s
m
a
x
-
δ
s
δ
s
<
H
ω
s
<
H
ω
s
m
a
x
1
H
ω
s
≤
δ
s
,where ω
s is the frequency in stop band; H
(
ω
s
) is the magnitude response in stop-band; H
s
m
a
x is maximum magnitude in the stop-band; δ
s is the maximum tolerable deviations in magnitude in stop-band; u
n
H
(
ω
s
) is the membership function of magnitude response in stop-band; n
∈
stop band.The main objective is to maximize the membership function of magnitude response in both pass-band and stop-band and is represented as below:(12)
μ
1
=
min
⁡
min
⁡
μ
n
H
ω
p
n
∈
pass band
,
min
⁡
μ
n
H
ω
s
n
∈
stop band
.The maximum value ofμ
1 provides the design of digital filter having minimum value of magnitude response error.


### 2.2. Phase Response Error

As described in (8), the phase response error E
2 is to be minimized. This minimization problem is converted into maximization problem as described as below:(13)
μ
2
=
1
1
+
E
2
,where μ
2 is the objective function of phase response. Maximum value of objective function μ
2 provides the minimum value of linear phase response error E
2.


### 2.3. Membership Function of Stability Constraints

The stability constraints for the design of digital IIR filter are obtained by using the Jury method [21] on the filter coefficients given in (9a). Using fuzzy set theory, the membership function of stability constraints given in (9a) is shown in Figure 3.Figure 3
Membership function of stability constraint.Mathematically, membership function of stability constraint given in (9a) is given below:(14)
μ
1
i
=
0
if
1
+
b
i
≤
0
1
if
1
+
b
i
>
0
i
=
1,2
,
…
,
M
.Similarly membership function for stability constraints given in (9b) to (9e) can be described. The overall objective to meet these above five stability constraints is mathematically given below and is to be maximized:(15)
μ
3
=
1
5
∑
k
=
1
2
1
M
∑
i
=
1
M
μ
k
i
+
∑
k
=
3
5
1
N
∑
j
=
1
N
μ
k
j
.If the valueμ
3 is 1 that means all constraints are satisfied. Otherwise it gives the percentage level of satisfaction of constraints.


### 2.4. Multiobjective Problem Formulation

The task of digital IIR filter design is to find an optimum structure with minimum magnitude response error and minimum linear phase response error while satisfying the stability constraints. By applying fuzzy set theory, multiobjective constrained problem is converted into multiobjective unconstrained optimization problem and is stated as below:(16)
Maximize
μ
1
X
,
μ
2
X
,
μ
3
X
T
,where μ
1
(
X
) is the membership function of magnitude response error given in (12); μ
2
(
X
) is the membership function of linear phase response error given in (13); μ
3
(
X
) is the membership function of stability constraints given in (15); X is a vector decision variable of dimensions D
×
1 with D
=
2
M
+
4
N
+
1.The objective is to find the value of filter coefficients being decision variables,X, which maximizes the entire objective functions simultaneously. The value of membership function indicates how much the solution satisfies the μ
i objective on the scale from 0 to 1. The maximum satisfaction of membership function for any filter coefficient combination is obtained by taking the intersection of the membership functions of participating objectives and is expressed as below:(17)
fit
j
t
=
min
⁡
μ
1
j
,
μ
2
j
,
μ
3
j
j
=
1,2
,
…
,
N
P
.


### 2.5. Optimal Order of Digital Filter

The expression for optimal order of digital IIR filter is described below:(18)
O
=
∑
j
=
0
M
α
j
+
2
∑
k
=
0
N
β
k
,where α
j and β
k are j
t
h and k
t
h control genes of corresponding first order and second order blocks, respectively, and the value of genes will be either 1 or 0. M and N are the maximum number of first order and second order blocks, respectively. Maximum order of the digital IIR filter will be M
+
2
N. The order of the filter has been determined by control genes as shown in Figure 4. The coding method has been taken from Yu and Xinjie [16]. The value of control gene determines whether the particular block will be considered for filter design or not. The block will be considered activated when the corresponding control gene is 1. The number of binary bits used to generate control genes depends upon the value of M and N. The decision vector x shown in (2) is modified as below:(19)
x
=
U
,
V
,
A
,
a
1
,
b
1
,
…
,
a
M
,
b
M
,
p
1
,
q
1
,
r
1
,
s
1
,
…
,
p
N
,
q
N
,
r
N
,
s
N
T
.Figure 4
Activation and deactivation of filter coefficients with control genes.The variableU is a positive integer with maximum value of (
2
M
-
1
) and the variable V is a positive integer with maximum value of (
2
N
-
1
). The variable U and V will also be optimized along with the coefficients of the filter.


## 2.1. Magnitude Response Error

The main objective to design digital IIR filter is to minimize the magnitude response error in a predefined pass-band and stop-band and within prescribed permissible ripples. In the prescribed pass-band the signal is allowed to pass and in the prescribed stop-band the signal is restricted.The membership function of magnitude response in pass-band is described in Figure1.Figure 1
Membership function of magnitude response during pass-band.Mathematical representation of membership function of magnitude response during pass-band is given as below:(10)
u
n
H
ω
p
=
0
H
ω
p
≤
H
ω
p
m
i
n
or
H
ω
p
≥
H
ω
p
m
a
x
H
ω
p
-
H
ω
p
m
i
n
1
-
δ
p
1
-
H
ω
p
m
i
n
H
ω
p
m
i
n
<
H
ω
p
<
1
-
δ
p
1
H
ω
p
m
a
x
-
H
ω
p
H
ω
p
m
a
x
-
1
+
δ
p
2
1
+
δ
p
2
<
H
ω
p
<
H
ω
p
m
a
x
1
1
-
δ
p
1
≤
H
ω
p
≤
1
+
δ
p
2
,where ω
p is the frequency in pass-band; H
(
ω
p
) is the magnitude response in pass-band; H
(
ω
p
)
m
i
n is the minimum magnitude in pass band; H
(
ω
p
)
m
a
x is maximum magnitude in the pass-band; δ
p
1 and δ
p
2 are the maximum tolerable deviations in magnitude in pass-band; δ
p
2 is kept almost nearer to zero; u
n
H
(
ω
p
) is the membership function of magnitude response in pass-band; n
∈
p
a
s
s
b
a
n
d.The membership function of magnitude response in stop-band is described in Figure2.Figure 2
Membership function of magnitude response during stop band.Mathematical expression of membership function of magnitude response during stop-band is given as below:(11)
μ
n
H
ω
s
=
0
H
ω
s
≥
H
ω
s
m
a
x
H
ω
s
m
a
x
-
H
ω
s
H
ω
s
m
a
x
-
δ
s
δ
s
<
H
ω
s
<
H
ω
s
m
a
x
1
H
ω
s
≤
δ
s
,where ω
s is the frequency in stop band; H
(
ω
s
) is the magnitude response in stop-band; H
s
m
a
x is maximum magnitude in the stop-band; δ
s is the maximum tolerable deviations in magnitude in stop-band; u
n
H
(
ω
s
) is the membership function of magnitude response in stop-band; n
∈
stop band.The main objective is to maximize the membership function of magnitude response in both pass-band and stop-band and is represented as below:(12)
μ
1
=
min
⁡
min
⁡
μ
n
H
ω
p
n
∈
pass band
,
min
⁡
μ
n
H
ω
s
n
∈
stop band
.The maximum value ofμ
1 provides the design of digital filter having minimum value of magnitude response error.


## 2.2. Phase Response Error

As described in (8), the phase response error E
2 is to be minimized. This minimization problem is converted into maximization problem as described as below:(13)
μ
2
=
1
1
+
E
2
,where μ
2 is the objective function of phase response. Maximum value of objective function μ
2 provides the minimum value of linear phase response error E
2.


## 2.3. Membership Function of Stability Constraints

The stability constraints for the design of digital IIR filter are obtained by using the Jury method [21] on the filter coefficients given in (9a). Using fuzzy set theory, the membership function of stability constraints given in (9a) is shown in Figure 3.Figure 3
Membership function of stability constraint.Mathematically, membership function of stability constraint given in (9a) is given below:(14)
μ
1
i
=
0
if
1
+
b
i
≤
0
1
if
1
+
b
i
>
0
i
=
1,2
,
…
,
M
.Similarly membership function for stability constraints given in (9b) to (9e) can be described. The overall objective to meet these above five stability constraints is mathematically given below and is to be maximized:(15)
μ
3
=
1
5
∑
k
=
1
2
1
M
∑
i
=
1
M
μ
k
i
+
∑
k
=
3
5
1
N
∑
j
=
1
N
μ
k
j
.If the valueμ
3 is 1 that means all constraints are satisfied. Otherwise it gives the percentage level of satisfaction of constraints.


## 2.4. Multiobjective Problem Formulation

The task of digital IIR filter design is to find an optimum structure with minimum magnitude response error and minimum linear phase response error while satisfying the stability constraints. By applying fuzzy set theory, multiobjective constrained problem is converted into multiobjective unconstrained optimization problem and is stated as below:(16)
Maximize
μ
1
X
,
μ
2
X
,
μ
3
X
T
,where μ
1
(
X
) is the membership function of magnitude response error given in (12); μ
2
(
X
) is the membership function of linear phase response error given in (13); μ
3
(
X
) is the membership function of stability constraints given in (15); X is a vector decision variable of dimensions D
×
1 with D
=
2
M
+
4
N
+
1.The objective is to find the value of filter coefficients being decision variables,X, which maximizes the entire objective functions simultaneously. The value of membership function indicates how much the solution satisfies the μ
i objective on the scale from 0 to 1. The maximum satisfaction of membership function for any filter coefficient combination is obtained by taking the intersection of the membership functions of participating objectives and is expressed as below:(17)
fit
j
t
=
min
⁡
μ
1
j
,
μ
2
j
,
μ
3
j
j
=
1,2
,
…
,
N
P
.


## 2.5. Optimal Order of Digital Filter

The expression for optimal order of digital IIR filter is described below:(18)
O
=
∑
j
=
0
M
α
j
+
2
∑
k
=
0
N
β
k
,where α
j and β
k are j
t
h and k
t
h control genes of corresponding first order and second order blocks, respectively, and the value of genes will be either 1 or 0. M and N are the maximum number of first order and second order blocks, respectively. Maximum order of the digital IIR filter will be M
+
2
N. The order of the filter has been determined by control genes as shown in Figure 4. The coding method has been taken from Yu and Xinjie [16]. The value of control gene determines whether the particular block will be considered for filter design or not. The block will be considered activated when the corresponding control gene is 1. The number of binary bits used to generate control genes depends upon the value of M and N. The decision vector x shown in (2) is modified as below:(19)
x
=
U
,
V
,
A
,
a
1
,
b
1
,
…
,
a
M
,
b
M
,
p
1
,
q
1
,
r
1
,
s
1
,
…
,
p
N
,
q
N
,
r
N
,
s
N
T
.Figure 4
Activation and deactivation of filter coefficients with control genes.The variableU is a positive integer with maximum value of (
2
M
-
1
) and the variable V is a positive integer with maximum value of (
2
N
-
1
). The variable U and V will also be optimized along with the coefficients of the filter.


## 3. Optimization Technique

### 3.1. Gravitational Search Algorithm

This optimization algorithm is based on the law of gravity. In this algorithm, each candidate of the population is considered as an object (mass). Each object has four specifications, that is, active gravitational mass, passive gravitational mass, inertia mass, and position. Gravitational and inertia masses are determined by the fitness function and position of the object corresponds to the solution of the problem. All the objects attract each other with a gravitational force. Candidates having good solutions have heavy masses. Heavy masses move slowly and attract lighter masses towards good solutions. At the end of iteration, masses are updated as per their new positions. With the lapse of time, we consider that all the masses be attracted by the heaviest mass and this mass will represent the optimum solution [22].In GSA, population is considered as an isolated universe of objects (masses). Value of gravitational mass and inertia mass of an object has been assumed equal. Each object obeys the Newtonian laws of gravitation and motion as described below.Law of Gravitation. Each object attracts every other object with gravitational force. The gravitational force between two objects is directly proportional to the product of their masses (M) and inversely proportional to the square of the distance (R) between them. But in GSA only distance R is taken instead ofR
2, because it gives better results. This is the deviation of GSA from the Newtonian laws of gravitation.Law of Motion. The current velocity of any object is equal to the fraction of previous velocity of the object and acceleration of the object. Acceleration of the object is equal to the force applied on the object divided by the inertia mass of the object.Algorithm has been initialized with random population ofD dimensional N
P objects (masses) and each object is describes as(20)
X
i
=
x
i
1
,
…
,
x
i
d
,
…
,
x
i
D
for
i
=
1,2
,
…
,
N
P
;
d
=
1,2
,
…
,
D
,where x
i
d represents the position of i
t
h object in the d
t
h dimension.At a specific iteration “t,” the force acting on i
t
h object by j
t
h object is described as below:(21)
F
i
j
d
t
=
G
t
M
i
t
×
M
j
t
R
i
j
t
+
ε
x
j
d
t
-
x
i
d
t
,where M
i and M
j are gravitational masses of i
t
h and j
t
h objects, respectively. G
(
t
) is the gravitational constant at iteration “t.” ε, a small constant, has been added to avoid extra ordinary high value of force between two objects and almost lies on the same place in the space. R
i
j is the Euclidian distance between i
t
h and j
t
h objects and described as below:(22)
R
i
j
t
=
X
i
t
,
X
j
t
2
.To provide stochastic characteristics to the algorithm, it is considered that total force acting on thei
t
h object in d
t
h dimension is a randomly weighted sum of d
t
h component of the force exerted by all other objects and is given by(23)
F
i
d
t
=
∑
j
=
1
,
j
≠
i
N
P
rand
⁡
j
F
i
j
d
t
,where r
a
n
d
j is a random variable in the interval [
0,1
], corresponding to the j
t
h object.As per the law of motion, the acceleration of thei
t
h object in the d
t
h direction at iteration “t” is defined as below:(24)
a
i
d
t
=
F
i
d
t
M
i
t
.The velocity of the object is calculated as sum of fraction of current velocity and acceleration of the object. Therefore, the velocity and position of the object is updated as described below:(25)
v
i
d
t
+
1
=
rand
⁡
i
×
v
i
d
t
+
a
i
d
t
,
(26)
x
i
d
t
+
1
=
x
i
d
t
+
v
i
d
t
+
1
,where r
a
n
d
i is a random variable in the interval [
0,1
]. The gravitational constant G has been initialized at the beginning and reduces with successive iterations. Gravitational constant G
(
t
) at iteration “t” is described as below:(27)
G
t
=
G
o
e
-
α
t
/
T
,where G
o is the initial value of gravitational constant. α is a predefined constant and T is the maximum number of iterations.Masses of the objects are updated as follows:(28)
m
i
t
=
fit
i
t
-
worst
t
best
t
-
worst
t
,
M
i
t
=
m
i
t
∑
j
=
1
N
m
j
t
for
i
=
1,2
,
…
,
N
,where fit
i
t represents the fitness value of the i
t
h object at iteration “t.” For minimization problem, best
t and worst
(
t
) are described as follows:(29)
best
t
=
min
j
∈
1
,
…
,
N
P
⁡
fit
j
t
,
worst
t
=
max
j
∈
1
,
…
,
N
P
⁡
fit
j
t
.To have good compromise between exploration and exploitation, the number of objects can be reduced with successive iterations. So it is supposed that only a set of objects with heavy masses will exert gravitational force on other objects. But to avoid trapping in local minima, algorithm should use exploration in the beginning. By lapse of iterations, the exploration should fade out and exploitation should fade in. So onlykbest objects exert gravitational force on other objects. Initiallykbest is a taken equal to number of objects and all objects exert force. With successive iterationskbest decreases linearly in such way that at the last iteration only one object has been left and this provides us the optimum result. Therefore (23) can be modified as follows:(30)
F
i
d
t
=
∑
j
∈
k
b
e
s
t
,
j
≠
i
rand
⁡
j
F
i
j
d
t
,wherekbest is the set of objects with best fitness value and maximum mass.If elements of object velocityv
i
d
t violate their limits, their values are updated as below:(31)
v
i
d
t
=
v
d
m
i
n
;
v
i
d
t
<
v
d
m
i
n
v
d
m
a
x
;
v
i
d
t
>
v
d
m
a
x
v
i
d
t
;
no violation of limits.Similarly if elements of object positionx
i
d
t violate their limits, their values are updated as below:(32)
x
i
d
t
=
x
d
m
i
n
;
x
i
d
t
<
x
d
m
i
n
x
d
m
a
x
;
x
i
d
t
>
x
d
m
a
x
x
i
d
t
;
no violation of limits.


### 3.2. Exploratory Move

In exploratory move, the current point is perturbed in all possible directions for each and every variable at a time and the best point is recorded. After each perturbation, the present point is changed to the best point. At the end of all variable perturbations, if the point found is different from the original point, the exploratory move is called a success; otherwise, the exploratory move is a failure. In any case, the best point arrived at the end of exploratory move.The evolutionary method is used to search for the optimal filter coefficients. In this method, forD number of filter coefficients, 2
D feasible solutions are generated. A D dimensional hypercube of side Δ is formed around the point. x
i
C represents the coefficients of IIR filter from the current point of hypercube. The better solution is obtained from the objective function of the IIR filter. Another hypercube is formed around the current better point and the iterative process is continued. All the corners of the hypercube, represented in the D binary bits code, are explored for better results simultaneously. Table 1 shows the coefficient pattern for 3-coefficient digital IIR filter where 3-bit binary code is used to represent the 8 corners of the three-dimensional hypercube (Figure 5). The decimal serial numbers of the hypercube are changed into their respective binary codes. The deviation from the current center point is obtained by replacing 1’s of the code with +Δ and 0’s with −Δ. As the number of coefficients of the IIR filter increased, the number of hypercube corners increases exponentially. So the process becomes time consuming [23].Table 1
Coefficient vector at the corners of hypercube.
Hyper cube corners
Possible combinations of 3-bit
C
2
C
1 C
0
Distance of hypercube from centre point
x
i
3
c
,
x
i
2
c
,
x
i
1
c
Possible coefficient pattern of the IIR filter at the corner of hypercube
0
000
−Δ
i
3  −Δ
i
2  −Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
-
Δ
i
1
1
001
−Δ
i
3  −Δ
i
2  +Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
+
Δ
i
1
2
010
−Δ
i
3  +Δ
i
2  −Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
-
Δ
i
1
3
011
−Δ
i
3  +Δ
i
2  +Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
+
Δ
i
1
4
100
+Δ
i
3  −Δ
i
2  −Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
-
Δ
i
1
5
101
+Δ
i
3  −Δ
i
2  +Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
+
Δ
i
1
6
110
+Δ
i
3  +Δ
i
2  −Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
-
Δ
i
1
7
111
+Δ
i
3  +Δ
i
2  +Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
+
Δ
i
1Figure 5
Three-dimensional hypercube representing corners in decimals.


### 3.3. BSA Strategy

To reduce the computational time, binary successive approximation (BSA) strategy is used to explore the optimal solution. BSA strategy to search for the optimal solution is explained in Figure6, where solution procedure moves towards the optimal solution by comparing two solutions at a time represented by the two corners of the hypercube [23].Figure 6
BSA for 3-bit code.The search process is started by initializing the coefficient vectorx
j
C
t giving objective functionF
t. To perform BSA strategy by the iterative process C
j
t is initially selected as below:(33)
C
i
j
t
=
1
;
for
j
=
1
0
;
for
j
=
2,3
,
4
,
…
,
D
.Two corners, with reference to above selected corner, are created for comparison as below:(34)
C
i
1
j
t
=
1
;
for
i
+
1
C
i
j
t
;
for
j
=
1,2
,
…
,
i
,
i
+
2
,
…
,
D
,
(35)
C
i
2
j
t
=
0
;
for
i
C
i
1
j
t
;
for
j
=
1,2
,
…
,
i
-
1
,
i
+
1
,
…
,
D
.In reference to these two corners, coefficient vectors are generated as shown in Table 1. Mathematically, it is represented in the generalized form: (36)
x
i
m
j
t
=
x
i
j
C
t
+
Δ
i
m
j
t
m
=
1,2
;
j
=
1,2
,
…
,
D
;
i
=
1,2
,
…
,
N
P
,where(37)
Δ
i
m
j
t
=
+
Δ
i
j
if
C
i
m
j
t
=
1
-
Δ
i
j
if
C
i
m
j
t
=
0
m
=
1,2
;
j
=
1,2
,
…
,
D
;
i
=
1,2
,
…
,
N
P
.The initial increment to coefficients is decided by(38)
Δ
j
=
x
j
m
a
x
-
x
j
m
i
n
δ
.Objective functions atx
1
j
k and x
2
j
k are evaluated using (17) as follows:(39)
F
i
m
t
=
f
x
i
m
j
t
m
=
1,2
.The minimum value of these two is selected to be computed with the rest of the corners, generated subsequently, and the selected corner for the generation of the next two corner is(40)
C
i
j
t
=
C
i
1
j
t
if
F
i
1
t
<
F
i
2
t
C
i
2
j
t
if
F
i
1
t
>
F
i
2
t
j
=
1,2
,
…
,
D
.This process is repeated till all the corners of the hypercube are explored and the overall minimum is selected to find the new centre point for the next iteration. When the last element ofC
i
j
t vector contains one of the last branches of BSA tree that is reached, which ensures that all corners of hypercube are explored, the procedure is terminated. In BSA method, the number of computations is reduced by large amount, as elaborated in Table 2.Table 2
Comparison of number of function evaluations.
Value ofM and N
Number of committed coefficients
D
Number of corners of hypercube    (2
D)
Number of comparisons by BSA method  (2
×
D)
1, 1
7
128
14
2, 2
13
8192
26
3, 3
19
524,288
38
4, 4
25
33,554,432
50
5, 5
31
2,147,483,648
62


### 3.4. Algorithms

The different steps of the proposed algorithms are shown in Algorithms1 and 2.Algorithm 1:Hybrid GSA.
(1) Search space identification.
(2) Randomized initialization of population.
(3) Fitness evaluation of objects using  (17).
WHILE (stopping criterion is not met).
(4) UpdateG
(
t
), b
e
s
t
(
t
), w
o
r
s
t
(
t
) and M
i
(
t
) (i
=
1,2
,
…
,
N
P).
(5) Calculation of the total force in different directions using (23).
(6) Calculation of acceleration and velocity using (24) and (25) respectively.
(7) Updating objects’ position using (26).
(8) Fitness evaluation of objects using (17).
(9) Apply exploratory move to improve the fitness value using Algorithm2.
ENDDO
STOP.Algorithm 2:BSA.
(1) Enter with filter coefficients being decision variable asx
i
0, F
0
=
f
x
i
0 and set k
=
0.
DO
(2.1)x
i
c
=
x
i
k  (
i
=
1,2
,
…
,
D
).
(2.2) SetC
i
k  i
=
1,2
,
…
,
D using (33), k
=
k
+
1.
(2.3) ComputeC
1
i
k  i
=
1,2
,
…
,
D using  (34).
(2.4) ComputeC
2
i
k  i
=
1,2
,
…
,
D using  (35).
(2.5) Obtainx
m
i
k and Δ
m
i
k
(
m
=
i
,
2
)
(
i
=
1,2
,
…
,
D
) from  (36) and  (37) respectively.
(2.6) Evaluate fitnessF
k from  (39).
DO
(2.7.1)j
=
j
+
1
(2.7.2) ObtainC
m
i
k, x
m
i
k and Δ
m
i
k
(
m
=
i
,
2
)
(
i
=
1,2
,
…
,
D
) from  (34),  (35),  (36), and  (37) respectively.
(2.7.3) Evaluate fitnessF
k from  (39).
(2.7.4) ObtainC
i
k  i
=
1,2
,
…
,
D from  (40).
WHILE
j
<
D
IF
F
k
<
F
k
-
1
THEN
(2.8.1)x
i
k
=
x
i
k
-
1  i
=
1,2
,
…
,
D and           F
k
=
F
k
-
1
ELSE
(2.8.2)Δ
i
=
Δ
i
/
σ  i
=
1,2
,
…
,
D
END
WHILE
Δ
≤
ϵ
RETURN


## 3.1. Gravitational Search Algorithm

This optimization algorithm is based on the law of gravity. In this algorithm, each candidate of the population is considered as an object (mass). Each object has four specifications, that is, active gravitational mass, passive gravitational mass, inertia mass, and position. Gravitational and inertia masses are determined by the fitness function and position of the object corresponds to the solution of the problem. All the objects attract each other with a gravitational force. Candidates having good solutions have heavy masses. Heavy masses move slowly and attract lighter masses towards good solutions. At the end of iteration, masses are updated as per their new positions. With the lapse of time, we consider that all the masses be attracted by the heaviest mass and this mass will represent the optimum solution [22].In GSA, population is considered as an isolated universe of objects (masses). Value of gravitational mass and inertia mass of an object has been assumed equal. Each object obeys the Newtonian laws of gravitation and motion as described below.Law of Gravitation. Each object attracts every other object with gravitational force. The gravitational force between two objects is directly proportional to the product of their masses (M) and inversely proportional to the square of the distance (R) between them. But in GSA only distance R is taken instead ofR
2, because it gives better results. This is the deviation of GSA from the Newtonian laws of gravitation.Law of Motion. The current velocity of any object is equal to the fraction of previous velocity of the object and acceleration of the object. Acceleration of the object is equal to the force applied on the object divided by the inertia mass of the object.Algorithm has been initialized with random population ofD dimensional N
P objects (masses) and each object is describes as(20)
X
i
=
x
i
1
,
…
,
x
i
d
,
…
,
x
i
D
for
i
=
1,2
,
…
,
N
P
;
d
=
1,2
,
…
,
D
,where x
i
d represents the position of i
t
h object in the d
t
h dimension.At a specific iteration “t,” the force acting on i
t
h object by j
t
h object is described as below:(21)
F
i
j
d
t
=
G
t
M
i
t
×
M
j
t
R
i
j
t
+
ε
x
j
d
t
-
x
i
d
t
,where M
i and M
j are gravitational masses of i
t
h and j
t
h objects, respectively. G
(
t
) is the gravitational constant at iteration “t.” ε, a small constant, has been added to avoid extra ordinary high value of force between two objects and almost lies on the same place in the space. R
i
j is the Euclidian distance between i
t
h and j
t
h objects and described as below:(22)
R
i
j
t
=
X
i
t
,
X
j
t
2
.To provide stochastic characteristics to the algorithm, it is considered that total force acting on thei
t
h object in d
t
h dimension is a randomly weighted sum of d
t
h component of the force exerted by all other objects and is given by(23)
F
i
d
t
=
∑
j
=
1
,
j
≠
i
N
P
rand
⁡
j
F
i
j
d
t
,where r
a
n
d
j is a random variable in the interval [
0,1
], corresponding to the j
t
h object.As per the law of motion, the acceleration of thei
t
h object in the d
t
h direction at iteration “t” is defined as below:(24)
a
i
d
t
=
F
i
d
t
M
i
t
.The velocity of the object is calculated as sum of fraction of current velocity and acceleration of the object. Therefore, the velocity and position of the object is updated as described below:(25)
v
i
d
t
+
1
=
rand
⁡
i
×
v
i
d
t
+
a
i
d
t
,
(26)
x
i
d
t
+
1
=
x
i
d
t
+
v
i
d
t
+
1
,where r
a
n
d
i is a random variable in the interval [
0,1
]. The gravitational constant G has been initialized at the beginning and reduces with successive iterations. Gravitational constant G
(
t
) at iteration “t” is described as below:(27)
G
t
=
G
o
e
-
α
t
/
T
,where G
o is the initial value of gravitational constant. α is a predefined constant and T is the maximum number of iterations.Masses of the objects are updated as follows:(28)
m
i
t
=
fit
i
t
-
worst
t
best
t
-
worst
t
,
M
i
t
=
m
i
t
∑
j
=
1
N
m
j
t
for
i
=
1,2
,
…
,
N
,where fit
i
t represents the fitness value of the i
t
h object at iteration “t.” For minimization problem, best
t and worst
(
t
) are described as follows:(29)
best
t
=
min
j
∈
1
,
…
,
N
P
⁡
fit
j
t
,
worst
t
=
max
j
∈
1
,
…
,
N
P
⁡
fit
j
t
.To have good compromise between exploration and exploitation, the number of objects can be reduced with successive iterations. So it is supposed that only a set of objects with heavy masses will exert gravitational force on other objects. But to avoid trapping in local minima, algorithm should use exploration in the beginning. By lapse of iterations, the exploration should fade out and exploitation should fade in. So onlykbest objects exert gravitational force on other objects. Initiallykbest is a taken equal to number of objects and all objects exert force. With successive iterationskbest decreases linearly in such way that at the last iteration only one object has been left and this provides us the optimum result. Therefore (23) can be modified as follows:(30)
F
i
d
t
=
∑
j
∈
k
b
e
s
t
,
j
≠
i
rand
⁡
j
F
i
j
d
t
,wherekbest is the set of objects with best fitness value and maximum mass.If elements of object velocityv
i
d
t violate their limits, their values are updated as below:(31)
v
i
d
t
=
v
d
m
i
n
;
v
i
d
t
<
v
d
m
i
n
v
d
m
a
x
;
v
i
d
t
>
v
d
m
a
x
v
i
d
t
;
no violation of limits.Similarly if elements of object positionx
i
d
t violate their limits, their values are updated as below:(32)
x
i
d
t
=
x
d
m
i
n
;
x
i
d
t
<
x
d
m
i
n
x
d
m
a
x
;
x
i
d
t
>
x
d
m
a
x
x
i
d
t
;
no violation of limits.


## 3.2. Exploratory Move

In exploratory move, the current point is perturbed in all possible directions for each and every variable at a time and the best point is recorded. After each perturbation, the present point is changed to the best point. At the end of all variable perturbations, if the point found is different from the original point, the exploratory move is called a success; otherwise, the exploratory move is a failure. In any case, the best point arrived at the end of exploratory move.The evolutionary method is used to search for the optimal filter coefficients. In this method, forD number of filter coefficients, 2
D feasible solutions are generated. A D dimensional hypercube of side Δ is formed around the point. x
i
C represents the coefficients of IIR filter from the current point of hypercube. The better solution is obtained from the objective function of the IIR filter. Another hypercube is formed around the current better point and the iterative process is continued. All the corners of the hypercube, represented in the D binary bits code, are explored for better results simultaneously. Table 1 shows the coefficient pattern for 3-coefficient digital IIR filter where 3-bit binary code is used to represent the 8 corners of the three-dimensional hypercube (Figure 5). The decimal serial numbers of the hypercube are changed into their respective binary codes. The deviation from the current center point is obtained by replacing 1’s of the code with +Δ and 0’s with −Δ. As the number of coefficients of the IIR filter increased, the number of hypercube corners increases exponentially. So the process becomes time consuming [23].Table 1
Coefficient vector at the corners of hypercube.
Hyper cube corners
Possible combinations of 3-bit
C
2
C
1 C
0
Distance of hypercube from centre point
x
i
3
c
,
x
i
2
c
,
x
i
1
c
Possible coefficient pattern of the IIR filter at the corner of hypercube
0
000
−Δ
i
3  −Δ
i
2  −Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
-
Δ
i
1
1
001
−Δ
i
3  −Δ
i
2  +Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
+
Δ
i
1
2
010
−Δ
i
3  +Δ
i
2  −Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
-
Δ
i
1
3
011
−Δ
i
3  +Δ
i
2  +Δ
i
1
x
i
3
c
-
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
+
Δ
i
1
4
100
+Δ
i
3  −Δ
i
2  −Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
-
Δ
i
1
5
101
+Δ
i
3  −Δ
i
2  +Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
-
Δ
i
2
x
i
1
c
+
Δ
i
1
6
110
+Δ
i
3  +Δ
i
2  −Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
-
Δ
i
1
7
111
+Δ
i
3  +Δ
i
2  +Δ
i
1
x
i
3
c
+
Δ
i
3
x
i
2
c
+
Δ
i
2
x
i
1
c
+
Δ
i
1Figure 5
Three-dimensional hypercube representing corners in decimals.


## 3.3. BSA Strategy

To reduce the computational time, binary successive approximation (BSA) strategy is used to explore the optimal solution. BSA strategy to search for the optimal solution is explained in Figure6, where solution procedure moves towards the optimal solution by comparing two solutions at a time represented by the two corners of the hypercube [23].Figure 6
BSA for 3-bit code.The search process is started by initializing the coefficient vectorx
j
C
t giving objective functionF
t. To perform BSA strategy by the iterative process C
j
t is initially selected as below:(33)
C
i
j
t
=
1
;
for
j
=
1
0
;
for
j
=
2,3
,
4
,
…
,
D
.Two corners, with reference to above selected corner, are created for comparison as below:(34)
C
i
1
j
t
=
1
;
for
i
+
1
C
i
j
t
;
for
j
=
1,2
,
…
,
i
,
i
+
2
,
…
,
D
,
(35)
C
i
2
j
t
=
0
;
for
i
C
i
1
j
t
;
for
j
=
1,2
,
…
,
i
-
1
,
i
+
1
,
…
,
D
.In reference to these two corners, coefficient vectors are generated as shown in Table 1. Mathematically, it is represented in the generalized form: (36)
x
i
m
j
t
=
x
i
j
C
t
+
Δ
i
m
j
t
m
=
1,2
;
j
=
1,2
,
…
,
D
;
i
=
1,2
,
…
,
N
P
,where(37)
Δ
i
m
j
t
=
+
Δ
i
j
if
C
i
m
j
t
=
1
-
Δ
i
j
if
C
i
m
j
t
=
0
m
=
1,2
;
j
=
1,2
,
…
,
D
;
i
=
1,2
,
…
,
N
P
.The initial increment to coefficients is decided by(38)
Δ
j
=
x
j
m
a
x
-
x
j
m
i
n
δ
.Objective functions atx
1
j
k and x
2
j
k are evaluated using (17) as follows:(39)
F
i
m
t
=
f
x
i
m
j
t
m
=
1,2
.The minimum value of these two is selected to be computed with the rest of the corners, generated subsequently, and the selected corner for the generation of the next two corner is(40)
C
i
j
t
=
C
i
1
j
t
if
F
i
1
t
<
F
i
2
t
C
i
2
j
t
if
F
i
1
t
>
F
i
2
t
j
=
1,2
,
…
,
D
.This process is repeated till all the corners of the hypercube are explored and the overall minimum is selected to find the new centre point for the next iteration. When the last element ofC
i
j
t vector contains one of the last branches of BSA tree that is reached, which ensures that all corners of hypercube are explored, the procedure is terminated. In BSA method, the number of computations is reduced by large amount, as elaborated in Table 2.Table 2
Comparison of number of function evaluations.
Value ofM and N
Number of committed coefficients
D
Number of corners of hypercube    (2
D)
Number of comparisons by BSA method  (2
×
D)
1, 1
7
128
14
2, 2
13
8192
26
3, 3
19
524,288
38
4, 4
25
33,554,432
50
5, 5
31
2,147,483,648
62


## 3.4. Algorithms

The different steps of the proposed algorithms are shown in Algorithms1 and 2.Algorithm 1:Hybrid GSA.
(1) Search space identification.
(2) Randomized initialization of population.
(3) Fitness evaluation of objects using  (17).
WHILE (stopping criterion is not met).
(4) UpdateG
(
t
), b
e
s
t
(
t
), w
o
r
s
t
(
t
) and M
i
(
t
) (i
=
1,2
,
…
,
N
P).
(5) Calculation of the total force in different directions using (23).
(6) Calculation of acceleration and velocity using (24) and (25) respectively.
(7) Updating objects’ position using (26).
(8) Fitness evaluation of objects using (17).
(9) Apply exploratory move to improve the fitness value using Algorithm2.
ENDDO
STOP.Algorithm 2:BSA.
(1) Enter with filter coefficients being decision variable asx
i
0, F
0
=
f
x
i
0 and set k
=
0.
DO
(2.1)x
i
c
=
x
i
k  (
i
=
1,2
,
…
,
D
).
(2.2) SetC
i
k  i
=
1,2
,
…
,
D using (33), k
=
k
+
1.
(2.3) ComputeC
1
i
k  i
=
1,2
,
…
,
D using  (34).
(2.4) ComputeC
2
i
k  i
=
1,2
,
…
,
D using  (35).
(2.5) Obtainx
m
i
k and Δ
m
i
k
(
m
=
i
,
2
)
(
i
=
1,2
,
…
,
D
) from  (36) and  (37) respectively.
(2.6) Evaluate fitnessF
k from  (39).
DO
(2.7.1)j
=
j
+
1
(2.7.2) ObtainC
m
i
k, x
m
i
k and Δ
m
i
k
(
m
=
i
,
2
)
(
i
=
1,2
,
…
,
D
) from  (34),  (35),  (36), and  (37) respectively.
(2.7.3) Evaluate fitnessF
k from  (39).
(2.7.4) ObtainC
i
k  i
=
1,2
,
…
,
D from  (40).
WHILE
j
<
D
IF
F
k
<
F
k
-
1
THEN
(2.8.1)x
i
k
=
x
i
k
-
1  i
=
1,2
,
…
,
D and           F
k
=
F
k
-
1
ELSE
(2.8.2)Δ
i
=
Δ
i
/
σ  i
=
1,2
,
…
,
D
END
WHILE
Δ
≤
ϵ
RETURN


## 4. Validation of Proposed Technique

The proposed GSA technique alone and then hybridized with BSA based evolutionary technique has been applied to standard unimodal test functions, multimodal test functions with large and small dimensional variables shown in Tables3, 4, and 5, respectively. The population has been taken as 50. Maximum number of iterations has been set to 200. In (27), G
o, the initial value of gravitational constant has been taken as 200. A predefined constant, α, has been taken as 20. The GSA is hybridized with BSA based on evolutionary search technique for searching the variable space locally. The exploratory move has been repeated 50 times. The algorithm is made to run 500 times independently to justify the global solution.Table 3
Unimodal test functions.
Name
Functions
Dimension
Search range
Sphere
f
01
=
∑
i
=
1
D
x
i
2
30
-
100,100
D
SumSquare
f
02
=
∑
i
=
1
D
i
x
i
2
30
-
10,10
D
SumPower
f
03
=
∑
i
=
1
D
x
i
(
i
+
1
)
30
-
1,1
D
Step
f
04
=
∑
i
=
1
D
x
i
+
0.5
2
30
-
100,100
D
Elliptic
f
05
=
∑
i
=
1
D
x
i
2
10
6
i
-
1
/
D
-
1
30
-
5,10
D
Exponential
f
06
=
exp
⁡
0.5
∑
i
=
1
D
x
i
-
1
30
-
1.28,1.28
D
Quartic
f
07
=
∑
i
=
1
D
i
x
i
4
+
r
a
n
d
o
m
0,1
30
-
1.28,1.28
D
Schwelfel_2.22
f
08
=
∑
i
=
1
D
x
i
+
∏
i
=
1
D
x
i
30
-
10,10
D
Schwelfel_1.2
f
09
=
∑
i
=
1
D
∑
j
=
1
i
x
j
2
30
-
100,100
D
Schwelfel_2.21
f
10
=
max
⁡
x
1
,
1
≤
i
≤
D
30
-
100,100
DTable 4
Multimodal test functions with large dimensional variables.
Name
Functions
Dimension
Search range
Ackly
f
11
=
20
exp
⁡
-
0.2
D
-
1
∑
i
=
1
D
x
i
2
-
exp
⁡
D
-
1
∑
i
=
1
D
cos
⁡
2
π
x
i
+
20
+
e
30
-
32,32
D
Greiwangk
f
12
=
∑
i
=
1
D
x
i
2
4000
-
∏
i
=
1
D
cos
⁡
x
i
i
+
1
30
-
600,600
D
Penalized-1
f
13
=
π
D
∑
i
=
1
D
-
1
y
i
-
1
2
1
+
10
sin
2
⁡
π
y
i
+
1
+
10
sin
2
⁡
π
y
1
+
y
D
-
1
2
+
∑
i
=
1
D
u
x
i
,
10,100,4
30
-
50,50
D
Penalized-2
f
14
=
0.1
∑
i
=
1
D
-
1
x
i
-
1
2
1
+
sin
2
⁡
3
π
x
i
+
1
+
sin
2
⁡
3
π
x
1
+
x
D
-
1
2
1
+
sin
2
⁡
2
π
x
D
+
∑
i
=
1
D
u
x
i
,
5,100,4
30
-
50,50
D
Alpine
f
15
=
∑
i
=
1
D
x
i
sin
⁡
x
i
+
0.1
x
i
30
-
10,10
D
Bohachevsky_2
f
16
=
∑
i
=
1
D
-
1
x
i
2
+
2
x
i
+
1
2
-
0.3
cos
⁡
3
π
x
i
cos
⁡
3
π
x
i
+
1
+
0.3
30
-
100,100
DTable 5
Multimodal test functions with small dimensional variables.
Name
Functions
Dimension
Search range
Foxholes
f
17
=
1
500
+
∑
j
=
1
25
j
+
∑
i
=
1
2
x
i
-
a
i
j
6
-
1
-
1
2
-
65.536
,
65.536
D
Kowalik
f
18
=
∑
i
=
1
11
a
i
-
x
1
b
i
2
+
b
i
x
2
b
i
2
+
b
i
x
3
+
x
4
4
-
5,5
D
Six-hump Camel-back
f
19
=
4
x
1
2
-
2.1
x
1
4
+
1
3
x
1
6
+
x
1
x
2
-
4
x
2
2
+
4
x
2
4
2
-
5,5
D
Branin
f
20
=
x
2
-
5.1
4
π
2
x
1
2
+
5
π
x
1
2
+
10
1
-
1
8
π
cos
⁡
x
1
+
10
2
[
-
5,10
], [
0,15
]
Goldstien-price
f
21
=
1
+
x
1
+
x
2
+
1
2
19
-
14
x
1
+
3
x
1
2
-
14
x
2
+
6
x
1
x
2
+
3
x
2
2    ×
30
+
2
x
1
-
3
x
2
2
18
-
32
x
1
+
12
x
1
2
+
48
x
2
-
32
x
1
x
2
+
27
x
2
2
2
-
2,2
D
Hartman-1
f
22
=
-
∑
i
=
1
4
c
i
exp
⁡
-
∑
j
=
1
4
a
i
j
x
j
-
p
i
j
2
4
0,1
D
Hartman-2
f
23
=
-
∑
i
=
1
4
c
i
exp
⁡
-
∑
j
=
1
6
a
i
j
x
j
-
p
i
j
2
6
0,1
D
Shekel’s-1
f
24
=
-
∑
i
=
1
5
X
-
a
i
X
-
a
i
T
+
c
i
-
1
4
0,10
D
Shekel’s-2
f
25
=
-
∑
i
=
1
7
X
-
a
i
X
-
a
i
T
+
c
i
-
1
4
0,10
D
Shekel’s-3
f
26
=
-
∑
i
=
1
10
X
-
a
i
X
-
a
i
T
+
c
i
-
1
4
0,10
DThe desired and obtained values of test functions are shown in Table6. It is observed that for test functions with large dimensional variables, hybridization of GSA with BSA technique gives better results as compared to GSA in terms of achieved function values.Table 6
Results for test functions.
Function
Desired optimal value of test function
GSA
GSA-BSA
Obtained value of test function
Standard deviation
Obtained value of test function
Standard deviation
f
01
0
6.15
×
10
-
15
3.38
×
10
-
01
9.94
×
10
-
16
4.98
×
10
-
01
f
02
0
1.26
×
10
-
13
2.91
×
10
-
05
1.27
×
10
-
14
1.35
×
10
-
05
f
03
0
1.57
×
10
-
13
1.32
×
10
-
10
1.10
×
10
-
13
1.33
×
10
-
10
f
04
0
0
1.56
×
10
+
00
0
8.59
×
10
-
01
f
05
0
1.24
×
10
+
02
1.18
×
10
+
02
4.24
×
10
+
01
1.08
×
10
+
02
f
06
0
0
7.86
×
10
-
08
0
2.02
×
10
-
08
f
07
0
1.25
×
10
-
06
1.41
×
10
-
04
8.15
×
10
-
17
4.57
×
10
-
06
f
08
0
3.44
×
10
-
07
5.15
×
10
-
04
1.34
×
10
-
07
2.27
×
10
-
08
f
09
0
1.13
×
10
-
12
4.27
×
10
+
00
2.30
×
10
-
14
1.03
×
10
+
01
f
10
0
2.28
×
10
-
02
9.58
×
10
-
01
1.73
×
10
-
01
1.11
×
10
+
00
f
11
0
7.36
×
10
-
01
7.04
×
10
-
02
2.87
×
10
-
01
1.73
×
10
-
01
f
12
0
1.24
×
10
+
00
3.81
×
10
-
01
1.10
×
10
+
00
2.15
×
10
-
01
f
13
0
1.65
×
10
-
02
1.60
×
10
-
01
1.18
×
10
-
02
1.10
×
10
-
01
f
14
0
2.36
×
10
-
01
1.00
×
10
+
00
2.20
×
10
-
01
9.32
×
10
-
01
f
15
0
7.67
×
10
-
04
1.21
×
10
-
02
1.71
×
10
-
08
1.21
×
10
-
02
f
16
0
0
4.40
×
10
+
00
0
2.49
×
10
+
00
f
17
1.0
0.9980
6.28
×
10
-
01
0.9980
6.17
×
10
-
01
f
18
0.000308
0.000370
7.35
×
10
-
05
0.000364
7.55
×
10
-
05
f
19
−1.0316
−1.0316
6.60
×
10
-
05
−1.0316
1.28
×
10
-
04
f
20
0.3979
0.3979
5.58
×
10
-
17
0.3979
5.58
×
10
-
17
f
21
3.0
3.0
1.34
×
10
-
15
3.0
1.34
×
10
-
15
f
22
−3.8628
−3.8628
8.93
×
10
-
15
−3.8628
8.93
×
10
-
15
f
23
−3.3224
−3.3852
1.25
×
10
-
02
−3.3852
1.25
×
10
-
02
f
24
−10.1532
−10.1532
2.09
×
10
+
00
−10.1532
2.48
×
10
+
00
f
25
−10.4029
−10.4029
2.81
×
10
-
01
−10.4029
1.79
×
10
-
14
f
26
−10.5364
−10.5364
1.61
×
10
-
14
−10.5364
1.61
×
10
-
14Wilcoxon Signed Rank test has been applied to validate the performance of hybrid algorithm at a significance level 0.05.Result 1 (Z-value).
Z-values are expressed in terms of mean and standard deviations of test statics. It is observed that Z-value is −2.5424 and its corresponding p value is 0.01108. So the result is significant at p
≤
0.05.Result 2 (W-value).
TheW-value is 12. The critical value of W for N
=
14 at p
≤
0.05 is 21. Therefore, the result is significant at p
≤
0.05.From the above results, it is concluded that proposed GSA and hybrid GSA can be applied to unimodal as well as multimodal test functions. Since error surface of IIR filter is a multimodal function, so the above proposed technique can be applied to design the digital IIR filters.


## 5. Design of Digital IIR Filter and Comparisons

The design of digital IIR filter in cascaded form has been implemented using proposed GSA and hybrid GSA (HGSA) techniques by searching the filter coefficients in such a manner so that the fitness value (17) approaches 1. The performance of the designed digital IIR filter is measured based on pass-band and stop-band ripples, phase response error, and order of the filter. Low-pass (LP), high-pass (HP), band-pass (BP), and band-stop (BS) IIR filters have been considered for the design. In this paper, the order of digital IIR filter is a variable in the optimization process and is optimized simultaneously along with μ
1, μ
2, and μ
3 objective functions. The maximum order for LP, HP, BP, and BS filters has remained 12 as shown in Table 7. Hence the maximum value of M and N is kept as 4 for LP and HP filters and 0 and 6, respectively, for BP and BS filters. The design conditions for pass-band and stop-band normalized frequencies of LP, HP, BP, and BS filters are also shown in Table 7, where ω is the normalized frequency of the signal and varies from 0 to π. The results of the digital IIR filter design given by Yu and Xinjie [16], Wang et al. [17], and Kaur et al. [18] are referred to, to compare with design obtained by proposed HGSA approach.Table 7
Prescribed design conditions for LP, HP, BP, and BS filters.
Filter type
Maximum order of filter
Pass band
Stop band
LP
12
0
≤
ω
≤
0.2
π
0.3
π
≤
ω
≤
π
HP
12
0.8
π
≤
ω
≤
π
0
≤
ω
≤
0.7
π
BP
12
0.4
π
≤
ω
≤
0.6
π
0
≤
ω
≤
0.25
π
0.75
π
≤
ω
≤
π
BS
12
0
≤
ω
≤
0.25
π
0.75
π
≤
ω
≤
π
0.4
π
≤
ω
≤
0.6
πThe phase response error (8) pass-band as well as stop-band ripples (4a) and (4b) obtained from the proposed GSA and HGSA techniques for LP, HP, BP, and BS filters are compared with CCGA, NSGA-II, LS-MOEA, and RCGA in Table 8. From Table 8, it is concluded that the proposed GSA and HGSA techniques offer better performance in terms of phase response error, pass-band ripples as well as stop-band ripples for LP, HP, BP, and BS filters. It is also revealed HGSA designs have better IIR filter than GSA. The filter design is performed for 500 independent trial runs to achieve minimum, maximum, average, and standard deviation of fitness function described in (17) and is depicted in Table 9. Very small value of standard deviation proves the robustness of the proposed hybrid search technique to achieve global solution.Table 8
Comparison of design results for LP, HP, BP, and BS Filters.
Technique
Lowest order of filter
Pass-band magnitude ripple
Stop-band magnitude ripple
Phase response error
LP filter
CCGA
3
0.9034
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1699
1.4749
×
10
-
4
NSGA-II
3
0.9117
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1719
1.2662
×
10
-
4
LS-MOEA
3
0.9083
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1586
1.0959
×
10
-
4
RCGA
3
0.9141
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1556
1.1788
×
10
-
4
GSA
3
0.9201
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1471
1.0391
×
10
-
4
HGSA
3
0.9208
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1475
1.0375
×
10
-
4
HP filter
CCGA
3
0.9044
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1749
9.7746
×
10
-
5
NSGA-II
3
0.8960
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1769
9.1419
×
10
-
5
LS-MOEA
3
0.9004
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1746
9.6251
×
10
-
5
RCGA
3
0.9004
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1742
9.5757
×
10
-
5
GSA
3
0.9052
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1674
7.5431
×
10
-
5
HGSA
3
0.9051
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1698
7.3734
×
10
-
5
BP filter
CCGA
4
0.8920
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1654
8.1751
×
10
-
5
NSGA-II
4
0.9100
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1771
3.6503
×
10
-
4
LS-MOEA
4
0.9285
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1734
6.0371
×
10
-
5
RCGA
4
0.9333
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1641
5.9070
×
10
-
5
GSA
4
0.9354
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1587
6.0121
×
10
-
5
HGSA
4
0.9398
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1634
5.6440
×
10
-
5
BS filter
CCGA
4
0.8966
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1733
1.6198
×
10
-
4
NSGA-II
4
0.8917
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1770
1.5190
×
10
-
4
LS-MOEA
4
0.8967
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1725
1.5084
×
10
-
4
RCGA
4
0.8975
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1708
1.5144
×
10
-
4
GSA
4
0.9004
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1481
1.4103
×
10
-
4
HGSA
4
0.9009
≤
H
(
e
j
ω
)
≤
1.0
H
(
e
j
ω
)
≤
0.1459
1.3827
×
10
-
4Table 9
Variation of fitness function for LP, HP, BP and BS filters.
Filter
Minimum
Maximum
Average
Standard Deviation
LP
0.9948596
0.9998963
0.9998357
3.6140
×
10
-
4
HP
0.7615849
0.9999232
0.8609298
8.4818
×
10
-
2
BP
0.9900417
0.9999436
0.9998430
7.3827
×
10
-
4
BS
0.9993995
0.9998618
0.9998390
2.2359
×
10
-
5The designed filters obtained by the proposed HGSA technique for LP, HP, BP, and BS are given by (41), (42), (43), and (44), respectively.The magnitude response and phase response of LP, HP, BP, and BS are shown in Figure7. The pole-zero plots for LP, HP, BP, and BS are shown in Figure 8. It is observed that the designed filters follow the stability constraints applied during designing process as all the poles lie inside the unit circle.Figure 7
Magnitude and phase responses of LP, HP, BP, and BS filters.Figure 8
Pole and zero plots of LP, HP, BP, and BS filters.From the above results, it is concluded that the proposed hybrid heuristic search technique is effective and robust technique for design of digital IIR filter with better phase response error, pass-band ripples, and stop-band ripples:(41)
H
LP
z
=
0.201016
×
z
+
0.216875
z
-
0.406875
×
z
2
-
0.949082
z
+
0.959715
z
2
-
1.230432
z
+
0.649521
,
(42)
H
HP
z
=
0.149117
×
z
-
0.597477
z
+
0.402310
×
z
2
+
0.877094
z
+
0.918969
z
2
+
1.220268
z
+
0.634219
,
(43)
H
BP
z
=
0.238985
×
z
2
-
1.608518
z
+
0.968449
z
2
+
0.622686
z
+
0.530151
×
z
2
+
1.616206
z
+
0.984543
z
2
-
0.638647
z
+
0.527260
,
(44)
H
BS
z
=
0.455623
×
z
2
+
0.358866
z
+
0.913907
z
2
-
0.738038
z
+
0.471956
×
z
2
-
0.363914
z
+
0.943718
z
2
+
0.741348
z
+
0.475340
.


## 6. Conclusion

This paper proposes the design of digital IIR filter with conflicting objectives by using hybrid GSA where GSA and BSA based evolutionary search algorithms are applied as global and local search optimization techniques, respectively. HGSA maintains the balance between exploration and exploitation. Lesser parameters are required to tune. The results reveal that the proposed HGSA techniques offer better performance in terms of phase response error and pass-band as well as stop-band ripples for LP, HP, BP, and BS filters as compared to CCGA and NSGA-II, LS-MOEA, and RCGA methods. It is also revealed that HGSA outperforms GSA to optimize unimodal, multimodal higher, and lower dimensions standard test problems as well as for the design of digital IIR filters. The proposed method independently designs digital filter, whether it is LP, HP, BP, or BS filter.



---

*Source: 282809-2015-11-05.xml*