A Fuzzy Gravitational Search Algorithm to Design Optimal IIR Filters

: The goodness of Inﬁnite Impulse Response (IIR) digital ﬁlters design depends on pass band ripple, stop band ripple and transition band values. The main problem is deﬁning a suitable error ﬁtness function that depends on these parameters. This ﬁtness function can be optimized by search algorithms such as evolutionary algorithms. This paper proposes an intelligent algorithm for the design of optimal 8th order IIR ﬁlters. The main contribution is the design of Fuzzy Inference Systems able to tune key parameters of a revisited version of the Gravitational Search Algorithm (GSA). In this way, a Fuzzy Gravitational Search Algorithm (FGSA) is designed. The optimization performances of FGSA are compared with those of Differential Evolution (DE) and GSA. The results show that FGSA is the algorithm that gives the best compromise between goodness, robustness and convergence rate for the design of 8th order IIR ﬁlters. Moreover, FGSA assures a good stability of the designed ﬁlters.


Introduction
The design of optimal Infinite Impulse Response (IIR) digital filters is a very interesting challenge. The main techniques to design IIR filters are traditional design technique and optimization techniques. The first method is commonly known as bilinear transformation approach [1]. The second approach regards the applications of optimization techniques to design optimal filters. Among linear optimization algorithms, the steepest-descent and quasi-Newton (QN) algorithms are used for IIR filters design [1,2]. QN algorithms offer the advantages of robustness and fast convergence. Moreover, because QN optimization approach is very flexible, it can be used to design filters with arbitrary amplitude and/or phase responses. QN algorithms have also been used to design linear-phase IIR filters [3]. Chen et al. [4] proposed a technique for IIR filters design based on the minimization the error between the order-reduced filter's response and the desired one in the Hankel-norm sense. In the optimization algorithm proposed by Lu and Hinamoto [5] for the design of optimal IIR filters, the coefficients of all sub-filters are jointly optimized through a sequence of linear updates with each update carried out using second-order cone programming.
IIR filters are used in a wide range of applications where a high-selectivity processing of discrete signals is needed [6]. Lai and Lin [7] imposed two elliptic constraints on the frequency response of an IIR filter: the first one to minimize the maximum phase error, whereas the second one to constrain the maximum magnitude error. Another constrained optimization was introduced by Nongpiur et al. [8] for the design of IIR Digital Differentiators. The method in [8] minimizes the group-delay deviation under the constraint that the maximum amplitude-response error must be below a fixed level. Constraints on the magnitude and phase responses for the design of nearly linear-phase IIR filters was the main contribution of [9]. Lang [10] presented a method with the Algorithm (FGSA) with dynamic parameter adaptation. This algorithm is applied to design 8th order IIR filters.
The paper is organized as follows. Section 2 contains the description of IIR filters design. The designed algorithm is presented in Section 3. Section 4 illustrates the achieved results. Section 5 contains the paper conclusions.

