Particle Swarm Optimization of Highly Selective Digital Filters over the Finite-Precision Multiplier Coefficient Space

Digital filters find wide variety of applications in modern di gital signal processing systems [ 1, 2]. As a result of the recent progress in such systems, there is an ev er growing demand for sharp transition band digital filters. These narrow transition bandwidth dig ital filters are usually designed by using the frequency response masking (FRM) approach [ 3]. The computational efficiency of the FRM technique makes it suitable for different applications, e.g. in audio signal processing and data compression [ 4].


Introduction
Digital filters find wide variety of applications in modern digital signal processing systems [1,2]. As a result of the recent progress in such systems, there is an ever growing demand for sharp transition band digital filters. These narrow transition bandwidth digital filters are usually designed by using the frequency response masking (FRM) approach [3]. The computational efficiency of the FRM technique makes it suitable for different applications, e.g. in audio signal processing and data compression [4].
Practical design of digital filters is based on optimization for satisfying the given design specifications together with the hardware architecture. However, the optimization may be carried out in terms of fixed configurations but variable multiplier coefficient values. On the other hand, the problem may concern the optimization of the hardware architecture without taking the multiplier coefficient values into consideration.
In order to optimize the given design specifications, the multiplier coefficient values can be determined in infinite precision by using hitherto optimization techniques. However, in an actual hardware implementation of the digital filters, the infinite precision multipliers should be quantized to their finite precision counterparts, but these finite precision multiplier coefficients may no longer satisfy the given design specifications. Consequently, from a hardware implementation point of view, there is a need for finite precision optimization techniques, capable of finding the optimized digital filter rapidly while keeping the computational complexity at a desired level. In principle, there exist two different techniques for the optimization of digital filters, namely, gradient-based and heuristic optimization approaches.
space. In [6], a Remez exchange algorithm was used for the optimization of FRM finite impulse response (FIR) digital filters and it was shown that this algorithm may provide a speed advantage over the linear programming approach. However, both these techniques suffer from sub-optimality problems. In [7], unconstrained weighted least-squares criterion was used to develop another technique for the optimization of digital filters. Convex optimization approaches such as semi-definite programming [8] and second-order cone programming [9] have also been applied to the optimization of digital filters. However, if a large number of constraints are present, these optimization techniques may become computationally inefficient in terms of time consumption and speed.
Heuristic optimization algorithms have emerged as promising candidates for the design and discrete optimization of digital filters, particularly due to the fact that they are capable of automatically finding near-optimum solutions while keeping the computational complexity of the algorithm at moderate levels. Simulated annealing (SA) and genetic algorithms (GAs) were widely used in the design and optimization of digital filters [10][11][12]. Particle swarm optimization (PSO) and seeker optimization algorithm (SOA) are two newly developed algorithms suitable for the optimization of various digital filters due to their few number of implementation parameters and high speed of convergence [13,14]. It was shown that SOA has advantages over PSO in terms of the speed of convergence and global search ability [15]. Tabu search (TS) [16], ant colony optimization (ACO) [17], immune algorithm (IA) [18] and differential evolution (DE) [19,20] are alternative candidates for the optimization of digital filters. All the foregoing techniques allow a robust search of the solution space through a parallel search in all directions without any recourse to gradient information. However, the aforementioned techniques were developed for infinite precision optimization of digital filters which require the user to perform a quantization step for a hardware implementation.
In [21][22][23], a technique was developed for finite-precision design and optimization of FRM digital filters using GAs. finite-precision optimization of FRM FIR digital filters using PSO was studied in [24,25] and finite-precision optimization of infinite impulse response based (IIR-based) FRM digital filters was studied in [26,27]. PSO was originally proposed by Kennedy and Eberhart in 1995 as a new intelligent optimization algorithm which simulates the migration and aggregation of a flock of birds seeking food [28]. It adopts a strategy based on particle swarm and parallel global random search, that may exhibit superior performance to other intelligent algorithms in computational speed and memory. In PSO, a potential candidate solution is represented as a particle in a multidimensional search space, where each dimension represents a distinct optimization variable. The particles in the multidimensional search space are characterized by corresponding fitness values. They make movements in the search space towards regions characterized by high fitness values.
The conventional FRM digital filters incorporate FIR interpolation digital subfilters. These digital subfilters are usually of high orders, rendering the resulting overall FRM digital filters as not economical, since the resulting digital filters occupy large chip areas and consume high amounts of power in their VLSI hardware implementations. In general, the multiplication operation is the most cost-sensitive part in such an implementation. Therefore, there is every incentive to reduce the number of multiplication operations in the digital filter realization. This problem may be circumvented by employing IIR interpolation digital subfilters [29,30].
There is a vast body of literature available for the design and optimization of digital IIR filters [31][32][33]. However, all the aforementioned designs are based on the exact transfer function coefficients which leads to an uneconomical hardware realization of such filters. In order to realize the constituent IIR interpolation digital subfilters on a hardware platform, the bilinear-lossless-discrete-integrator (bilinear-LDI) digital filter design approach is employed [34]. These digital subfilters are realized as a sum/difference of a pair of bilinear-LDI digital allpass networks. The salient features of the bilinear-LDI digital filters are that they lend themselves to fast two-cycle parallel digital signal processing speeds, while being minimal in the number of digital multiplication operations (and, practically, minimal in number of digital addition and unit-delay operations).
The starting point in the design of FRM digital filters is to find the multiplier coefficients constituent in the FRM digital filter in infinite precision by using the hitherto gradient-based optimization techniques (e.g. Parks-McClellan approach [35] for FIR digital filters) followed by a quantization step. The quantization can be performed by constraining the multiplier coefficients values to conform to certain number systems such as the signed power-of-two (SPT) system. SPT is a computationally efficient number system which can further reduce the hardware complexity of the FRM IIR digital filters. In this number system, each multiplier coefficient is represented with only a few non-zero bits within its wordlength, permitting the decomposition of the multiplication operation into a finite series of shift and add operations. Digital filters incorporating SPT multiplier coefficient representation are commonly referred to as multiplierless digital filters [36]. However, the SPT representation of a given number is not unique, resulting in redundancy in the multiplier coefficient representation. This redundancy can adversely affect the corresponding computational complexity due to recourse to compare operations repetitively.
The canonical signed digit (CSD) number system is a special case of the SPT number system which circumvents the above redundancy problem by limiting the number of non-zero bits in the representation of the multiplier coefficients. It is usually used in combination with subexpression sharing and elimination, which in turn results in substantial reduction in the cost of the VLSI hardware implementation of the digital filters [37]. In CSD number system, no two (or more) non-zero bits can appear consecutively in the representation of the multiplier coefficients, reducing the maximum number of non-zero bits by a factor of two in terms of shift and add operations [38].
After multiplier coefficient quantization, the resulting FRM digital filter may no longer satisfy the given target design specifications. Therefore, the next step in the design of FRM digital filters is to perform a further optimization to make the finite precision FRM digital filter to conform to the design specifications. This can be achieved by resorting to a finite-precision optimization technique such as PSO.
A direct application of the conventional PSO algorithm to the optimization of the above FRM digital filters gives rise to three separate problems: • The first problem arises because in the course of optimization, the multiplier coefficient update operations lead to values that may no longer conform to the desired CSD wordlength, etc. (due to random nature of velocity and position of particles). This problem is resolved by generating indexed look-up tables (LUTs) of permissible CSD multiplier coefficient values, and by employing the indices of LUTs to represent FRM digital filter multiplier coefficient values.
• The second problem stems from the fact that in case of FRM IIR digital filters, the resulting FRM IIR digital filters may no longer be bounded-input-bounded-output (BIBO) stable. This problem can be resolved by generation and successive augmentation of template LUTs until the BIBO stability constraints remain satisfied [23].
• Finally, the third problem arises because even in case of having indexed LUTs, the particles may go over the boundaries of LUTs in course of optimization (due to the inherent limited search space). This can be resolved by introducing barren layers. A barren layer is a region, with a certain width and certain entries, which is added to the problem space such that the particles tend to shy away from such a region. The width of the barren layers is calculated based on a worst case scenario that may happen in the particles movements in the search space. However, the entries of barren layers are different for different problems and depend on the topology of the search space and the fitness function used in the problem.
This chapter discusses in detail the design, realization and discrete PSO of FRM IIR digital filters. FRM IIR digital filters are designed by FIR masking digital subfilters together with IIR interpolation digital subfilters. The FIR filter design is straightforward and can be performed by using hitherto techniques. The IIR digital subfilter design topology consists of a parallel combination of a pair of allpass networks such that its magnitude-frequency response matches that of an odd order elliptic minimum Q-factor (EMQF) transfer function. This design is realized using the bilinear-LDI approach, with multiplier coefficient values represented as finite-precision CSD numbers.
The above FRM digital filters are optimized over the discrete multiplier coefficient space, resulting in FRM digital filters which are capable of direct implementation in digital hardware platform without any need for further optimization. A new PSO algorithm is developed to tackle three different problems. In this PSO algorithm, a set of indexed LUTs of permissible CSD multiplier coefficient values is generated to ensure that in the course of optimization, the multiplier coefficient update operations constituent in the underlying PSO algorithm lead to values that are guaranteed to conform to the desired CSD wordlength, etc. In addition, a general set of constraints is derived in terms of multiplier coefficients to guarantee that the IIR bilinear-LDI interpolation digital subfilters automatically remain BIBO stable throughout the course of PSO algorithm. Moreover, by introducing barren layers, the particles are ensured to automatically remain inside the boundaries of LUTs in course of optimization.

