Optimized fractional-order Butterworth filter design in complex F-plane

This paper introduces a new technique to optimally design the fractional-order Butterworth low-pass filter in the complex F-plane. Design stability is assured by incorporating the critical phase angle as an inequality constraint. The poles of the proposed approximants reside on the unit circle in the stable region of the F-plane. The improved accuracy of the suggested scheme as compared to the recently published literature is demonstrated. A mixed-integer genetic algorithm which considers the parallel combinations of resistors and capacitors for the Valsa network is used to optimize the frequency responses of the fractional-order capacitor emulators as part of the experimental verification using the Sallen–Key filter topology. The total harmonic distortion and spurious-free dynamic range of the practical 1.5th-order Butterwoth filter are measured as 0.13% and 62.18 dBc, respectively; the maximum and mean absolute relative magnitude errors are 0.03929 and 0.02051, respectively.


Introduction
The Riemann-Liouville definition for the fractional derivative of arbitrary order υ for a function f (t) is given by where n ∈ Z + ; n −1 < υ ≤ n; a and t are the limits of the operation; and (·) denotes the gamma function [13,20,23]. For zero initial conditions, the Laplace transform of (1) is given by s υ F(s).
Applications of fractional calculus are pervading in various engineering disciplines [15], including fractional-order (FO) circuit theory and filter design [19]. Due to the presence of additional design parameters or 'tuning knobs' (viz., υ), FO filters exhibit amplitude, delay, and transient characteristics, which are not possible to achieve using the classical ones [2,11]. Generalization of integer-order filters to the FO-domain allows the exact order (typically, a real number) to be realized; thereby, providing precise roll-off and strict matching of design specifications [9]. Implementations of analog FO filters involve the fractance devices (e.g., FO capacitor/inductor) [4], whose immittance relationship is where in the case of FO capacitor, it is appropriate to consider I (s) as an admittance and k is the pseudo-capacitance expressed in Farad/second 1−β . In the case of FO inductor, I(s) can be considered as an impedance and k is the pseudo-inductance expressed in Henry/second 1−β . In both these cases, it holds β ∈ (0, 1). The frequency-domain expression for s β is given by ( jω) β = ω β ∠(βπ/2), where the analog frequency variable ω is expressed in radians per second (rad/s). Various approaches have been reported to emulate the frequency response of s β , [5]. The squared magnitude function of the normalized FO Butterworth filter (FOBF) [3] is given by where M = N + α, N ∈ Z + and α ∈ (0, 1). In the literature, two different approaches are reported to approximate the magnitude-frequency characteristic of the FOBF. The first technique involves mod-eling in the complex s-plane using the FO [3,10,25] or the integer-order transfer functions (see [16] and the references cited therein). Alternatively, FOBFs were also designed in the complex W -plane (W = s 1/q , where α = r /q; r , q ∈ Z + ) using the pole-placement algorithms [1,18].
A recent work has demonstrated the effectiveness of the F-domain-based approach in the reduced order modeling of FO commensurate system [17]. While realizing FO filters with optimization techniques exist in the literature [10,16,25], this paper is the first attempt to design the analog FOBF directly in the complex F-plane optimally. A constrained optimization strategy is formulated to guarantee design stability. Relative to the existing W -plane based techniques [1,18], the proposed method provides the following contributions: 1. Unlike [1], the symmetric distribution of coefficients in the denominator polynomial for the integer-order Butterworth filter is also exhibited by the proposed FOBF, 2. Higher-order FOBFs are practically realized in [1] by cascading Nth-order Butterworth filter with an αth-order one. In contrast, the proposed approach can directly generate the optimal FOBF model for any value of M, 3. Unlike [18], the all-pole structure of the Mth-order proposed model always requires N+2 terms, which is the same as that of the (N+1)th order Butterworth filter, and 4. Although the W -plane-based techniques [1,18] generate a stable design, these methods do not provide optimal pole placement in the W -plane. Consequently, the retro-fitted model in the s-plane is also not an optimal one. In contrast, the proposed constrained optimization approach guarantees both stability as well as optimal placement of poles in the unit circle of the complex F-plane for the normalized FOBF approximant. The corresponding s-domain-based proposed model achieves lower error as compared to the cited literature in approximating the theoretical magnitude response.
In the rest of the paper, Sect. 2 presents the proposed scheme; stability and modeling performances are analyzed in Sect. 3. Experimental results are presented in Sect. 4. Finally, conclusions are drawn in Sect. 5.