IIR Filter Design
The relation between inputs and outputs of IIR filters is given by Equation (1) [44]: where x(p) is the filter input, whereas y(p) is the output. The order of filters is defined by n with n ≥ m. By assuming that a 0 = 0, the transfer function of IIR filter can be expressed as in Equation (2): Assuming that z = e jω , it follows that IIR filter frequency response becomes as in Equation (3) H where ω ∈ [0, π] is the digital frequency. Generally, the used approach to design IIR filters is to consider a Mean Square Error optimization problem [45][46][47]. MSE fitness function can be expressed as in Equation (4) where N s is the number of frequency points used for the computation of the error fitness function; d(p) and y(p) are the filter's desired and actual responses, respectively. The actual response y(p) is calculated through Equation (1), whereas the values of d(s) are set to be very close to ideal filters' values. The difference between d(p) and y(p) is the error between the desired and the actual filter responses. The design goal is to minimize the MSE J1(ω) with proper adjustment of filter coefficients b 0 , ..., b m , a 0 , ..., a n . IIR filters' optimization problem depends on the choice of transfer function coefficients b 0 , b 1 , ..., b m and a 0 , a 1 , ..., a n . Because the quality of a IIR filter depends on pass band ripple, stop band ripple and transition band, we propose a new fitness function that takes into account these three parameters. In particular, a good IIR filter has small pass and stop band ripple, and narrow transition band. In order to assure such constraints, the fitness function J 2 (ω) in Equation (5) is defined. In Equation (5), n s is the number of samples, δ p is the pass band ripple, δ s is the stop band ripple and |H(ω)| is the absolute value of H(ω) in Equation (3): In particular, the absolute value of H(ω) is calculated for each i = 1, ..., n s with (6) where the values of ω i are spaced frequency points between 0 and the pass band normalized edge frequency ω p . Similarly, for each l = 1, ..., n s where the values of ω l are spaced frequency points between 0 and the stop band normalized edge frequency ω s . Note that, when i = n s , it follows that ω i = ω n s = ω p . In the same way, the stop band normalized edge frequency ω s is achieved when l = n s ; thus, ω l = ω n s = ω s . In order to design optimal IIR filters, a constrained minimization of the error fitness function defined in Equation (5) is needed. On the other hand, the stability is an important issue for IIR digital filters design [10,48,49]. Jiang and Kwan [50] proposed a stability constraint with a prescribed pole radius derived from the argument principle of complex analysis. The optimization of the proposed fitness function defined in Equation (5) follows the stability constraints in [50].