The conventional PSO algorithm
Let us consider an optimization problem consisting of N design variables, and let us refer to each solution as a particle. Let us further consider a swarm of K particles in the N-dimensional search space. The position of the k-th particle in the search space can be assigned a N-dimensional position vector X k = {x k1 , x k2 , . . . , x kN }. In this way, the element x kj (for j = 1, 2, . . . , N) represents the j-th coordinate of the particle X k .
The PSO optimization fitness function maps each particle X k in the search space to a fitness value. In addition, the particle X k is assigned a N-dimensional velocity vector The PSO optimization search is directed towards promising regions by taking into account the velocity vector V k together with the best previous position of the k-th particle X best k = {x best k1 , x best k2 , . . . , x best kN }, and the best global position of the swarm G best = {g best 1 , g best 2 , . . . , g best N } (i.e. the location of the particle with the best fitness value).
The conventional PSO is initialized by spreading the particles X k through the search space in a random fashion. Then, the particles make movements through the search space towards regions characterized by high fitness values with corresponding velocities V k . The movement of each particle is governed by the best previous location of the same particle X best k , and by the global best location G best . The velocity of particle movement is determined from the previous best location of the particle, the global best location, and the previous velocity.

Figure 1. Movement of Particles in PSO Algorithm
The velocity and position of each particle in the i-th iteration throughout the course of PSO are updated in accordance with the equations: The parameter w represents an inertia weight; c 1 and c 2 are the correction (learning) factors, and r 1 and r 2 are random numbers in the interval [0, 1]. The velocity is limited between v min and v max to avoid very large particle movements in the search space, where v min < 0 and v max > 0. Fig. 1 illustrates how the particles move in a two-dimensional search space (N = 2). In this figure, two particles are present in the swarm, i.e. K = 2.
The first term in the right hand side of movement update Eqn. (1), weighted by w, signifies the dependence of the current particle velocity on its value in the previous iteration. The second term, weighted by c 1 , signifies an attractor to pull the particle towards its previous best position. The third term, weighted by c 2 controls the movement of the particle towards the global best position.
In addition to the update Eqns. (1) and (2), one can limit the coordinates in a particle between two user defined values x j min and x j max in order to limit the search space. However, This operation increases the complexity and consumes time.

