IIR filter optimization using improved chaotic harmony search algorithm

Abstract Due to the fact that the error surface of adaptive infinite impulse response (IIR) systems is generally nonlinear and multimodal, conventional derivative-based techniques fail when used in adaptive Filter design. In this sense, global optimization techniques are required in order to avoid local minima. Harmony search (HS), a musical inspired metaheuristic, is a recently introduced population-based algorithm that has been successfully applied to global optimization problems. In the present paper, adaptive IIR filtering is formulated as a nonlinear optimization problem and then an improved version of HS incorporating chaotic search (CIHS) is introduced to solve the identification problem of three benchmark IIR systems. Furthermore, the performance of the proposed methodology is compared with HS and two well-known metaheuristic algorithms, genetic algorithm (GA) and particle swarm optimization (PSO) and a modified version of PSO called PSOW (Particle Swarm Optimization with weight Factor). The results demonstrate that the proposed method has superior performance over the other above-mentioned algorithms in terms of convergence speed and accuracy.


Introduction
Adaptive filtering techniques have been significantly advanced in recent years and have been successfully applied in a variety of fields in digital signal processing, communication and control. Adaptive filtering when used for system identification tends to provide a model that represents the best fit to an unknown plant based on finite impulse response (FIR) or infinite impulse response (IIR) structures. While theory and design of adaptive filters based on FIR filter structure is a well-developed subject, the same is not true for linear IIR systems. This stems from the fact that IIR structures tend to produce multimodal error surfaces with respect to filter coefficients. This fact leads conventional derivative-based learning algorithms such as least mean square easily gets stuck in local minima when solving such optimization problems. This is because they try to find the global minima by moving only in the direction of negative gradient. In this sense, each local minima becomes a potential trap that prevents algorithm from reaching the global extreme point. Classical recursive methods, in addition to convergence to local minimum, pose stability problem, slow convergence and also their performance substantially deteriorate when reducedorder adaptive models are used [1]. On the other hand, an IIR structure due to having both poles and zeros can give a better approximation of real-world systems. Besides, to achieve a particular level of performance, an IIR filter requires less number of coefficients than the FIR filter which corresponds to less computational burden. A number of classical adaptive system identification and filtering techniques have been reported in the literature [1][2][3][4]. Traditionally, least square techniques have been well studied for the identification of static and linear systems [5]. For nonlinear system identification, different algorithms have been used in the past including neural networks [6][7][8] and gradient-based search techniques such as least mean square [9]. In order to alleviate these deficiencies, population-based search algorithms such as genetic algorithm has also been used. But its effectiveness was affected by convergence time. After that, population-based stochastic optimization algorithms have been discussed in various literatures for design and identification of IIR filters which had the capability of faster convergence and global search of solution space in comparison to conventional methods [10][11][12][13]. Application of particle swarm optimization (PSO) and its variants could be found in [14][15][16]. Recently, a new approach based on artificial bee colony optimization for digital IIR system identification is proposed [17]. Seeker optimization, Cat swarm optimization and harmony search (HS) have been also proposed in [18][19][20].
When the complexity of the problem increases or where time allowed for convergence is limited (in dynamic systems), finding the global optimum becomes challenging. In these cases, hybrid algorithms are introduced to improve the performance by combining the best feature of each participating algorithms [14].
The goal of this paper is to introduce an enhanced HS algorithm by combining chaotic search and concepts from Swarm intelligence to avoid local minima and to enhance global convergence characteristics of the algorithm. Although, HS itself can produce good solutions at a reasonable time for a complex optimization problem, researchers are still trying to improve the fine-tuning characteristics and convergence rate of HS algorithm [21][22][23][24].
Chaos is one of the characteristics of nonlinear systems which include infinite non-periodic bounded motions. Nonlinear dynamic systems could be iteratively used to generate chaotic sequences of numbers. Many chaotic maps in the literature possess certainty, ergodicity and the stochastic property. As a novel optimization technique, chaos has gained much attention. For a given cost function, by following chaotic ergodic orbits, chaotic dynamic system may eventually reach the global optimum or its good approximation. Recently, chaotic sequences have been adopted in place of random sequences [25][26][27]. They also have been combined with some metaheuristic optimization algorithms to improve performance of these algorithms by chaotic evolution of variables [10,28,29]. This evolution includes two main steps: firstly, mapping from the chaotic space to solution space and then searching optimal regions using chaotic dynamics instead of random search [30].
In this paper in order to enhance the global convergence of HS algorithm, firstly we proposed a new variant of HS, called improved HS (IHS) by adopting concepts from swarm intelligence. Furthermore, sequences generated from logistic chaotic map substitutes random numbers for two key parameters of HS. Finally, to enhance the fine-tuning characteristic, top solutions found by HS are sent to a chaotic local search (CLS) based on the logistic map, where the best solution will be replaced if the result of CLS is better than HS. The simulation results pertaining to identification of three IIR systems and one nonlinear systems with reduced-order models show the superior performance of our proposed method compared to HS and two other well-known metaheuristic algorithms GA and PSO and a variant of PSO, called PSOW.
The remaining of this paper is organized as follows. Section 2 considers the mathematical formulation of IIR system identification. Section 3 discusses the HS algorithm. The improved HS with chaos is introduced in Section 4. The proposed CIHS-based IIR system identification method is provided in Section 5. The simulation results and discussions are given in Section 6. The paper ends with conclusions in Section 7.