The Fuzzy Gravitational Search Algorithm
The proposed error fitness function in Equation (5) has to be minimized through an optimization algorithm. The design of our algorithm starts with a suitable definition of GSA parameter, which supplies the number of agents that apply the force to other individuals [30]. Such GSA parameter is referred as Kbest and it decreases linearly to 1 over the increment of iterations. The idea is to increase the convergence speed by defining Kbest as in Equation (8) Kbest where i is the i-th iteration, n a is the number of agents and β a parameter.
A key parameter in GSA is the gravitational constant G(t) [30], which depends on the initial value G 0 , the number of iterations N and the value of parameter α (see Equation (9)): The next step is to design two Fuzzy Inference Systems (FIS) able to adjust β and α parameter in Equation (8) and Equation (9), respectively. The tuning of these parameters must assure a good trade-off between exploration and exploitation of the search process. For this aim, we define a quantity P p ∈ [0, 1], which gives a measure of the population progress (see Equation (10)), where J 2 (ω) (i) is the error fitness mean value calculated on the n a agents. This computation is accomplished for each iteration i = 1, ..., N, where N represents the iterations number. Generally, the accuracy of an FIS depends on the number of membership functions (MF): a higher membership functions number tends to cause an increase in FIS. Moreover, the number of MFs has a huge impact on the system complexity; therefore, an optimal trade-off between accuracy and complexity is needed. Thus, the fuzzy inputs definition for the first FIS depends on this issue. We define the iterations number N and the population progress P p as fuzzy inputs with membership functions Low (L), Medium (M) and High (H) (see Figures 1 and 2). In order to achieve more refined values of α, nine membership functions for the fuzzy output are defined (see Figure 3). The choice of triangular/trapezoidal MFs depends on the performance of FIS for a generic system [51]. The definition of fuzzy rules for FIS-α depends on GSA behavior. To assure exploration at the beginning of iterations, the agents must have a huge acceleration: this condition is achieved with low values of α (see Equation (9)). A lack of improvement means premature convergence: supplying a lower value of α, the agents escape from local optima. Low values of α are required when GSA lies in the middle of the procedure and there are no improvements. Very very high values of α are necessary with improvements at the end of iterations. Table 1 shows FIS-α rules. The architecture of FIS-α is shown in Figure 4.
The main step in the design of an FIS is the definition of fuzzy rules. Table 1 shows the fuzzy rules for FIS-α. The rules are based on the behavior of GSA. At the beginning of iterations, i.e., when N = L, to assure exploration, the agents must have a big acceleration; therefore, a low value of α is necessary. In fact, if α is low, then G increases (see Equation (9)), and then the acceleration tends to increase [30]. On the other hand, an early lack of improvement, that is P p = H and N = L, is a sign of premature convergence: with a very very low value of α, the individuals can escape from local optima. When GSA is at the middle of the procedure (N = M) and there is lack of improvement (P p = H), the values of α must be basically low. At the end of the iterations (N = H), if there is an improvement in the optima research (P p = L), then α must be very very high. High values of α tend to decrease the value of the gravitational constant and therefore the acceleration. Figure 4 shows the architecture of FIS-α to adjust the parameter α.
Pp (3) α (9) FuzzySystem (mamdani) 9 rules The second FIS (referred as FIS-β) has the task of adjusting the parameter β of quantity Kbest defined in Equation (8). This FIS has one input and one output as shown in Figure 5. The fuzzy input is the number of iterations N that has three triangular/trapezoidal membership functions: Low (L), Medium (M) and High (H) as in FIS-α (see Figure 1). In this way, FIS-β has three rules and thus we define three possible values of the fuzzy output β. Figure 6 shows the three triangular/trapezoidal membership functions Low (L), Medium (M) and High (H) of β. The definition of the fuzzy rules is based on the fact that there is a tendency for Kbest to decrease with higher iteration number. The base rule is shown in Table 2 and Figure 5 illustrates the architecture of FIS-β.
The target of the search procedure is to find the best coefficients b 0 , b 1 , ..., b m and a 0 , a 1 , ..., a n of H with m = n, where n is the order of IIR filter. In order to do this, Algorithm 1 is proposed. Both of the designed FISs are used in the optimization process.   S1. Initialize randomly the coefficients (b 0 , b 1 , ..., b n , a 0 , a 1 , ..., a n ) for each agent in the search space, where n is the order of IIR filter. Note that the number of coefficients is equal to 2(n + 1) and n a is the number of agents. Let up and low be the extrema of search interval, the matrix of coefficients for all agents X (n a ,2(n+1)) is computed by for each agent i ∈ [1, n a ] ∩ II N and coefficient j ∈ [1, 2(n + 1)] ∩ II N; rand(i, j) is a function that generates random numbers between 0 and 1. S2. Fix the extrema value α in f and α sup of fuzzy output α and the initial value for the first two iterations α(1) and α (2). S3. For each iteration from 1 to n a , execute the steps from S3.1 to S3.9. S3.1. Check the search space boundaries for agents according to the coefficients in X; agents that go out of the search space, are reinitialized randomly, i.e., for the i-th agent such that in the j-th coefficient x(i, j) > up or x(i, j) < low, compute x(i, j) by using Equation (11). S3.2. Evaluate the agents; compute the fitness J 2 (ω) of each agent by passing the coefficients of X to the test function and select the minimum fitness value among agents. S3.3. Calculate the masses of each agent, i.e., Equations (12)- (16), where M aj and M pi are the active and passive gravitational masses related to agent j and i, respectively, and M i (t) is the inertial mass of i-th agent at iteration t [30]: best(t) = min j∈1,...,n a S3.4. Compute the gravitational constant G (see Equation (9) Therefore, the new value of Kbest is computed by Equation (8). S3.6. According to Kbest, calculate the acceleration of each agent in gravitational field, see Equations (9) and (18)-(22), where d = 2(n + 1), with n filter order. S3.7. Compute the population progress P p ∈ [0, 1] as defined in Equation (10). S3.8. The values of normalized iteration number t n and population progress P p are passed as inputs to FIS-α, which gives a value of α between 0 and 1, denoted by α out . This value is normalized between α in f and α sup with the formula S3.9. Update the velocity v and coefficients in X of i-th agent, with the Equations (24) and where S4. Give in output the values of best coefficients (b 0 , b 1 , ..., b n , a 0 , a 1 , ..., a n ) opt .
FGSA designed in [35,36] have a more simple fuzzy rule base than FIS-α: the just one fuzzy input is the number of iterations. A similar approach has been introduced by Sombra [42] with a fuzzy system characterized by three fuzzy rules to adjust the parameter α of GSA. However, Algorithm 1 guarantees a high accuracy in the computation of α, with a minor adding of complexity compared to the mentioned approaches. On the other hand, Khabisi and Rashedi [43] designed an FGSA with 36 fuzzy rules, without achieving relevant results. The dynamic Type-2 fuzzy logic α adaptation proposed in [39] improves the convergence performances with an increase of the system complexity. The parameter Kbest has been adjusted by Olivas [38] directly as output fuzzy, whereas the proposed algorithm tunes the parameter β with a new definition of Kbest (see Equation (8)). This fact assures a better regulation of Kbest during the steps of our FGSA. Moreover, GSA with Wavelet Mutation proposed in [32] has the drawbacks of the rigorous trials required for the tuning of control parameters for the wavelet mutation method. Finally, Algorithm 1 supplies a good trade-off between accuracy and complexity compared with the mentioned approaches.

Experimental Results
Algorithm 1 is compared with the application of GSA and DE for the design of IIR filters. All optimization algorithms are run in the MATLAB environment (R2016a, MathWorks, Natick, MA, USA) on 2.20 GHz-speed processor. In particular, the fuzzy toolbox of MATLAB is exploited. This tool gives the possibility of designing Fuzzy Inference Systems based on the Mamdani inference method [52]. Because the MATLAB environment allows for referring to objects created by tools with the MATLAB code, the settings of the designed FISs can be dynamically modified. Moreover, the parameters of the fuzzy part of Algorithm 1 are computed by using the Center of the Mass defuzzification method. The MATLAB environment has been also exploited for its workspace data storage capability.
Referring to FGSA, the number of agents n a is set to 50 and the iterations number N is equal to 500. Moreover, α in f = 1, α sup = 20 and β in f = 2, β sup = 4 in Equations (23) and (17), respectively. The value of α in GSA is 20 (as in [30]), n a = 50 and N = 500. In DE, the population size is 50, the crossover probability is 0.2 and the maximum number of iterations is 500.
Equation (11) computes the values of coefficients in H: we assume a range from low = −2 to up = +2 with a filter order n = 8, which is b 0 , . .., b 8 , a 0 Table 3, with w frequency width. Moreover, the frequency range from 0 to π is divided into n s = 256 equally spaced sample points. DE, GSA and FGSA are run for 30 times to get the best solutions and the results in Tables 4-8 are the average (first sub-row) and standard deviation (second sub-row) on 30 experiments.       Table 4 shows that the proposed algorithm gives good results for LP filters because pass band ripple, stop band ripple and transition band are less than DE and GSA values. This fact is confirmed by comparing the magnitude response over normalized frequency of DE, GSA and FGSA (see Figure 7) for LP filters. By comparing FGSA and GSA, the pass band ripple (PB-ripple in the table) is improved by 3%, the stop band ripple (SB-ripple) is reduced by 7% and the transition band (tb) is improved by 16%. Moreover, the robustness of the proposed approach is analyzed by computing the standard deviation on 30 experiments. In fact, the robustness depends on standard deviation values: a result is more robust if the data have a smaller standard deviation. Observing Table 4, we can note that standard deviation values of GSA and FGSA are about the same, whereas DE results show less robustness than GSA and FGSA. Finally, the proposed approach shows a good robustness for the design of LP filters.  Table 5 shows the average pass and stop band ripples and transition band for HP filters. In this case, FGSA is better than DE. Moreover, FGSA has about the same value of stop band ripple of GSA, but pass band ripple and transition band are greater than GSA. However, FGSA gives a reasonable trend of magnitude response over frequency (see Figure 8). As in LP filters' results, the outcomes on 30 experiments show a good robustness of FGSA (see the standard deviation values in Table 5). For BP filters, FGSA and GSA are better than DE (see Table 6). The stop band ripple of FGSA is less than GSA, whereas pass band ripple and transition band are greater than GSA. Figure 9 shows a symmetric trend of GSA magnitude response. Moreover, good results of robustness are achieved (see Table 6).  Table 7). In particular, FGSA reduces the pass band ripple by 15% with respect to GSA. Moreover, Figure 10 shows that FGSA has a trend very close to an ideal SB filter. Table 7 shows low standard deviation values for FGSA: this fact assures a good robustness of the developed method. Usually, the search algorithm computation time is a measure of the procedure convergence speed. On the other hand, the computation time is not good to evaluate the convergence speed because it depends on hardware performances, programming language and designer skills. A good way to evaluate the convergence speed is to consider the objective function evaluations number up to the minimum value of the function (see [23,[31][32][33][34]53]). Finally, we consider the ratio between the fitness function calculations number n f and the evaluations number N: such ratio defines the convergence rate denoted by c r (see Equation (26)): Table 8 contains the convergence profile results of DE, GSA and FGSA for 8th order LP, HP, BP and SB IIR filters on 30 experiments. Each row contains the best value of fitness error function and standard deviation obtained by using the specified algorithm. In the table, f LP denotes the optimal value of fitness function for the filter LP, f HP denotes the optimal value of fitness function for the filter HP, f BP denotes the optimal value of fitness function for the filter BP and f SB denotes the optimal value of fitness function for the filter SB. Note that FGSA shows a better design performance for LP and SB filters, whereas there is a certain equilibrium between FGSA and GSA for HP and BP filters. Moreover, DE is the worst IIR filter design algorithm in terms of fitness function minimization. Referring to robustness, DE, GSA and FGSA have about the same low values of standard deviation of order 10 −3 (see Table 8). This fact assures a very good robustness of the proposed approach. Table 9 shows the convergence rates of DE, GSA and FGSA for 8th order LP, HP, BP and SB IIR filters. They are referred to as c r LP , c r HP , c r BP and c r SB for LP, HP, BP and SB filters, respectively. The analysis on the convergence rate results shows that FGSA has a convergence rate better than GSA for LP and HP filters, whereas GSA is better than FGSA for BP and SB filters. However, FGSA gives the best results because the improvements on convergence rate are 6% and 10% for LP and HP, respectively, versus 7% and 1% of GSA for BP and SB. Moreover, DE gives the worst values of convergence rate. Figure 11 shows the trend of the error fitness function over the number of iteration N for LP filters by using DE, GSA and FGSA. Note that FGSA has a better convergence rate than GSA and DE.
Finally, FGSA achieves the best compromise between IIR filter design performances and convergence rate with a good robustness. These facts make FGSA better than DE and GSA for the optimal design of 8th order IIR filters.
Stability analysis of the designed IIR filters is shown in Figures 11-15, where the circle markers represent the zeros, whereas the cross markers are the poles. The pole-zero plots demonstrate the existence of poles within the unit circle, which assures the Bounded Input Bounded Output (BIBO) stability condition. However, adding constraints to optimization algorithms may cause an increase of computational complexity. Recently, Pelusi et al. [54] have proposed a Neural and Fuzzy Gravitational Search Algorithm (NFGSA) able to search local optima with low complexity. The future challenge will be the application of NFGSA for designing optimal IIR filters with fuzzy stability constraints.

Conclusions
An intelligent algorithm able to optimize the design of 8th order IIR filters has been described. Because the quality of a filter depends on pass band ripple, stop band ripple and transition band, the target of the paper is the optimization of an error fitness function that depends on these parameters. Such task is accomplished through a suitable optimization algorithm. The proposed algorithm is a combination between fuzzy techniques and GSA. In particular, two fuzzy systems able to adjust some parameters of GSA have been designed. Moreover, to improve GSA, one of these parameters has been re-defined. Our algorithm has been compared with DE and GSA for the design of IIR filters. The results show that FGSA is the best algorithm to design 8th order IIR filters in terms goodness, robustness and convergence rate. Moreover, the proposed algorithm always gives a stable filter. Further research tasks will focus on: (1) the improvement of the fitness function definition; (2) the design of FISs for other GSA parameters assuring a good compromise between best solution and high convergence speed for the design of IIR filters; and (3) the comparison with other optimization algorithms such as Particle Swarm Algorithm and Genetic Algorithms. A future fascinating challenge will be the design of optimal IIR filters with fuzzy stability constraints.