Design of lowpass FRM digital filters
The block diagram in Fig. 2 Here, z represents the discrete-time complex frequency, and ω represents the corresponding (normalized) real frequency variable. Moreover, F 0 (z) and F 1 (z) represent FIR masking digital subfilters, while H a (z M ) and H b (z M ) represent M-fold interpolated versions of H a (z) and H b (z), respectively. In case of FIR digital interpolation subfilters, for a linear-phase filter H a (z) of order N FIR , the relationship between H b (z) and H a (z) is as follows: and hence H b (z) can be implemented by subtracting the output of H a (z) from the delayed version of the input, as shown in Fig. 3.
The FRM digital filter in Fig. 2 has an overall transfer function

Filter Passband Edge Stopband Edge
Case I The masking digital subfilters F 0 (z) and F 1 (z) are employed to suppress the unwanted image bands produced by the interpolated digital subfilters H a (z M ) and H b (z M ). The masking filters are made to have equal order (by zero padding) in order to ensure that their phase characteristics are similar. The corresponding interpolated digital subfilters H a (z M ) and H b (z M ) can realize transition bands which are a factor of M sharper than those of H a (z) and H b (z), without increasing the number of required non-zero digital multipliers. The magnitude frequency-response of the various subfilters incorporated by the FRM digital filter design approach are shown in Fig. 4.
Here, Case I design is when the transition band of H(z) is extracted from that of H a (z M ) and Case II design is when the transition band of H(z) is extracted from that of H b (z M ). The edge frequencies of the overall digital FRM filter and its constituent subfilters are given in Table 1, where I L represents the number of image lobes to be masked given by: Case II (6) where ⌊ ⌋ denotes the largest integer from the lower side, and ⌈ ⌉ signifies the smallest integer from the upper side.

Design of bandpass FRM digital filters
In general, it is possible to extend the conventional FRM approach for the design of bandpass or bandstop FRM digital filters. However, the resulting FRM digital filters are constrained to have identical lower and upper transition bandwidths. In [39], this restriction was relaxed by realizing the bandstop FRM FIR digital filter as a parallel combination of a corresponding pair of lowpass and highpass FIR for Case II [3].
where H l p (z) represents a lowpass and H hp (z) represents a highpass FRM digital filter. In this way, H l p (z) and H hp (z) can be obtained with the help of Eqn. (5) as The lower transition bandwidth is governed by the constituent transition bandwidth of the highpass FRM digital filter, while the upper transition bandwidth is governed by the constituent transition bandwidth of the lowpass FRM digital filter. The realization for bandpass FRM digital filter are as shown in Fig. 5.