IIR filtering and system identification
The application of IIR filter in system identification has been widely studied since many problems in signal processing can be characterized as a system identification problem ( Figure 1) [12,18]. The problem of determining a mathematical model for an unknown system by monitoring its input-output data is known as system identification [2]. The task of any given parametric system identification algorithm is to vary the model parameters until a predefined approximation criterion is satisfied. The block diagram of an arbitrary IIR system identification algorithm is shown in Figure 1. The adaptive algorithm essays to tune the adaptive filter coefficients such that the error between the output of the unknown system and the estimated output is minimized. In other words, the parameters of IIR filters are successively updated by the algorithm to solve a minimization-optimization problem.
An IIR system is described as: wherein, is the IIR transfer function, Y(z) is the z-transform of the output y(n), and U(z) denotes the z-transform of the input u(n). a i, i = 0, 1, 2, . . . , mand b i , i = 1, 2, 3, . . . , n are the feed-forward and feed-back coefficient of the IIR system, respectively. The IIR filter can be formulated as a difference equation As illustrated in Figure 1, v(n) is the additive noise in the output of the system. Combining v(n) and y(n), we get the overall output of the system d(n). Additionally, Figure 1. Block diagram of system identification using an IIR filter.
for the same set of inputs, the adaptive IIR filter block producesŷ(n). The estimated transfer function can be represented aŝ whereâ i andb i signify the approximated coefficients of the IIR system. In other words, the transfer function of the actual system is to be identified using the transfer function of the adaptive filter. The difference between d(n) andŷ(n) produces the input to the adaptive algorithm. The adaptive algorithm uses this residual to adjust the parameters of the IIR system. It can be concluded from the figure that The cost function (mean square error) to be minimized by the adaptive identification algorithm is given by where, N denotes the number of input samples and E(.) is the statistical expectation operator. The optimization algorithms employed in this paper search the solution space to locate those values of parameters, which contribute to the minimization of (7).

HS optimization
The HS metaheuristic is a novel optimization algorithm inspired by the underlying principles of music improvization that is successfully used in several sciences and engineering applications [31][32][33][34]. When musicians are improvizing, they usually test various pitch combinations to make up a harmony. In fact, the aim of the music is to search for a perfect state of harmony. In this sense, the process of searching for optimal solution in engineering is analogous to this search for a pleasing harmony in memory. In a real optimization problem, each musician is replaced by a decision variable and favourite pitches are equivalent to favourite variable values. Table 1 presents a comparison between music improvization and optimization. In order to explain the HS algorithm in more detail, it is required to idealize the improvization process done by an expert musician. There are three possible choices for a musician, (1) to play any famous pitch from memory. (2) to execute a pitch adjacent to any other in his memory (3) execute a random pitch from the range of all possible pitches [32]. Geem et al. [35] in 2010 formulized these three options to create a new metaheuristic. The corresponding components of these options  [32]. The last component is the randomization, which is to increase the diversity of the solutions. Although pitch adjusting has a similar rule, it is limited to a local search. The use of randomization can provide exploration of various regions so as to find the optimum through global search. The steps of the HS algorithm can be followed in Figure 2.

Improved HS with chaos
In this section, the IHS algorithm will be introduced first, and then, after a brief description on chaotic behaviour, the proposed chaotic IHS algorithm will be discussed.

