Digital IIR Filters Design Using Differential Evolution Algorithm with a Controllable Probabilistic Population Size

Design of a digital infinite-impulse-response (IIR) filter is the process of synthesizing and implementing a recursive filter network so that a set of prescribed excitations results a set of desired responses. However, the error surface of IIR filters is usually non-linear and multi-modal. In order to find the global minimum indeed, an improved differential evolution (DE) is proposed for digital IIR filter design in this paper. The suggested algorithm is a kind of DE variants with a controllable probabilistic (CPDE) population size. It considers the convergence speed and the computational cost simultaneously by nonperiodic partial increasing or declining individuals according to fitness diversities. In addition, we discuss as well some important aspects for IIR filter design, such as the cost function value, the influence of (noise) perturbations, the convergence rate and successful percentage, the parameter measurement, etc. As to the simulation result, it shows that the presented algorithm is viable and comparable. Compared with six existing State-of-the-Art algorithms-based digital IIR filter design methods obtained by numerical experiments, CPDE is relatively more promising and competitive.


Introduction
Filtering problem [1,2] is a widely studied research topic in various fields of control and signal processing. The main objective of filtering is synthesizing and implementing a filter network [3] to modify, reshape, or manipulate the frequency spectrum of a signal according to some desired specifications. As one of the most successful filter networks, the well-known Digital infinite-impulseresponse (IIR) filter has been extensively used in many practical systems [4][5][6][7], such as engineering system, network system, nuclear reactor, biological system, chemical system and electrical networks system. However, it has been recognized now that the IIR filter will generally not guarantee satisfactory performance if its feedback coefficients are chosen not appropriately during the adaptation process [8]. Apart from this disadvantage, the possibility of having a multi-modal and nonlinear error surface is another important design challenge for recursive filters [9,10]. To improve the robustness, in recent years, many heuristic optimization design methods have been developed, such as simulated annealing (SA) [11], ant colony optimization (ACO) [12], particle swarm optimization (PSO) [13], seeker optimization algorithm (SOA) [14], artificial bee colony (ABC) [15,16] and differential evolution (DE) [17], etc.
A SA is usually sensitive to its starting point of the search and requires too many function evaluations to converge to the global minima. The ACO imitates the social behavior of real ant colonies and it has been originally developed for combinatorial optimiza-tion problems. But, it may occasionally be trapped into local stagnation or premature convergence resulting in a low optimizing precision or even failure [18]. What's more, the conventional PSO algorithm [19] as shown in several studies can easily fly into the local optima. It also lacks the ability to jump out of the local optima when solving complex multimodal tasks. The SOA simulates the act of human searching and has been widely developed for system identification [20]. Nonetheless, the performance of SOA is also affected by its parameters, and it could not easily escape from premature convergence. Differential evolution, proposed by Storn and Price [21], is a population-based heuristic search algorithm with dual features of reliability and flexibility. It implements the evolutionary generation-and-test paradigm for global optimization by using the current population information of distance and direction to guide the search. It has many advantages such as simplicity, reliability, high performance and easy implementation, which gives great potential application to IIR design. In the seminal DE algorithm, perturbation is operated by adding a weighted moving vector (the weight F is called scale factor) and modifying the values of some randomly selected coordinates. The perturbed solution, namely the offspring, is then evaluated by means of its objective function and compared with its corresponding parent. If the newly generated solution outperforms its parent, then a replacement occurs; otherwise the parent solution is retained. To provide a rigorous proof for its probabilistic convergence, [22] has modeled the population as a dynamical system in which the probability density function (PDF) of the population vectors changes with time. It was shown therein that the dynamics is asymptotically stable (which implies convergence) at the equilibrium PDF, which is a Dirac delta function placed at the global optima. Later on, various mutation strategies [23] were used for the generation of new solutions to augment the robustness of the underlying algorithm.
In DE, it is often the case that, for optimization problems such as single-objective, multi-objective, large scale, constrained and dynamic problems, the population size is naturally fixed on a constant value; see, e.g., [24]. Unfortunately, it is usually difficult to determine how large the population size is suitable for solving numerical optimization problems. For instance, a definite population size is given in [25] which increases linearly with the problem dimension; yet the sparse and noisy data makes it difficult to accurately estimate the maximum population size. Inspired by this fact, an efficient population utilization strategy for DE (DESAP) [26] was developed to automatically tune its population size from initialization to completion right through the evolutionary search process. Nevertheless, the population utilization method depends on its encoding methodology, which is a restriction for the current population with complex dynamical behaviors [27]. No significant advantages can be observed while using relative encoding. Subsequently, the idea of population adaptation has been applied in solving dynamic optimization problems [28], where multi-population approach (DynDE) is placed onto DE aimed at locating optima faster. Yet such a method requires a determinated topology that may be sensitive to the noise of measurement in some extent [29]. It is worth mentioning that although the population may not be as large as possible, it ought to meet the requirements of given engineering. Therefore, a new reduction method for the population size was shown up for the jDE in order to enhance algorithmic performance [30], where the population size was progressively declining until the arrival of the final budget during the optimization process. Unfortunately, this method can not keep track of the progress of individuals in the sustainable reduction.
Many studies have indicated that various computational predictors or models developed in biology and biomedicine, such as those in identifying DNA-binding proteins [31], predicting Gprotein-coupled receptors (GPCRs) and their types [32], identifying nuclear receptors and their subfamilies [33,34], identifying the subcellular localization of proteins from various organisms [35][36][37][38][39], can timely provide very useful insights and informations for both basic research and drug development. These predictors all use the methods of digital signal processing. In view of this, the present study is attempted to addresses an important problem in designing digital IIR filters in hopes that it may become a useful tool for the related information-treating areas. However, most of the developed adaptive population methods have their advantages and disadvantages. So far, it remains open that how to utilize the dynamic population strategy to solve real-world practical problem. We aim at employing a Markov jumping (switching) population updating DE for digital IIR filter, so that the dynamic population can quickly converge to the potential global optimum by taking advantage of the current search information. Thus, the CPDEbased evolutionary method is simulated in digital IIR filter design, and its performance is compared to that of three versions of DE, CMA-ES, GL-25 and SOA. In the community of six digital IIR filter design problems, it is shown empirically that CPDE is capable of producing highly competitive results compared with other EAs.