Design of FRM digital filters incorporating IIR interpolation digital subfilters
In the case of FRM IIR digital filters, H a (z) and H b (z) (in section 3) act as IIR interpolation digital subfilters. The masking filters F 0 (z) and F 1 (z) are not changed (i.e. they are still equal order FIR digital filters). Therefore, Eqn. (5) is still valid for the FRM IIR digital filter.
The IIR interpolation digital subfilter H a (z) is chosen to have an odd order N I IR . Odd-ordered elliptic transfer functions can be represented as a sum of or difference between two allpass transfer functions [40]. Therefore, H a (z) can be realized as the addition of two allpass digital networks G 0 (z) and G 1 (z) as follows: where G 0 (z) is odd-ordered and G 1 (z) is even-ordered. The interesting fact is that the difference between G 0 (z) and G 1 (z) results in a filter that is power complementary to H a (z), and can subsequently be used as the power complementary interpolation digital subfilter H b (z) as in the following:  It can be easily verified that H a (z) and H b (z) are power complementary digital filters [29], i.e. they satisfy Eqn. (3). In addition, it is well known that this structure halves the number of multiplier coefficients required for the implementation of FRM digital filters and therefore is the most economical realization since it requires a total of only N I IR multiplier coefficients to realize both H a (z) and H b (z).
The overall transfer function of H(z) given by Eqn. (5) can be expressed as: The block diagram in Fig. 6 shows the IIR interpolation digital subfilters H a (z) and H b (z) realized as a parallel combination of two allpass networks. It should be noted that if H a (z) is a lowpass filter, H b (z), which is the power complementary of H a (z), is a highpass filter. Fig. 7 shows an overall FRM IIR digital filter realization.
One may rearrange the structure in Fig. 7 by using Eqns. (10)(11). This can be performed by defining two digital subfilters as follows: Then H(z) in Eqn. (12) simplifies to: The advantage of realizing the FRM IIR digital filter as shown in Fig. 8 is that two adders shown in Fig. 7 are removed and they are no longer required. This subsequently simplifies the hardware implementation of the overall FRM IIR digital filter. However, it should be noted that the FIR masking digital subfilters F 0 (z) and F 1 (z) are made to be equal order using zero padding, and this results in the masking filters being moderately sparse. This is not the case when A(z) and B(z) are used instead. Therefore, the gain in hardware that could be achieved by using the realization in Fig. 8 is offset by a greater number of non-zero multiplier coefficients required in the realization of FRM IIR digital filters.

Realization of IIR interpolation digital subfilters using Elliptic Filters with Minimum Q-factor (EMQF)
Bilinear-LDI transformation falls into the category of digital filter realization techniques that transform an analog reference filter to its digital counterpart. Despite all the advantages of EMQF filters, they suffer from not being able to independently specify passband and stopband ripples [41], [42] of the filter. Additionally, EMQF filters have exceedingly low passband attenuation.
All the poles of an EMQF transfer function reside on a circle in the s domain rendering them to have equal magnitudes. Given a squared passband and stopband tolerance of δ p and δ a , respectively, for an EMQF filter, the passband ripple ∆ p and minimum stopband attenuation ∆ a can be obtained as follows [43]: The required passband and stopband edge frequencies for the analog reference filter H a (s) can be determined using design specifications along with Table 1. Frequency wrapping from digital to analog domain, and vice versa, has to be taken into account in accordance with: where Ω A is the analog frequency variable, where ω d is the digital frequency variable, and where T is the sampling period.
Once the transfer function of the analog reference filter H a (s) is determined, it is represented as a sum of two allpass analog filters G 0 (s) and G 1 (s). In addition, H b (s), which is the power complementary of H a (s) is represented as the difference of G 0 (s) and G 1 (s). The poles of G 0 (s) and G 1 (s) are determined by cyclically distributing the poles of the reference filter H a (s) [43]. In the next section, belinear-LDI design technique is used to transform the two allpass networks G 0 (s) and G 1 (s) into digital domain.