Proposed technique
It is well-known that for a function F = s β , where β ∈ (0, 1), the region of instability in the s-plane (|θ S | < 0.5π) is mapped into |θ F | < 0.5βπ in the complex F-plane [21]. Therefore, the region of stability in the F-plane is given by |θ F | > 0.5βπ.
Since s = F 1/β , then the squared magnitude function in the F-plane for the FOBF (normalized) can be obtained by substituting ω = F 1/β / j in (3), as represented by Selecting β = M/(N + 1) and substituting in (4) leads to the magnitude function as given below.
Thus, the F-domain-based M th -order Butterworth filter (normalized) is modeled as an all-pole polynomial function in F, as given by where u k (k = 1, 2, …) denotes the coefficients of D(F). Recall that the classical P thorder (P = 1, 2, …) Butterworth filter has an analogous representation in the s-domain [24], as given by H P (s) = 1 s P + a 1 s P−1 + a 2 s P−2 + · · · + a 2 s 2 + a 1 s + 1 .
The objective function for the constrained optimization problem to achieve a stable FOBF model in the F-plane is formulated as (8) Subject to: θ D(F) > θ C , where the frequency variable ω is varied with logarithmic spacing for L = 100 sample points between 10 −2 − 10 2 rad/s; θ C = 0.5βπ denotes the critical phase (in rad); θ D(F) represents the minimum absolute phase of the roots of D(F); and the decision variables vector is given by Thus, the dimension (d) of the problem is N/2 for even values of N and (N+1)/2 if N is an odd number. To maintain conformity with the classical Butterworth filter, the lower bound of X is set as 0 and the upper bound (Ub) is chosen as the values of the coefficients for the P th -order Butterworth filter, where P = N + 1, as given by Finally, the s-domain transfer function of the proposed FOBF (normalized) can be obtained as