Improved HS
Compared to the other metaheuristics like GA or PSO, HS offers some advantageous: It has fewer mathematical requirements and thus easier implementation. In addition, there is evidence to suggest that HS is less sensitive to its parameters than PSO [36]. Although the basic HS is efficient but improvable, it can be seen from the simulation results that the solutions are still changing as the optima are approaching [20,37,38]. In order to improve the convergence of the HS, some works have been performed in [21][22][23][24]. In this paper, an IHS will be proposed which is inspired by the concept of global best in PSO. We will modify the HM consideration step such that the new harmony can mimic the best harmony in HM. When using HM, it is desirable to pick up the best harmony from the memory. This will happen with a probability of 1/HMS in the basic HS since a random selection is applied to choose the new harmony from the HM. Our proposed HIS has exactly the same steps as the HS with the exception that HM consideration step is replaced by HM (gbest) instead of HM (fix(U(0,1) *HMS). Where gbest refers to the best harmony among all harmonies in terms of minimum fitness function.

Chaos
Chaos is a deterministic, pseudorandom dynamic behaviour in nonlinear dynamical systems that are non-periodic, non-converging & bounded. It exhibits sensitivity dependence on initial conditions. Mathematically, chaos is randomness of a simple deterministic dynamical system and chaotic system may be considered as source of randomness [22] and [23].
Although, it appears to be stochastic, it occurs in a deterministic nonlinear system under deterministic conditions. A chaotic map is a discrete-time dynamical system running in the chaotic state.
The chaotic sequence z k : k = 0, 1, 2, . . . can be used as spread-spectrum sequence and as a random number sequence. The logistic map, circle map, and sinusoidal map are among the well-known one-dimensional maps used in the chaotic search. The logistic map is represented by the following equation [39].
In recent years, chaos has been extended to various optimization areas like in [40,41]. In random search optimization algorithms, the method using chaotic variable instead of random variable are called chaotic optimization algorithms. In these algorithms due to the non-repetition and ergodicity of chaos, it can carry out global searches at higher speed than stochastic searches that depends on probabilities [42]. Many have used chaos for multimodal optimization problems in the past. The effective role of chaotic search in finding global minima are reported in [40][41][42][43].

Chaotic harmony search
It was mentioned in Section 3 that one of the drawbacks of the HS is its premature convergence, especially while handling with more than one local optima. HCMR, PAR, FW and the initialization of HM are key factors to affect the convergence of HS. In classical HS, the abovementioned parameters are adjusted as fixed values in initialization step. In this method, the number of iterations plays an important role to find an optimal solution. For example, for a small PAR and large FW, much number of iterations are required to find optimum solutions. Small FW values in final iterations increase the fine-tuning of solution vectors by local exploitation and in early iterations bigger FW value can increase the diversity of solutions vector for global exploration [43]. In [21,23], some efforts have been done to dynamically update the PAR and FW. These parameters could be selected chaotically by using chaotic maps. In 2010, Alatas proposed different chaotic HS algorithms by applying different chaotic maps to the HMCR, PAR, FW and the initialization of HM [43]. In this paper, PAR and FW values have not been fixed and they have been modified as follows to improve the global convergence by escaping the local optima.
PAR (t + 1) = f (PAR (t)) , 0 < PAR (t) < 1, t = 0, 1, 2, . . . (11) where f corresponds to logistic map and t denotes iteration number. Furthermore, the algorithm is hybridized with a CLS. For this purpose, when the three steps of HS algorithm are performed, the top HMS /5 of harmonies in the HM are passed through a chaotic map (logistic map in this paper). This step is performed by (1) mapping the top decision variables in solution space to chaotic variable in the interval (0,1), (2) determining the new chaotic variable for each of them using the logistic equation and (3) converting new chaotic variables to decision variables in the range of the solution space. Then the best harmony of this chaotic search is chosen and compared to the new harmony resulted from the main algorithm and the new solution will be updated if the result of CLS is better.

CIHS-based IIR system identification
The identification algorithm can be summarized in the following steps: where, Y = [y(1) y (2) . . . y(p)] T , is the output of the plant contaminated by measurement noise, and that corresponds to the least fitness (best attainable match between the IIR model and the actual system in the sense of MSE) shows the estimated parameters.

Results and discussions
In this section, three benchmark linear IIR systems are considered for the case study. Parameter identification burden is carried out using a model having less order than of the actual system. These reduced-order cases pose challenge to the optimization algorithm since they produce highly multimodal error surfaces. In addition, as the number of coefficients decreases, the degree of freedom reduces and it becomes more difficult to identify the actual system. In order to ensure the validity of the results, each experiment is repeated in 20 consecutive trials and the resultant (best, worst, standard deviation, and mean) values of the minimum MSE's of each run are given in corresponding tables. To provide a more comprehensive comparison, other than the proposed algorithm (CIHS), same simulations are repeated using standard versions of GA, PSO and HS and a modified version of PSO introduced as PSOW, as well. Each simulation is carried out in MATLAB v.7.5. In all cases, the population size is set to 50, the maximum number of iterations (NI) is set to 1000 and the input data is a Gaussian white noise with zero mean and unit variance. The output data is contaminated with a Gaussian random noise with zero mean and a variance of 0.001.Parameter adjustment is as the following: In CIHS PAR and FW values are determined by a logistic map in each iteration and HMCR = 0.95. In basic HS parameters are set as PAR = 0.5, HMCR = 0.95, and FW is bounded to 1% of each variable range. In GA algorithm, the bit number is set to 16, mutation probability is 0.1, and crossover step is of single-point type with a probability of 0.6. In PSO, acceleration constants are set to 2 and inertia weight is linearly decreased from 0.9 to 0.4. In PSOW, the adaptive inertia weight factor w, is determined as follows: where w min and w max denote the maximum and minimum w, respectively. f is the current objective function value of the particle, f avg and f min are the average and minimum objective values of all particles, respectively [30]. In this approach, w is varied based on the objective function value so that particles with low objective values can be protected while particles with objective value greater than average will be disrupted. Hence, it provides a good way to maintain population diversity and sustain good convergence capacity. The constant parameters for PSOW are: w max = 1.2, w min = 0.2, c 1 = c 2 = 2 and v max is limited to the 15% of the search   space. In all algorithms, random numbers take values between 0 and 1.
Example 9.1: Consider the following IIR system System (14) is modelled using the following IIR structure The simulation results related to the reduced-order model (15) are given in Figure 3. The analysis of the figure shows that the utilization of CIHS has resulted in greater estimation accuracy and higher convergence speed. HS converges to its minimum MSE level after about 700 numbers of iterations, which demonstrate its poor convergence speed. Table 2 shows that the CIHS provides the best average result in terms of the MSE. Table 3 demonstrates that HS requires much less computational time than the other algorithms. CIHS needs higher computational time than HS, mainly because it has to go through more number of fitness evaluations. Each iteration of HS corresponds to q number of fitness evaluations while that of CIHS corresponds to q + HMS /5.

Example 9.2:
The transfer function of the plant is given by the following equation This third-order plant is modelled using a second-order IIR filter. Hence the transfer function of the model is given by The convergence characteristic shown in Figure 4 demonstrates that GA, PSO and PSOW converge to a suboptimal solution. However, HS and CHS do not stagnate and reach their minimum noise floor level. The CHS takes 250 generations to reach its minimum MSE level, whereas HS takes 700 generations to reach its corresponding value. CHS outperforms the HS with higher convergence speed and lower MSE value. Table 4 shows that the CHS provides the best average result in terms of MSE. Table 5 indicates that CHS needs higher computation time than the HS. Hence, there is a tradeoff between the quality of solution and computational time.
The IIR filter transfer function used for the identification purpose is given below: Using Equation (19), the same set of simulations have been executed as in Examples 9.1 and 9.2. Figure 5 represents the average MSE graphs from 20 simulation tests. It is shown in Figure 5 and Table 6 that the convergence speed of CIHS is much greater than HS and the minimum MSE level obtained using CIHS is much smaller than that of HS, PSO, PSOW and GA. The CHS takes 250 generations to reach its minimum MSE level, whereas HS solutions are still changing as the number of generations reaches 1000. The respective computation time of the algorithms are listed in Table 7.
The computational time of an adaptive filtering algorithm is a critical issue in real-time applications.   The results demonstrate that HS-based adaptive IIR filtering algorithms are much faster than the GA and PSO. Altogether, the simulation results given in this section reveals that CIHS has minor chance of premature convergence and hence, it is a promising optimization tool in IIR adaptive filtering. According to , HS solutions tends to reach lower errors as the number of iterations reaches 1000. To investigate the impact of higher iteration numbers on convergence characteristics of HS-based algorithms, we repeated the experiment of Example 9.3 with HS and CIHS algorithms. Equation (18) is used as the transfer function of the plant and Equation (19) as the IIR filter. All the parameters of algorithms are the same as the parameters of Example 9.3 except the iteration number which is set to 10,000 this time. Figure 6 shows the results. Both algorithms have reached better results in terms of accuracy but CIHS again showed better accuracy and speed.

Conclusion
In this paper, a new version of Harmony Search algorithm (CIHS) is developed using social component of PSO and chaotic search combined to the original algorithm to enhance exploration and exploitation capability during the search. The new algorithm is outlined and has been applied to identification of three benchmark IIR plants. The performance assessment of the CIHS in identifying an unknown system with a reduced-order IIR model in comparison to those obtained by the HS, PSO, PSOW and GA clearly exhibits faster convergence and lower values of MSE for CIHS which makes it the best algorithm among the five, for adaptive system identification. In addition, it has been shown that HS-based IIR system identification methods would result in a much less computational complexity. Therefore, the proposed method can be employed in real-time tasks. To confirm the robustness of the proposed algorithm, CIHS needs to be applied to more complex and real-world optimization applications. For future work of the authors, investigating the result of various chaotic maps on the algorithm's performance and real-time hardware implementation of the algorithm will be considered. Also, comparison of the performance between the proposed algorithm with other recently introduced nature-inspired algorithms can be a subject of future works.