Implementation of EMQF interpolation subfilters using bilinear-LDI design approach
In this section, the design procedure in [34,44] is briefly explained to design and implement digital filters G 0 (z) and G 1 (z) using the the bilinear-LDI approach. This approach transforms analog reference filters G 0 (s) and G 1 (s) to obtain their digital filter counterparts G 0 (z) and G 1 (z).
The bilinear frequency transformation maps the analog frequency variable s to its digital domain counterpart z in accordance with: where T represents the sampling period, for mapping the transfer function of a prototype reference filter from the analog domain to the digital domain. The bilinear transform maps the left half of the complex s-plane to the interior of the unit circle in the z-plane. Therefore, BIBO stable filters in the s domain are converted to filters in the z domain which preserve that stability. Similarly, if the analog reference filter is minimum-phase, the previous characteristic of bilinear transform guarantees that the resulting digital filter is also minimum-phase. It also preserves the sensitivity properties of the analog reference filter. However, bilinear transform may result in a digital filter that has delay-free loops in its implementation. Unfortunately, delay-free loops prevent the implementation of a digital filter to be realizable in hardware platform.
The LDI frequency transformation ensures the absence of delay-free loops in the digital implementation and is given by The LDI frequency transformation maps the hardware implementation of the analog reference filter to digital domain. While the LDI frequency transformation guarantees that there are no delay-free loops in the implementation of the digital filter, it does this to the cost of resulting in a digital filter having poor magnitude-frequency responses. Moreover, it is incapable of preserving the BIBO stability properties of the analog reference filter.
The bilinear-LDI approach is a combination of the two above mentioned realization techniques. In bilinear-LDI transform, a precompensation is performed to the reference analog filter. Then, the conventional LDI design technique is applied to a network resulting from the precompensated analog prototype filter. The precompensation is such that the application of the LDI design technique results in a filter that exactly matches the bilinear frequency transform of the uncompensated analog prototype filter.
The resulting bilinear-LDI digital filters have several desirable features from a hardware realization point of view. They are minimal in the number of digital multiplication operations. Although they are not minimal in the number of digital adders and unit-delays, the additional adders and the additional unit delay lead to certain advantages when the concept of generalized delay unit is used for the realization of the network [34]. Moreover, The bilinear-LDI digital filters lend themselves to fast two-cycle parallel digital signal processing speeds and they exhibit exceptionally low passband sensitivity to their multiplier coefficient values, resulting in small coefficient wordlengths.
As discussed in Section 5, the analog reference filter H a (s) is decomposed into two allpass analog networks G 0 (s) and G 1 (s). The digital allpass networks G 0 (z) and G 1 (z) are obtained from G 0 (s) and G 1 (s) using the bilinear-LDI design approach.
It should be pointed out that G 0 (s) is an odd-ordered allpass function. Therefore, it has a pole on the real axis in the s domain. On the other hand, G 1 (s) ends up having an even-ordered allpass function. It is well known that an allpass transfer function can be written in the general form [34]: where P(s) is a Hurwitz polynomial of order, say,ñ . Moreover, P(s) can be expressed as: where EvP(s) denotes the even and OdP(s) denotes the odd part of P(s).  By simple manipulation of Eqns. (21) and (22) one can get Here,K = 1 or −1, and Z(s) is a realizable reactive impedance given by

OdP(s) EvP(s) for evenñ
EvP(s) OdP(s) for oddñ (24) whereñ is the order of G(s) (odd when realizing G 0 (s) and even when realizing G 1 (s)). The impedance Z(s) has a zero at s = 0 for evenñ and a pole at s = 0 for oddñ, while having a zero at s = ∞ both for evenñ and for oddñ.
The bilinear-LDI digital realization of G(s) is achieved by using the following steps: • The transfer function G(s) is decomposed in the form where Here, G(s) can be realized as the transfer function of the signal-flow graph in Fig. 9. Furthermore, g(s) represents a lowpass or highpass analog filter that can be realized as the transfer function of the voltage divider network in Fig. 10.

Figure 11. Realization of Impedance Z(s)
Finally, Z(s) represents realizable reactances (consisting of capacitors and inductors only) and can be decomposed into its Foster II canonical form, as in Fig. 11, in accordance with where m =ñ/2 for evenñ and m = (ñ + 1)/2 for oddñ, and where C i represent capacitances and L i represent inductances (for i = 1, 2, . . . , m), and inductor L 1 is only present for evenñ.
• The impedance Z(s) in Fig. 11 is substituted into Fig. 10 and the precompensation is applied to the resulting network. This amounts to a modification of circuit elements in accordance with: The resistance in r 0 in Fig. 10 is modified to: and with r 0 = 1Ω and for i = 2, 3, ..., m.
• Since the voltage/current signal-flow graph of the precompensated network [34] consists of analog integrators only and it has no analog differentiators, it can be used for bilinear-LDI realization method. Therefore, the analog integrators in the signal-flow graph of the precompensated network are replaced by LDI digital integrators, and by impedance-scaling, the resulting network is scaled by z − 1 2 to eliminate any half-delay elements. The resulting digital network is displayed in Fig. 12. The multiplier coefficients in Fig. 12 are as follows: for i = 1, 2, ..., m.