Performance analysis using MATLAB simulations
The applicability of metaheuristic algorithms for solving constrained optimization problems is well established in the literature [14]. The search process of these algorithms involves a combination of stochastic and deterministic rules. These methods employ multiple agents that explore the hyper-dimensional problem search space in the initial stages, followed by local search in the later stages. This is unlike the 'fminsearch' function in MATLAB employed in [25], where the unconstrained optimization technique involves a direct search strategy based on the Nelder-Mead simplex algorithm. The iterative nature of the metaheuristic search process may lead to higher computational time as compared to the traditional optimization algorithms. It is worth noting that metaheuristic algorithms can only achieve a near-global optimal solution. Furthermore, a closed-form solution is also not attained by a metaheuristic technique. However, these methods have established a niche due to their ease of implementation and ability to effectively handle linear or non-linear, convex or non-convex, multimodal, multidimensional, unconstrained, and constrained optimization problems. The proposed optimization technique is implemented in MATLAB using the constrained composite differential evolution (C 2 oDE) algorithm [27]. C 2 oDE has demonstrated effective performance in solving the benchmark optimization problems provided in IEEE CEC2006 and CEC2010 test suites. In this work, the basic parameters for C 2 oDE are set as population size = 100 and termination condition = 10000d function evaluations, while the internal parameters are the same as specified in [27]. The accuracy and stability of the proposed FOBFs are investigated by considering various design orders. Table 1 presents the optimal value of X obtained for design cases such as M = 1.5, 2.5, 2.8, 3.2, 3.6, and 4.2. The following observations can be made about the pole locations and stability of the designed FOBFs based on the results shown in Table 2: 1. The roots of D(F) exhibit a magnitude of 1, which implies that all the poles of the proposed FOBF reside on the unit circle in the F-plane. A comparable situation exists for the pole locations of the classical Butterworth filter (normalized), where the poles are located on the unit circle in the left-half s-plane [24]. 2. The region of stability in the F-plane is larger than in the s-plane for β < 1. All the proposed models achieve stability in the F-plane since all the poles of H(F) have a negative real part. Note also from Table 2 that θ D(F) > θ C for all the    Table 3 shows the proposed models after transformation from F-domain to the s-domain. If M is an irrational number or it leads to the occurrence of a recurring decimal for β, then β may be truncated to three digits after the decimal point (viz., H(s) for M = 2.5, 2.8).
The absolute relative magnitude error (ARME) metric, as defined below, is used to evaluate the modeling accuracy.   [18,25]. The reported model in [18] exhibits peaking near the cut-off frequency, whereas the rolloff behavior of the FOBF in [25] deviates from the theoretical one. These findings are further confirmed from the ARME comparison plots presented in Fig. 2 (bottom). The magnitude response of the proposed 3.2th-order Butterworth filter is compared with that of the published literature [18], as shown in Fig. 3. The reported design exhibits peaking in the magnitude response near the normalized cut-off frequency, whereas, the proposed model stays in proximity to the theoretical characteristic throughout the design range.
The maximum (max) and mean (evaluated using 1000 sampled frequency points with logarithmic spacing in the interval [10 −3 , 10 3 ] rad/s) ARME comparisons with   Table 4. The proposed FOBFs achieve the least error for all the considered cases, thus demonstrating superior accuracy. The ARME comparison plots with [18] for M = 1.5, 2.5, 3.2, 3.6, and 4.2 are presented in Fig. 4 to justify the improved accuracy of the proposed approach graphically. Further investigations are conducted by considering the design of (2 + α)-order Butterworth filter with α varying from 0.01 to 0.99 in steps of 0.01. Results shown in  [25] for the proposed (2 + α)-order Butterworth filters Fig. 5 demonstrate that: (i) coefficient u 1 increases as α is increased and u 1 approaches 2 as M approaches 3. It is noteworthy that the third-order Butterworth polynomial is given by s 3 + 2s 2 + 2s + 1 [24], (ii) θ D(F) exceeds θ C throughout the design range, which justifies the stability criterion in the F-plane for the proposed FOBFs, and (iii) the proposed designs outperform the results in [25] regarding the mean ARME for M ∈ [2.10, 2.99].