Illustration
Application of the IIR filter in system identification has been widely studied since many problems encountered in signal processing can be characterized as a system identification problem ( Figure 1) [40,41]. Therefore, in the simulation study, IIR filters are designed for the system identification purpose. In this section, we will utilize a modified DE to adjust the parameters of the filters until the error between the output of the filter and the unknown system is minimized. Subsequently, we provide an overall comparison between the performance of CPDE and several other State-of-the-Art algorithms to verify the effectiveness and usefulness of the proposed method.
Six typical system identification problems [14] make up the test suite used for this comparative study, which are listed in Tables 1  and 2. H s (z) and H f (z) specify the system and filter transfer functions, respectively; x(k) indicates the system input; SNR is the Signal to Noise Ratio; v(k) is the system noise, which is independent of x(k); N (0,1) presents the white Gaussian noise (WGN) in zero-mean normal distribution with variance 1. In Table 2, w records all coefficients of six digital IIR filters; Search space is the predefined boundary constraints, that will be analyzed in Section 05; N denotes the data length used in calculating the mean-square-error (MSE). The examples were selected so as to include problems with the following characteristics: unimodal/ multi-modal, no noise/noisy. For each algorithm and each test function, 30 independent runs are conducted with 100,000 FES as the termination criterion.
Traditionally,''generation'' is a natural form of computational cost for statistical comparison [11][12][13][14]17]. However, the population may not be the same in different algorithms. The algorithm with a larger population may obtain a better performance together with much more function evaluations in every generation. Thus, in this paper, the function evaluations (FES) are conducted here to represent its computational cost for algorithm comparison.
In all simulations, the population size of the most EAs is 100 with the exception of EPSDE and SaDE. As suggested in Ref. [42,43], the population size of EPSDE and SaDE is chosen to be 50. Seven existing EA algorithms are shown in Table 3 in detail. CMA-ES represents the state of the art of Evolution Strategies and it is a referent in the continuous optimization field. GL-25 is a hybrid real-coded genetic algorithm which combines the global and local search. EPSDE is an adaptive DE with ensemble of parameters which incorporates a self-organizing method. jDE is a standard DE with adapted parameter setting. SaDE delivers a mutation strategy pool where strategy is self-adapted based on their previous performance. SOA is a novel heuristic stochastic optimization algorithm based on the simulation of the act of human searching. The parameters for these EAs are provided in Table 3.