Constraints for guaranteed BIBO stability
In order for the FRM digital filter consisting of CSD multiplier coefficientsm FRM to be BIBO stable, it is both necessary and sufficient for the bilinear-LDI IIR interpolation digital subfilters H a (z) and H b (z) to be BIBO stable. Likewise, in order for the interpolation digital subfilters H a (z) and H b (z) to be BIBO stable, it is both necessary and sufficient for the bilinear-LDI allpass digital networks G 0 (z) and G 1 (z) to be BIBO stable. In this way, it is required that the bilinear-LDI digital allpass networks G 0 (z) and G 1 (z) remain BIBO stable throughout the course of the PSO algorithm.
In the course of PSO algorithm, the infinite-precision multiplier coefficients m L i and m C i can only take quantized valuesm L i andm C i that belong to CSD(L, l, f ). In order for the bilinear-LDI digital allpass networks G 0 (z) and G 1 (z) to remain BIBO stable, it is required that the values of the corresponding quantized reactive elementsL i andĈ i remain positive [45] in the course of optimization. This is due to the properties of the bilinear frequency transformation from analog to digital domain. In order to find the conditions for BIBO stability and in accordance with Eqns. (35) and (36), one has: Moreover, in accordance with Eqns. (31)(32)(33)(34), one has: whereL 1 = ∞ for odd-ordered allpass network G 0 (z).
By substituting Eqns. (37) and (38) into Eqns. (39)(40)(41)(42), and by solving the resulting equations for the reactive elementsL i andĈ i , one can obtainL From Eqns. (43)(44)(45)(46),L i > 0 andĈ i > 0 provide that Then, in order to make the CSD FRM digital filter BIBO stable, it is necessary and sufficient to choose the values of the multiplier coefficientsm FRM ∈ CSD(L, l, f ) such that the inequality constraints (47-50) are satisfied. The equations and corresponding condition required for BIBO stability are summarized in Table 2.
In order to make the CSD lowpass digital IIR FRM filter BIBO stable, it is necessary and sufficient to choose the values of the multiplier coefficientsm L i ,m C i ∈ CSD(L, l, f ) such that the inequality constraints of Table 2 are satisfied.
It should be pointed out that constraint (49) is most stringent whenm L i is at its largest possible value. Similarly, constraint (50) is most stringent whenm L 1 ,m L i andm C i are all at their largest possible values (whilem L i andm C i still adhere to constraintm C i < 4 (m L i ) −1 ).

Proposed PSO of FRM IIR digital filters
The proposed particle swarm optimization of BIBO stable FRM IIR digital filters is carried out over the CSD multiplier coefficient space CSD(L 0 or 1 , l 0 or 1 , f 0 or 1 ), where L 0 or 1 represents the multiplier coefficient wordlength, where l 0 or 1 represents the maximum number of non-zero digits, and where f 0 or 1 represents the number of fractional part digits (for FIR or IIR digital subfilters, respectively).

Element Equation Inequality
Constraintŝ The starting point of any stochastic algorithm plays an important role in the convergence behavior of the optimization algorithm [46]. Therefore, it is important to generate the initial swarm in proper positions in the search space rather than complete random generation of the initial population. In order to achieve this, the following technique is employed:

Initiation of PSO
To start the PSO algorithm from a good position in the search space the infinite precision multiplier coefficient values of the seed particle are generated by using classical techniques as discussed in previous sections. These infinite precision multiplier coefficient values are turned into their finite precision counterparts by simply rounding them to their closest CSD values. This seed particle is used as the center of the swarm and a cloud of particles are generated randomly around the seed particle. It should be noted that the distance of the randomly generated particles should not be far from the seed particle. In this way, the initial swarm contains particles which have high chances of being near the optimal solution. The multiplier coefficient values of the swarm are taken from a set of CSD LUTs which are constructed as follows:

FRM IIR digital filter template LUTs
It is necessary and sufficient to choose the values of the multiplier coefficients, such that the inequality constraints (47-50) are satisfied. In order to achieve this, the LUTs are constructed as follows: • One LUT is constructed for all multiplier coefficient valuesm FIR ∈ CSD(L 0 , l 0 , f 0 ) for the masking digital subfilters F 0 (z) and F 1 (z). The values of L 0 , l 0 and f 0 are determined empirically based on the amplitude frequency-response of the masking digital subfilters F 0 (z) and F 1 (z).
• A LUT is constructed for all multiplier coefficient valuesm I IR ∈ CSD(L 1 , l 1 , f 1 ) for the digital allpass networks G 0 (z) and G 1 (z). Once again, the values of L 1 , l 1 and f 1 are determined empirically. Also, it is expedient to assume thatm I IR have only positive values.
• The above CSD LUT is used to form one size-reduced LUT per the multiplier coefficient for digital allpass networks G 0 (z) and G 1 (z), where each size-reduced LUT initially includes CSD values bounded from below by the smallest representable value belonging to CSD(L 1 , l 1 , f 1 ), and from above by the corresponding value of the finite-wordlength coefficients for the seed particle. The size-reduced LUTs are augmented before PSO process commences. The purpose of this augmentation is to ensure that the exploration space include as many of those CSD multiplier coefficientsm L 1 ,m C 1 ,m L i andm C i which still satisfy the BIBO stability constraints (47-50).
The above constructed LUTs are used as template LUTs. There are two problems concerning the PSO of FRM IIR digital filters over the CSD multiplier coefficient space. To overcome these problems, the template LUTs must be further processed. These two problems and the way to solve them are discussed in the following.