Measurement results
The magnitude-frequency characteristic of the proposed 1.5th-order Butterworth filter with the transfer function as shown in Table 2 and a design cut-off frequency of 10,000 rad/s (f 0 = 1.5915 kHz) is verified experimentally. The filter is realized using the unity-gain Sallen-Key topology no. 1 of catalog in [22]. The operational amplifier (op-amp) based proposed circuit (see Fig. 6) employs fractional-order capacitors (FOCs) of an order β = 0.75 (ϕ = −67.5 • ) with pseudo-capacitances of C β1 = 100 nF·s −0.25 and C β2 = 1 µF·s −0.25 in place of the conventional capacitors.
A mixed-integer genetic algorithm [12], upgraded to consider parallel combinations of resistors and capacitors, was used for optimizing the phase and magnitude responses of the FOCs. The corresponding 6 branch Valsa RC ladder network-based FOC emulator [26] is also illustrated in Fig. 6. It provides a 7th-order approximation of the FOC, which has an admittance in the following general form Electronic Industries Alliance (EIA) standard E96 (1% tolerance) and E12 (10% tolerance) compliant resistor and capacitor values, respectively, which are commercially available in 0402 size RC kits [7,8], were used to construct the FOCs.  Table 5. The photographs of the fabricated FOC emulators with dimensions of 15 mm × 17 mm are depicted in Fig. 7. The ideal (target), simulated (nominal), and measured phase, magnitude, and pseudo-capacitance responses vs. frequency are shown in Fig. 8. Normalized histograms (%) of phase angle and pseudo-capacitance errors in the desired CPZ and the zoomed-in magnitude plot around the 10 3 Hz region are given as an inset. Simulations and Monte Carlo (MC) statistical analysis were performed in OrCAD Capture CIS 10.3 software. The MC was performed with passive element tolerances based on used datasheets for Gaussian distribution, with fixed random number seed value 17533 and 1000 runs to observe potential effects due to manufacturing processes.  The impedance measurements were carried out using the Digilent Analog Discovery 2 device (see Fig. 7) in the frequency range of 3 Hz to 300 kHz (1001 logarithmically spaced points per five decades) with a sinusoidal input signal of 100 mV rms . The performances of the FOCs in terms of relevant criteria, such as absolute phase angle deviation (APAD), phase band (PB), absolute pseudo-capacitance deviation (APCD), pseudo-capacitance band (PCB), etc., are reported in Table 6  The comparison of measured FOC behavior with the simulated (nominal) and MC indices highlight excellent matching in the whole desired CPZ.
The proposed FOBF was assembled on a breadboard and EIA standard-compliant E96 series (1% tolerance) through-hole TESLA Lanškroun resistors (R F1 = 3.6 k and 49.9 in series, R F2 = 2.74 k ) were used. The photograph of the hardware setup is presented in Fig. 9. The supply voltage of 5 V for the op-amp LT1360 [6] was provided from the Agilent E3630A power supply. The magnitude data points for 641 logarithmic spaced frequency points in the range [10 1 , 10 5 ] Hz were measured using the OMICRON Lab Bode 100 network analyzer and displayed using the Bode Analyzer Suite software. The level of the testing harmonic signal was set to 1.25 V PP . The receiver bandwidth of the analyzer was fixed at 300 Hz. The Table 6 Performance characteristics of the Valsa network in CPZ = 10 Hz to 100 kHz   time-domain behavior of the filter was observed on Agilent InfiniiVision DSO-X 2002A digital storage oscilloscope. A peak-to-peak sinusoidal voltage (V PP ) of 1.25 V was applied to the filter circuit from the Agilent 33521A function/arbitrary waveform generator.
The magnitude-frequency response measurement of the proposed Sallen-Key filter is compared with the ideal characteristic in Fig. 10 (top). It may be observed that the magnitude response of the practical filter conforms with the ideal behavior for the considered bandwidth. From the zoomed-in plot, it is confirmed that the magnitude of the proposed implementation at f 0 also stays in proximity with the theoretical (a) (b) Fig. 11 Experimental a input-output waveforms observed in an oscilloscope and b Fast Fourier Transform spectrum for the proposed 1.5th order Butterworth filter with an input frequency of f H = 1.59 kHz value (−3.0 dB). The ARME response of the constructed filter is presented in Fig. 10 (bottom). The mean ARME value is obtained as 0.02051, whereas the max ARME of 0.03929 occurs at 434 Hz. The time-domain response of the practical filter circuit is measured at the halfpower frequency f H = 1.59 kHz, as shown in Fig. 11a. For the applied V PP of 1.25 V, the peak-to-peak output voltage at f H was obtained as 0.88 V, which attains an excellent agreement with the theoretical anticipation of 0.88375 V. The Fourier spectrum of the measured output signal displayed up to the seventh harmonic is shown in Fig. 11b. Achieving a large Spurious-Free Dynamic Range (SFDR) and a small Total Harmonic Distortion (THD) for a practical filter circuit is desirable. The SFDR and the THD, evaluated from the plotted harmonics, are obtained as 62.18 dBc and 0.13%, respectively, which justifies the effectiveness of the proposed implementation.

Conclusions
An F-plane-based constrained optimization approach is presented for the first time to design the FOBF with guaranteed stability. Similarities with the classical Butterworth low-pass filter and advantages over the W -plane-based FOBFs are highlighted. Superior accuracy over the recent literature about the maximum and mean ARME is demonstrated through several examples. The proposed FOBF of order 1.5 was realized using the Sallen-Key filter employing fractional-order capacitor emulators. Experimental results exemplified an excellent agreement with the theoretical magnitude response. However, the proposed metaheuristic optimization approach based on C 2 oDE algorithm has two limitations: (i) it does not guarantee the generation of global optimal solution, and (ii) a closed-form solution is not attainable. Further research from both the optimization and the fractional calculus community is needed to alleviate these limitations.