Comparison on the Solution Accuracy
In this section, an overall comparison of the performance is provided between the CPDE variant and other six State-of-the-Art EAs (i.e., CMA-ES, GL-25, EPSDE, jDE, SaDE and SOA). We evaluate the performance of seven heuristic algorithms over six typical nonlinear uncertain discrete-time problems. Fig. 2 illustrates the cost function value versus number of evaluations averaged over 30 random runs for the seven algorithms. The subfigures amplify the convergence graphs in clarity. Table 4 reports the experimental results of Examples 1-6, averaged over 30 independent runs with 100,000 FES.
From the Table 4 and Fig. 2, the CPDE provides the best performance on the ex 2 {ex 3 , ex 5 and ex 6 , then ranks the second on the ex 1 and ex 4 . SOA gives the best performance on the ex 1 and ex 4 . The results show that GL-25 and SOA have good ability of convergence speed. Fig. 3 also shows instance of evolution of the parameters of two filters for CPDE.
To be specific, on the three multimodal problems (ex 1 {ex 2 and ex 6 ), CPDE performs much better than CMA-ES, GL-25, EPSDE, jDE, SaDE and SOA on the latter two functions. SOA delivers good accuracy and the highest convergence rate on ex 1 , while CPDE outperforms other five methods. To sum up, CPDE is the winner on multimodal functions. This might be due to the fact that CPDE implements the overall adaptive variable population size method, which can help the DE search the optimum as well as maintain a higher convergence speed when dealing with multimodal rotated functions. On the remaining three unimodal functions (ex 3 {ex 5 ), CPDE performs significantly better than six others for ex 3 and ex 5 . SOA can provide good accuracy on ex 4 , while CPDE achieves the highest convergence rate. The outstanding performance of CPDE is due to its dynamic PDS, which leads to very fast convergence. Overall, the CPDE is the best among the seven methods in the comparison conducted on unimodal functions and expanded multi-modal functions. For a thorough comparison, the t-test has also been carried out in this paper. Table 4 presents the total score on every function of this two-tailed test with a significance level of 0.05 between the CPDE variant and other heuristic algorithms. Rows ''+ (Better),'' '' = (Same),'' and ''-(Worse)'' give the number of functions that the CPDE performs significantly better than, almost the same as, and significantly worse than the compared algorithm on fitness values in 30 runs, respectively. As confirmed in t-test, the CPDE in general offers more improved performance than the other six State-of-the-Art EAs.
a unit-variance white Gaussian pseudonoise sequence 0 and As mentioned above, the CPDE has shown a very competitive performance in the six filtering problems. In practical engineering, noise exist universally in nature [44]. Therefore, in the past few decades researchers have witnessed significant progress on filtering and control for linear/nonlinear systems with various types of noises among which Gaussian noise is one of the most general signals that has been widely studied [45,46]. Here, we further evaluate the proposed framework on the six expanded stochastic systems, where a zero-mean Gaussian white-noise is added. The maximum number of FES is set to be 100,000 in all runs. Table 5 summarizes the experimental results.
From the Table 5, the CPDE provides the best performance on the ex 2 , ex 3 {ex 6 , then ranks the second on the ex 1 . SOA offers the best performance on the ex 1 . The results show that GL-25 and jDE have a good ability of convergence speed.
More specifically, on the three multimodal problems (ex 1 {ex 2 and ex 6 ), although it worked slightly weaker on some functions, the CPDE in general offered more improved performance than all the EAs compared. It performs much better on the ex 2 and ex 6 , and attains slightly worse performances than the best solutions on the ex 1 . On the remaining three unimodal functions (ex 3 {ex 5 ), CPDE performs much better than other EAs. Hence, CPDE exhibits the highest performance in noise-expanded filtering problems, which can efficiently adjust the population structure and guide the evolution process toward more promising solutions.
The t-test is also summarized in Table 5. In fact, CPDE performs better than CMA-ES, GL-25, EPSDE, jDE, SaDE and SOA on 6, 6, 6, 6, 6 and 5 out of 6 test functions, respectively. Thus, CPDE is better than other six competitors in filtering system identification problems. Comparing the results and the cost function graphs, among these seven EA algorithms, the GL-25 can converge to the best solution found so far very quickly though it is easy to stuck in the local optima. The SOA has good global search ability and slow convergence speed. The jDE have good search capability on noise-expanded filtering problems. The CPDE has good local search ability and global search ability at the same time.

Comparison on Convergence Rate and Successful Percentage
The convergence rate for achieving the global optimum is another key point for testing the algorithm performance. Note that in solving real-world optimization problems, the ''function evaluation'' overwhelms the algorithm overhead [47]. Hence, the computation times of these algorithms are not provided here. Table 6 shows that CPDE needs the least FES to achieve the acceptable solution on ex 3 {ex 4 and ex 6 , which reveals that CPDE has a higher convergence rate than other algorithms. Though SOA or GL-25 might outperform CPDE on the other functions, SOA and GL-25 have much worse successful ratio and accuracy than CPDE on the problems for comparison. In addition, CPDE can achieve accepted value with a good convergence speed and accuracy on most of the problems, as shown in Tables 4, 5 and Fig. 2. Tables 4, 5 also show that CPDE yields a highest percentage for achieving acceptable solutions in 30 runs. According to the no free lunch theorem [48], any elevated performance over one class of problems is offset by performance over another class. Hence, one algorithm cannot perform better on convergence speed and accuracy than the others on every optimization problem.
In summary, the CPDE performs best on unimodal problems with or without noise and has a good search ability of multimodal problems. Owing to the controllable probabilistic technique, the CPDE processes capabilities of fast convergence speed, highest successful ratio and the best search accuracy among these EAs.

Performance of Controllable Probabilistic Approach
In this section, the controllable probabilistic (CP) approach is used to test the search performance of CPDE. In all the experiments, threshold f is adjusted in the following. Moreover, the parameter p 2 is also considered, which denotes the number of potential candidates for perturbation.
In this paper, f indicates the trigger thresholds, which is used to control the sensitivity of the dynamic CPDE. While f is set as one, the population size will be adjusted in each generation. Setting a higher f value will result in a lower sensitivity of the CPDE, while a lower f value will lead to a higher efficiency of the population adjustment. Notice that the parameter f should set to be larger than one. Failure to do this will result in an instant elimination of a newborn individual with poor performance, which may provide some degree of diversity preservation. On the other hand, coefficient p 2 also influences the perturbation process substantially. Table 7 shows the comparisons between CPDE with other three parameter settings of CPDE over Examples 1-6 with noise perturbation. It indicates that CPDE is not sensitive to the adjustment of parameters. In order to make a balance of the search accuracy and robustness, f~4 and p 2~5 are used as a representative parameter setting in our paper. This setting will prevent the instant elimination of a newborn individual and keep the CP approach high sensitivity.

Discussion
The CPDE is an improved differential evolution algorithm with a controllable probabilistic population size. When particles are clustered together in a region and trapped into the local basin, CPDE perturbs the population and generates the necessary ''fine'' individuals to share their up-to-date information. In addition, CPDE removes redundant individual with its entropy and ranking metrics to save computational load. In this paper, a CPDE-based digital filter design method has been proposed, and the benefits of CPDE for designing digital IIR filters have been studied.
An overall comparison between the performance of the CPDE and other six State-of-the-Art EAs (i.e., jDE, SaDE, EPSDE, CMA-ES, GL-25 and SOA) was provided over 6 typical robust  system identification problems, and the result clearly indicated the CPDE achieved a substantially significant improvement on the performance. Furthermore, convergence rate was also validated that the CPDE has good convergence performance to achieve the fixed accuracy level with acceptable generations. Thus, it is believed that the proposed CPDE is capable of rapidly and efficiently adapting the parameters of a wide variety of IIR structures and will become a promising candidate for digital filter design. Previous work has shown the importance of system identification on digital IIR filter design. Furthermore, CPDE has effectively been applied to estimate the structure of nonlinear uncertain discrete-time system. Therefore, our method is possible to be used to reconstruct the topology structure for on-line  adaptive filtering applications [49]. Another possible application is to identify topology and parameters of complex networks [50][51][52] and biological time series [53] by dynamic population strategy, provided online measurement and increasing/decreasing techniques are feasible. Generally, the suggested technique enables us to identify the unknown parameter of real networks which allows the required control applications (perturbations). Some possible experimental research is now under our investigation in controller design and tuning.

Description of the Problem
Consider the digital IIR filter with the input-output relationship governed by the difference equation: where x(k) and y(k) are the filter's input and output, respectively, and M( §L) is the filter order. The transfer function of this IIR filter can be written in the following general form: Hence, the design of this filter can be formulated as an optimization problem with the cost function J(w) stated as follows: where d(k) and y(k) are the desired and actual responses of the filter, respectively; and e(k)~d(k){y(k) is the filter's error signal; N is the number of samples used to calculate the objective function.
w~½a T b T T~½ a 0 a 1 :::a L b 1 :: denotes the filter coefficient vector. The aim is to minimize the cost function J(w) by adjusting w. An important consideration during the adaptive process is to maintain the stability of the IIR filter. Not all filters defined by Eq. (1) are feasible or implementable.  An efficient way of maintaining stability is to convert the direct form to a lattice form and make sure that all reflection coefficients k i , 0ƒiƒM{1, have magnitudes less than one. We will adopt a similar approach as in [11] to guarantee the stability of the IIR filter during adaptation. Thus, the actual filter coefficient vector used in optimization is w~½a 0 a 1 :::a L k 0 :::k M{1 T . In the circumstances, the coefficient space Y is formed by the constraints of a i [½{2,2(i~1,:::,L) and the magnitudes of k i that are less than one. For the sake of simplicity, we adopt the predefined boundary constraints as ½{2,2 to compare other existing EAs fairly.

A Controllable Probabilistic DE
The stochastic system is iterated forward in time using a synchronous DE updating scheme. In our work, a modedependent population updating equation with Markovian switching parameters is introduced with the hope to keep track of the progress of individuals and further improve the search abilities. A detailed algorithm design of CPDE can be found in Documentation S1.
In CPDE, there are two levels of sub-optimizers, population decreasing strategy (PDS) and population increasing strategy (PIS). The probability of selecting different sub-optimizer to improve the online solution-searching status is completely up to its nonhomogeneous Markov chain. For choosing required sub-optimizers adaptively, consider the following probability transition matrix in Eq. (5): Here, M(1)~a, M(2)~b and M(3)~c stand for the population maintaining state, population increasing state and population decreasing state, respectively. The expectations of Markov chain P are automatically updated by the search environment. It is worthwhile to mention that p 22 and p 33 are set to be 0. In this case, the increasing and decreasing operators will not be performed in consecutive generations.
If few trial vectors can outperform the corresponding parent in selection operation, particles may be clustered together and trapped into the local basin. In such a case, the PIS is employed to add new individuals into the population and share their up-to-date information to help the individuals escape the local basin.

5, ð6Þ
where d 2 is the number of dimensions for perturbation, which is monotonically decreasing through the evolutionary search. Parameters a and b are the magnification coefficient. Parameters l and G are here considered as the generation variables. During early stages of the optimization process, much more reproductions will be generated to spread out its particles within the decision space. Nevertheless, solutions of the population tend to concentrate in specific parts of the decision space during the later period of optimization process. However, in case most of individuals can spawn new promising offspings in the evolutionary process, it then signifies that redundant intermediate particles exist. In this case, we introduce the PDS to remove poor particles to avoid undesirable computational cost and excessive search complexity.
where t denots an overall deletion indicator. The variable rank i and H i (D) indicate the rank metric and the entropy metric for individual i, respectively. It can be observed that from Eq. (7), for the individuals that have high rank values (i.e., away from the global best solution) or low entropy values (i.e., located in the crowded regions); these particles will have a higher probability of elimination.

Supporting Information
Documentation S1 A material algorithm design of CPDE.

(PDF)
Author Contributions