PSO indirect search method
In PSO, the required new particle position is obtained from the previous position of the particle through the addition of a random (normalized) velocity value. However, by directly applying the conventional PSO to the above optimization over the CSD multiplier coefficients, one may obtain new particle positions whose coordinate values are no longer in CSD(L 0 or 1 , l 0 or 1 , f 0 or 1 ). In order to overcome this problem, the optimization search is carried out indirectly via the indices to the LUT CSD values (as opposed to LUT CSD values themselves). In this way, the CSD coordinate values for each particle position are obtained by integer indices to the CSD LUTs. The key point in the indirect search rests with ensuring that the index set is closed, i.e. by ensuring that each index points to a valid CSD value in the LUT, and that the resulting particle in the course of PSO adheres to the prespecified CSD number format.
If the velocity values are replaced by their closest integer values, the update equations become modified tov Here,x kj ,v kj ,x best kj ,ĝ best j ,v min andv max are all integer values wherev min < 0 andv max > 0. In addition, w is limited in the interval [0, 0.5) (as discussed shortly).

Barren layers
Due to their finite length, the template LUTs inevitably lead to a bounded optimization search space. In order to ensure that the particles do not cross over to the outside of the search space in the course of PSO, the search space is constructed as a combination of two regions, namely the interior and barren layers. The barren layer is constructed to yield relatively low fitness values, and is represented as header and footer in the template LUT. There are two problems concerning the construction of the barren layers:

barren layer entries
The first problem in the construction of barren layers concerns how to make the fitness values in the barren layer relatively low. This problem can be resolved by filling the header part by unrealistically large, and the footer part by unrealistically small CSD multiplier coefficient values.

barren layer width
The second problem, on the other hand, concerns how to determine the width of the barren layer such that the particles do not cross over to the outside of the search space even under the worst case scenario. These two problems relate to the number of entries and the CSD values of the entries in header and footer parts of the template LUTs. To overcome this problem, let us consider the j-th variable in the k-th particle is in the boundaries of one of the template LUTs in iteration i − 1. The worst case scenario occurs when x i−1 kj moves toward the barren layer with the peak permissible velocities (v max for the header, and v min for the footer). If in the i-th iteration x i kj is in the footer: and if it is in the header: Eqns. (53-56) show that the velocity of the particle in iteration i + 1 tends to move the particle in a direction opposite to the direction of the barren regions. Here, the worst case happens when r 1 = r 2 = 0. In this way, the number of entries L f in the footer part, and the number of entries L h in the header part is determined in accordance with Let us recall that since 0 ≤ w < 0.5, In addition, after some iterationsv i+1 kj = 0. Otherwise, if w ≥ 0.5,v i+1 kj can never become zero, and the width of the barren layer will be infinity.
The augmented LUTs remains fixed in the course of PSO, restricting automatic particle movement inside the limited search space. Modifying the index values inside each particle by adding the current indices to the length of the footer barren region, L f , PSO algorithm is ready to start the optimization of FRM digital filters.

Design methodology
The design methodology for the proposed PSO of BIBO stable bilinear-LDI based FRM IIR digital filters over the CSD multiplier coefficient space can be summarized as follows: 1. Designing the interpolation digital subfilter: the first step in determining the interpolation subfilter specifications is to fix the interpolation factor M from a pre-specified range. This is done in a way that the order of the FIR masking filters is kept minimal. Using the passband edge frequency ω p and stopband edge frequency ω a and the expressions for boundary frequencies given in Table 1, one can determine the filter case and calculate the approximate passband edgeθ and stopband edgẽ φ of the digital interpolation lowpass subfilter H(e jω ), for every value of the user specified range of interpolation factors M. The order of the FIR masking filters depends on the minimum distance between consecutive image replicas of either the interpolated subfilter H a (e jMω ) or its complement H b (e jMω ). Then, displacement λ M and distanceD M for each interpolation factor M are given as: To minimize the length of FIR-masking filters, the value of M that results in the largest value ofD M is chosen. This determines the optimal interpolation factor M as well as the approximate passband edgeθ and stopband edgeφ of the digital interpolation subfilter H(e jω ). EMQF filters have the property of equal square magnitude ripple size in the passband and stopband. Therefore, of the two ripple specifications, whichever gives the smallest tolerance in the squared magnitude response determines both the passband ripple R p and stopband attenuation R a of the interpolation digital subfilter H a (e jω ). The interpolation digital subfilter order N I IR is then determined using R p , R a , θ andφ. N I IR must be rounded to the nearest larger odd integer so that it can be implemented by a parallel combination of two allpass networks. With the order N I IR , and passband and stopband ripples R p and R a fixed, the ratio of the analog passband edge θ A and stopband edge φ A is a constant k given by [47] In order to satisfy the passband edge specification, the digital passband edge ω p =θ for Case I filters. The digital stopband edge ω a is then determined using the analog ratio k. (Here, frequency warping from digital to analog domain, and vice versa, given by Eqn. (18) needs to be taken into account.) Similarly, ω a =φ for Case II filters, and ω p can be determined by using ratio k. Also, using given ripple specifications along with the boundary frequencies described in Table 1, one can determine the transfer function of the FIR masking filters F 0 (e jω ) and F 1 (e jω ).
2. Generation of seed FRM digital filter particle: The seed FRM digital filter particle is formed as follows: • A particle with B 1 coordinates is formed in which each coordinate serves as an index of the corresponding CSD LUT for each multiplier coefficient constituent in the interpolation digital subfilters. For FRM IIR digital filters, the multiplier coefficients correspond to the bilinear-LDI allpass digital networks G 0 (z) and G 1 (z). • A particle with B 2 coordinates is formed in which each coordinate serves as an index of the corresponding CSD LUT for each multiplier coefficient in the FIR masking digital subfilters F 0 (z) and F 1 (z).
3. Generation of Initial Swarm: An initial swarm of K particles is formed by generating a random cloud around the seed particle as discussed in section 8.1.

Fitness Evaluation:
The fitness function for CSD FRM IIR digital filters is defined in accordance with f itness group−delay = ς p (63) where with ∆ω p representing the passband frequency region(s), with ∆ω a representing the stopband frequency region(s), and with τ(ω) representing the group-delay frequency response of the FRM IIR digital filter. Here, W p , W a , and W gd represent weighting factors for the passband and stopband magnitude responses, and for the group-delay response, respectively. Moreover, µ τ represents the average group-delay over the passband region.
In [48], a convenient way to represent digital networks in terms of matrix representation is presented. This technique can be used to find the magnitude and group delay frequency response of the digital network in Fig. 12. Let us consider the input to the digital network in Fig. 12 to be x D and the output of it to be y D . In addition, let the output of the i-th time delay in Fig. 12 to be x i and the input to the i-th time delay to be y i . The transfer function matrix of the network, T, can be found as where y = [y D , y 1 , y 2 , . . . , y 2m+1 ] t 2 and x = [x D , x 1 , x 2 , . . . , x 2m+1 ] t , and T is a (2m + 2) × (2m + 2) matrix with the entries obtained as Eqn. (69).
Since x i = z −1 y i , the transfer function G(z) = y D x D can be found as where e is a row vector and c is a column vector of length 2m + 1, and where I is the identity matrix and D is a (2m + 1) × (2m + 1) matrix in accordance with  L 0 l 0 f 0 L 1 l 1 f 1 11 3 10 12 3 7

Bandpass FRM IIR digital filter design example
Consider the design of a bandpass FRM IIR digital filter satisfying the magnitude response design specifications given in Table 3 over the CSD multiplier coefficient space.
The parameters for the PSO of bandpass FRM IIR digital filter is shown in Table 4 and the CSD parameters are presented in Table 5.
Given the design specification in Table 3, The order of the digital allpass networks G 0 lp (z), G 1 lp (z), G 0 h p (z) and G 1 h p (z) are found to be 3, 4, 3 and 4, respectively. In addition, the digital masking subfilters F 0 lp (z), F 1 lp (z), F 0 h p (z) and F 1 h p (z) have the same length as the previous example, i.e. 24, 42, 25 and 35 respectively, resulting in N = 140. In this example a set of fifteen CSD LUTs are required, fourteen LUTs for the multiplier coefficients m C 0,1 , m C 0,2 , m C 0,3 , m L 0,2 , m L 0,3 , m C 1,1 , m L 1,1 , m C 1,2 and m L 1,2 constituent in the digital allpass networks G 0 lp (z), G 1 lp (z), G 0 h p (z) and G 1 h p (z), and one template LUT for all the multiplier coefficients constituent in the masking digital subfilters F 0 lp (z), F 1 lp (z), F 0 h p (z) and F 1 h p (z).
Finally, by using Parks McClellan approach, the subfilters F 0 lp (z), F 1 lp (z), F 0 h p (z) and F 1 h p (z) can be designed. Also, by using the EMQF technique, the digital allpass networks G 0 lp (z), G 1 lp (z), G 0 h p (z) and G 1 h p (z) can be designed. Consequently, the magnitude and group delay frequency responses of the overall infinite-precision bandpass FRM IIR digital filter H(z) is obtained as shown in Figs. 13 and 14.  By applying the proposed PSO to the initial FRM IIR digital filter and after about 160 iterations, the PSO converges to the optimal bandpass FRM IIR digital filter having a magnitude frequency response as shown in Fig. 17. In addition, Fig. 18 gives us a closer look to the magnitude frequency response of the passband region of the bandpass FRM IIR digital filter. Fig. 19 illustrates the group delay frequency response of the optimized bandpass FRM IIR digital filter. The values of the multiplier coefficients for the lowpass and highpass sections of the bandpass FRM IIR digital filter are obtained as summarized in Tables 6 and 7.    Table 8 represents the comparison of the CSD bandpass FRM IIR digital filters before and after PSO.