On Digital Realizations of Non-integer Order Filters

Typical approach to non-integer order filtering consists of analogue design and implementation. Digital realization of non-integer order systems is susceptible to problems such as infinite memory requirement and sensitivity to numerical errors. The aim of this paper is to present two efficient methods for digital realization of non-integer order filters: discrete time-domain Oustaloup approximation and Laguerre impulse response approximation. Properties of both methods are investigated with use of non-integer low-pass filter. Filters realized with presented methods are then used for filtering of EEG signal. Paper concludes with discussion of merits and flaws of both methods.


Introduction
Non-integer (fractional) order systems are a rapidly growing field of interest for both mathematicians and engineers. One of the intensively analyzed aspects of this domain is non-integer order filters. One can observe that the theory of such filters is relatively well grounded; however, many problems of implementation are still open. A need for efficient implementation is obvious as potential applications are numerous in areas such as telecommunication, biomedical engineering, control and many others.
In this paper, the authors focus on efficient methods of creating discrete realizations of filters that do not have the burden of infinite memory and provide good representation of non-integer order filter frequency characteristics. In order to do so, two methods are proposed resulting in approximations in the form of discrete statespace systems. These methods are discrete time-domain Oustaloup approximation and Laguerre impulse response approximation (LIRA).
The problem of approximating the non-integer order system with an integer order one is being analyzed from many years. Most popular approaches are based on using Oustaloup transfer function approximation or Continuous Fraction Expansion (CFE) in the domain of transfer functions. Both of these methods were used for approximating integrators. The problem is that Oustaloup method is very sensitive to high discretization frequencies and rounding errors. It was observed in earlier authors' works [7] that it can become destabilized very easily. On the other hand, CFE method shows inferior quality in frequency characteristic representation [25]. Detailed analysis of CFE approximation in discrete time can be found in [12,44,45].
The first method used in this article, the Oustaloup method, is used in the literature in two different, yet equivalent, versions. The original approach, developed by Oustaloup [26,27], is based on approximation of fractional systems in frequency domain. This approach is widely used, e.g., [11,22,25,32] and many others. However, this method has some flaws which cannot be neglected-when discretized, it does not guarantee stability of the system (the poles of discrete system are outside unit circle) (see, e.g., [30]). In order to avoid, i.e., numerical issues induced by this method, another kind of approximation was proposed-instead of transfer function, the state-space approach is considered. State-space realization of Oustaloup transfer function was considered in, e.g., [23,36,39,40]; however, its discrete properties were not considered in this context. In this paper, the method from [7] is used. This approach allows avoiding numerical problems observed in practical realization. The proposed approach is to realize every block of the transfer function in form of a state-space system. Those first-order systems will be then collected in a single matrix resulting in full matrix realization. The results were first presented in authors' earlier works [7] and [6]. This approach proved to be more robust to different discretization schemes. The superiority of this approach was then validated with real-time control experiments-magnetic levitation control [10] and air heater control [9].
The other method considered in this paper, Laguerre impulse response approximation (LIRA), chosen by authors to analyze fractional filters uses approximation with (orthonormal) Laguerre functions. Some works, concerning this type of approximation, are, e.g., [3,24]. The authors' approach was developed independently in [5]. It introduces substantial improvements such as L 1 convergence, estimation of approx-imation error and choice of optimal approximation basis. For examples of use, see [46][47][48].
The most important results on non-integer order filters come from analog filtering theory. In [34], Radwan et al. generalize standard second-order filters. In particular, they analyze filters of various characteristics (i.e., low pass, high pass, band pass, notch) constructed using two fractional-order capacitors of order α. In [2], Soltan Ali et al. describe a transformation of Butterworth filter to non-integer order domain. They present both analytical and numerical results combined with Advanced Design System simulation. The result was a filter with more degrees of freedom than classical filter and therefore with better design flexibility. In [17], Freeborn et al. presented fractional two-step Tow-Thomas biquad filter. They used an approximation of fractional-order capacitors in form of integer order transfer function. The results for low-pass and asymmetric band-pass filters are verified both with MATLAB and PSPICE, and with physical experiments using non-integer order capacitors. Another example of secondorder filter generalization was described by Tripathy et al. [42]. Unlike earlier works, this paper uses two fractional-order capacitors with different orders. Experimental results were also compared with simulation (MATLAB, PSPICE), along with some analysis concerning stability of the filter and its sensitivity. In [35], the non-integer RLC filter is considered. The authors, A. G. Radwan and M. E. Fouda, use a gradientbased optimization technique to find optimal parameters in the frequency domain. They describe the performance of such filters along with some new issues which need resolving when designing non-integer order filters. Important results come from attempts to directly implement the non-integer order filters. For example, in [33], we can find an analysis of non-integer order oscillators. Radwan et al. derive there the Barkhausen condition for non-integer order system to oscillate. They also verify numerically the results as examples taking, i.e., Wien oscillator, Colpitts oscillator and phase-shift oscillator. In [1], we see new techniques for implementing continuous-time secondorder band-pass filters with high-quality factors and asymmetric slopes using multiple amplifier biquads and frequency-dependent negative resistor. Ahmadi et al. present there four possible realizations of the filters: one based on a frequency-dependent negative resistor (FDNR), another based on an inductor and two based on multiple amplifier biquads (MABs). In [41], there are another examples of design and implementation of fractional-order filters, along with performance analysis-in both simulation and experiment. This research stated clearly the superiority of non-integer order systems over classical systems. Tsirimokou et al [43] analyzed the behavior of so-called companding filters in low-voltage environment also using alternative domain, sinh instead of classic log-one. In [16], Freeborn et al. propose the use of field-programmable analogue array hardware to implement an approximated fractional step transfer function. They compare the experimental results with simulation performed in MATLAB. In order to implement the filters, they approximate the fractional operator s α with integer order transfer function.
Another approach consists of designing non-integer order filters directly in discrete domain. For example, in Stanislawski et al. [37,38] use Grünwald-Letnikow difference to design discrete filters. Another proposition of non-integer order discrete filter can be found in [13].
Main contribution of this paper is the use of efficient and numerically robust methods for implementation of non-integer order filters. The methods use the state-space formalism allowing easy implementation in the form of both MATLAB/Simulink blocks and direct coding. Methods have different underlying principles which allow covering wider area of filters allowing choice of proper method for appropriate task. Both methods allow very fast sampling frequencies without loss of stability. Methods are illustrated with use of low-pass generalization of classical second-order filter (also known as bi-fractional filter).
Rest of the paper is organized as follows. Firstly, certain preliminaries are presented, containing definitions and properties of necessary elements. This section includes also a brief description of analyzed filter, along with certain results on its stability and general remarks on discrete filter realization. Next, the first of considered methods is presented, the Oustaloup method, and general discussion on its digital implementations is given, including the description of some issues observed during discretization. The next section provides theoretical background for another method used for approximation based on impulse response of analyzed filter. Then, these two methods are compared basing on proposed performance indicator in H ∞ space. As an example of implementation, the results of filtering an EEG signal are presented. Paper ends with conclusions and some propositions on further works.

Remark 1
In this paper, we will use the following notation: Scalars will be marked with lowercase slanted letters (e.g., α ∈ R); vectors with lowercase bold-upright (e.g., x ∈ R n ) and matrices with uppercase bold-upright (e.g., A ∈ R n×n ).

Preliminaries
This section includes some theoretical preliminaries of the paper. In particular, necessary definitions of non-integer derivatives and function spaces are given.
Also the filter used to illustrate the operation of considered methods is described. Finally, general considerations regarding discrete realization of non-integer order filters are presented.

Non-integer Order Calculus
There are at least three widely used definitions of non-integer order derivatives [31]. In this paper, we will use the Caputo derivative of order α: where Γ (·) denotes the gamma function Another definition is useful for practical and numerical implementations-Grünwald-Letnikow definition of fractional derivative of order α It should be noted that Caputo and Grünwald-Letnikow derivatives used in non-integer differential equations give the same results only for 0 initial conditions. The third popular definition-Riemann-Liouville-is not used in this paper, and it is, however, equivalent to Grünwald-Letnikov one. Also we will denote by L 1 the Banach space of absolute integrable functions, and by L 2 the Hilbert space of square integrable functions with scalar product and norm x 2 2 = x|x . Additionally, H ∞ will denote Hardy space of bounded holomorphic functions, with the norm

Non-integer Low-Pass Filter
In their celebrated paper [34], Radwan et al. have presented a non-integer generalization of classical second-order filter section. In this paper, we focus on the low-pass variant (known also as bi-fractional filter). Such filters are a class of non-integer filters fully characterized by three parameters, α, ξ and ω 0 , and are given by the following transfer function Equivalent representation of (6) is the realization in the form of a system of differential equations of order α. This system can take form (see [19]) with the following matrices  (7) and zero initial conditions x(0) = 0 ∈ R 2 . This system can be realized with noninteger order integrators as in Fig. 1. The conditions for stability of filter (6) are presented in Theorem 1. (6) or (7)- (8) for α ∈ (0, 1) is asymptotically stable if and only if one of the following conditions holds

Theorem 1 System in the form
Proof Proof can be found in [29].

Discrete Realization of Non-integer Order Filters
In order to implement filters digitally, one needs to take under consideration that non-integer order systems have infinite memory. Their natural discrete realization in the form of Grünwald-Letnikov derivative is not available as it would fill any computer's memory. Grünwald-Letnikov derivative can be, however, used for filtering of finite length signals, e.g., in post-processing. It can be done using numerical method proposed in [28] based on Grünwald-Letnikov derivative (3) in the following form where u(t) is the original signal and u f (t) is the filtered signal. This direct approach will be used in this paper as a comparison with proposed approximations which are presented in the following sections.

Oustaloup Approximation
In this section, we will present well-known Oustaloup approximation, problems with its direct discretization and proposed time-domain realization.

Continuous Oustaloup Approximation
Continuous Oustaloup filter approximation to a fractional-order differentiator G(s) = s α is widely used in applications [25]. An Oustaloup filter can be designed as where: Approximation is designed for frequencies ω ∈ [ω b , ω h ], and N denotes the order of the approximation. As it can be seen, its representation takes form of a product of a series of stable first-order linear systems. As one can observe, choosing a wide band of approximation results in large ω u and high-order N results in spacing of poles from close to −ω h to those very close to −ω b . This spacing is logarithmic with a grouping near −ω b and causes problems in discretization. Wide band of approximation is desirable, because approximation behaves the best in the interior of the interval and not at its boundary, so certain margins need to be kept.

Direct Discretization of Transfer Function
In applications, especially when using MATLAB/Simulink environment for real-time control or filtering, two methods of direct implementation of continuous-time systems are established: 1. Use them as continuous blocks in Simulink, then use a fixed step explicit solver (e.g., RK4) to evaluate blocks in real time. 2. Create discrete transfer functions using some discretization scheme-usually Tustin because it should preserve stability.
The problem with first approach that uses explicit solvers for Oustaloup filters comes from the fastest pole (see [4]) of the transfer function. If the band of approximation is wide, and includes high frequencies, this pole can be located very far to the left of complex plane. It can be observed (see [18]) that Runge-Kutta-type algorithms preserve stability of linear differential equation for such eigenvalues λ and discretization step T that stability function of the algorithm R(T λ) is less than one. For classical Runge-Kutta algorithm, the stability function is That means that for λ ∈ R, T λ must belong to interval (−2.785293..., 0). In considered case and wide bands, the sampling frequency must be at least relatively close to the upper end of the approximation band, which can lead to errors in approximation of frequency characteristics near the Nyquist frequency. Problems with the second approach come from different reasons. Poles of continuous transfer function have a tendency to group near 0 (especially when lower frequency is small). Those poles will be mapped close to 1 during discretization. In that case, when discretizing every pole separately, denominator of the entire system will include a group of discrete poles similar to where ε i > 0 are distances of the pole from stability boundary (point (−1, j0) on complex plane) and they are usually very small (orders of magnitude from 10 −4 to 10 −9 are not uncommon). In that case, final denominator will include numbers that will be products of ε i with each other, resulting in numbers close to, or below 2.22 · 10 −16 which is a smallest number that can be added to another in MATLAB (so-called machine epsilon; in different computation systems, the number is of similar magnitude). It results in rounding error, which high sensitivity of polynomial roots to coefficient values leads to instability. These rounding errors are unavoidable, even when substituting to analytical formulas.
This phenomenon can be observed for example for approximation of order N = 20. Rounding errors cause Tustin approximation to map real stable poles into complex unstable poles (see Fig. 2).

Time-Domain Realization
This method was developed in order to counteract the rounding errors present when discretizing the Oustaloup transfer function.
One can easily observe that for zero initial condition Because of that (13) can be written in vector matrix notatioṅ or in briefẋ = Ax + Bu What can be immediately observed is that the matrix A is triangular. This is extremely important in this problem, as all its eigenvalues (poles of transfer function) are on its diagonal, so there is no need for computing products of eigenvalues, which would lead to rounding errors. Desired property of discretization would be that it should preserve triangular structure, so there would not be any need for eigenvalue multiplication.

State-Space Discretization
Desired discretization should preserve both traingular structure and stability of the system. Explicit numerical methods generally preserve structure (Euler forward, RK4), but they are still susceptible to problems listed earlier. That is why implicit methods have to be used. In this case, we have used state-space Tustin approximation [15]. Discretized system takes form (for simplicity, e.g., u(T k) is written as u(k)) where As it can be observed, in this case the matrix is a product of two lower triangular matrices, which preserves structure. Thus, the diagonal consists of discrete system eigenvalues which are discretized individually, so the field for rounding errors is severely reduced. It can be observed that all eigenvalues are in R and the system is controllable. Also because the method is implicit, the eigenvalues far from imaginary axis are mapped well into unit circle.

Remark 2
It should be noted that Tustin method is implemented in MATLAB c2d function; however, this implementation is not optimized for detecting matrix structure and resulting matrix is not fully triangular (there are nonzero elements with order of magnitude 10 −12 over diagonal). Implementation using efficient matrix inversion algorithms with structure detection solved the problem.

Filter Realization
The realization of considered filter with use of this variant of Oustaloup method is depicted in Fig. 3. This realization is difficult to analyze, as one has difficulty finding poles of closed-loop system. The closed-loop equivalent system in matrix notation is presented in Lemma 1.

Lemma 1 Time-domain discrete approximation of filter
is given by Fig. 3 Realization of filter (6) using discrete time-domain Oustaloup approximation where Proof The proof is tedious and not very innovative. That is why it is omitted. It requires a solution of a simple set of linear equations and observing which products of vectors result in scalars.

Laguerre Impulse Response Approximation Method
The method of Laguerre impulse response approximation LIRA was introduced by authors in [5]. It should not be confused with traditional discretization method. It is a method for approximating systems given by transfer functionŝ or more precisely defined by convolution where j ≤ σ j ≤ j +1, j = 1, 2, . . . , n, j ≤ γ j ≤ j +1, j = 1, 2, . . . , m, p j ,q j ∈ R.
The initial conditions are zero. It is also assumed that |u(t)| ≤ u max for t ≥ 0 and u(t) = 0 for t < 0 [5]. We assume thatĝ is a Laplace transform of a certain function g : [0, ∞) → R which fulfills g ∈ L 1 (0, ∞) ∩ L 2 (0, ∞). This function g is called the impulse response.
It can be shown that the convolution (28) for given input signal u can be approximated with a solution of a system of n linear ordinary differential equations. The approximation uses an orthonormal basis in L 2 (0, ∞) e k (θ, μ) = 2μe −μθ L k (2μθ ), k = 0, 1, 2, . . . (29) where μ is an arbitrary positive constant and L k is k-th Laguerre polynomial of form Theorem 2 gives the conditions that must be fulfilled in order to find the approximation with minimal error.
Theorem 2 If g ∈ L 1 (0, ∞) ∩ L 2 (0, ∞) and |u(t)| ≤ u max then: 1. Response of system defined by convolution (28) can be approximated with where set of functions ξ k (t) : [0, ∞) → R is solution of a system of linear differential equationsξ and 2. For every ε > 0, there exists a number n 0 dependent on g, ε and u max for which approximation error e n (t) = x(t) − x n (t) fulfills the inequality for all n ≥ n 0 and t ≥ 0 Proof For the proof, see [5].
The formula (33) for calculating the coefficients is not convenient for numerical implementation as it requires an analytical fomula for impulse response. In [5], the authors presented the following recurrence formula where Remark 3 Choice of parameter μ Method is convergent for any set of orthonormal Laguerre functions, parametrized by μ. However, performance of the method strongly depends on its value, especially for lower orders. In authors' earlier works, it was shown that μ chosen by maximization of function: introduces smallest approximation error in the sense of L 2 (0, ∞). Influence of choosing nonoptimal μ was investigated in [8,47].
In summary, the approximation of system defined by convolution (28) can be given in matrix notation as (e.g., for approximation order n + 1, ξ = ξ 0 ξ 1 ξ 2 ξ 3 ξ 4 . . . ξ n T ): As it can be seen, the state matrix in (38) is also lower triangular. Because of that use of Tustin discretization scheme is reasonable, as it preserves this structure. It should be, however, noted that similarities between time-domain discrete Oustaloup method and LIRA end here.
Instead of approximating the integrators, LIRA method approximates the entire system. As it can be seen in Lemma 1, the closed-loop realization of filter (6) with Oustaloup integrators is quite complicated and closed-loop eigenvalues are not easily available for analysis. In case of LIRA, eigenvalues of filter are not only available, but also state matrix in (38) has n identical eigenvalues and their value can be regulated (moved away from right half plane) and their discrete form is naturally more robust toward rounding errors. It should be also noted that this approximation leads to effective approximation of any order, whereas Oustaloup approximation allows only even values (because two integrators have to be realized).
It should be also noted that explicit discretization schemes can be used, as there are no eigenvalues located farther to the left of imaginary axis than others, which is a Fig. 4 Choice of μ maximizing J (μ) for filter (6) and different values of ξ and α and ω 0 = 1, darker colors indicate the lower values natural cause of numerical instability. For further analysis, the same code for obtaining Tustin discretization was used as for the time-domain Oustaloup method.
One weakness of the method is the inability to easily determine the set of filter parameters for which the impulse response of filter fulfills the assumption g ∈ L 1 (0, ∞) ∩ L 2 (0, ∞). Experiments show that assumptions of appropriate integrability are fulfilled if the highest order of denominator of corresponding transfer function is at least 1/2 and the transfer function is stable. Fortunately, one can easily observe that the method does not converge for systems not fulfilling that assumption. In Fig. 4, the optimal values of parameter μ are computed for filter (6) assuming that ω 0 = 1 and other parameters are changed in viable range. In two regions, the convergence was not obtained. First one is the area of unstable pairs (ξ, α) (see theorem 1), and the other one is located in the area of very low orders α. What is interesting is that in the considered range of convergence all μ are of the same order of magnitude (around 1) and the area boundary is not smooth.
Other observed weakness of the method is visible for very high orders of approximation (over 20). In that situation, the factorials present in the formulas for β coefficients cause numerical errors of rising magnitude. More on that effect is given in [8,47].

Analysis of Filter Behavior
In this section, LIRA and time-domain discrete Oustaloup approximations will be analyzed in order to show merits and flaws of both methods. Whithout loss of generality, analysis is conducted for filter (6) with ω 0 = 1 and varying values of ξ ∈ (−1, 1) and α ∈ (0, 1). Sampling period was chosen as 1ms. For the Oustaloup method, the frequency band was chosen as [10 −6 , 10 3 ] and in that band the characteristics were compared. For LIRA method parameter, μ was chosen as maximum of (37) for every case.
The comparison will include: -analysis of approximation error expressed in the form of H ∞ norm for wide spectrum of parameters; -analysis of Bode plots for selected filters; -filtration of an EEG signal.

H ∞ comparison
This analysis was conducted in order to find an objective criterion for determining the quality of approximation. Because the filter realization is the main concern of this work, H ∞ norm was chosen as it represents the frequency characteristics. In Tables  1 and 2, the values of approximation error norm G(s) −G(s) ∞ are presented for approximation orders of 6 and 12, respectively, for both methods of approximation. G(s) denotes the approximation of filter. The values of H ∞ norm were calculated numerically by computing frequency response of approximations state-space representations in MATLAB (for considered frequency band) and subtracting it from analytically computed frequency response of original filter. Then, maximum of modulus was chosen. In both tables, columns correspond to values of α and rows correspond to ξ . For row of ξ , there are two sub-rows containing errors of time-domain Oustaloup    Tables 1 and 2 allows drawing the following conclusions: 1. For values of α ≥ 0.5, low-order approximations obtained by LIRA method were either better or at least comparable to those obtained by time-domain Oustaloup method. For higher orders, approximation errors are still comparable. 2. For α < 0.5 LIRA approximations are generally worse, except for some isolated cases. This behavior of LIRA is justified that only for those values of α the impulse responses are believed to be in the required space of L 1 ∩L 2 . It should be, however, noted that it was not yet formally proven. 3. Generally, for LIRA approximation, errors are decreasing with its order, but convergence (in the sense of H ∞ norm) is slow. Therefore, it can be observed that if approximation is correct it will be correct even for low orders. 4. Increasing order of Oustaloup approximation results in more evenly decreasing error for all admissible values of α. Remark 4 It should be noted that behavior described in the above conclusions is consistent for all investigated orders of approximation. We have included only sixth and 12th orders because all the important effects are visible.

Detailed Analysis of Selected Filters
In this section, we describe in detail frequency responses of selected filters and their approximations. They were chosen to represent significant effects visible in Tables 1 and 2. The first of the considered examples is the filter with α = 0.7 and ξ = 0, i.e., Results of approximation of sixth order for this filter are presented in the Fig. 5. As it can be seen consistently with conclusion of previous section, the LIRA method provides very good approximation, especially in the magnitude aspect. Also phase approximation is consistent for frequencies ω < 10ω 0 . In case of the Oustaloup approximation, oscillations in both magnitude and phase characteristics are observed. Substantial differences are especially visible in the phase characteristic for ω > ω 0 . Approximations of the same filter of 12th order are presented in the Fig. 6. As it can be seen increasing of order in LIRA leads to very similar results. Slight differences are only in the phase characteristic. Increasing order in the Oustaloup method leads to substantial improvement, as the oscillations are strongly reduced and phase characteristic is consistent for ω < 10ω 0 .
The second example considers the case of α < 0.5, in particular filter with parameters α = 0.3 and ξ = 0.6, i.e., Results of approximation of sixth order for this filter are presented in the Fig. 7. Approximation behavior is, as expected in previous section, in favor of Oustaloup method. Magnitude characteristic is oscillating around the ideal value visibly keeping the trend. Phase characteristic exhibits oscillatory behavior. Approximation obtained with LIRA method is consistent only for frequencies ω > 10ω 0 . For lower frequencies, it provides static damping close to 16 dB. Phase characteristic also is significantly different from the desired one. Increase in order to 12th leads to results presented in the Fig. 8. For Oustaloup method, a substantial improvement is observed, especially in the magnitude characteristic, which now very closely follows the characteristic of (40). In case of LIRA, the It is one of the mentioned before isolated cases, where for α < 0.5 the LIRA method provides substantially better low-order approximation than Oustaloup method. It can be observed in Fig. 9. Approximation behavior is very similar to the one characteristic for α > 0.5. It should be noted that increase in order for Oustaloup method improves the quality of approximation (not presented in detail, but visible in Table 2).

Example of Filtering
In this section, a brief example of filter behavior in time domain will be presented. A set of EEG data is being considered for analysis of delta waves which are located in the frequencies below 3 Hz. Authors were tasked with providing a non-integer low-pass filter which would introduce slower damping than 40 dB/dec with standard second-order filter [20,21]. Filter (6) was used with parameters: α = 0.7, ξ = 0.5, ω 0 = (6π) α . Signals were sampled with frequency of 50 Hz. This filter was originally used in postprocessing with non-integer filtration given by (9)- (12). Here we present it as a comparison with approximated filters. Filter approximations were of sixth order. As one can see in the Fig. 10, the results are qualitatively identical. Analysis of frequency spectrum of filtered and original signals allows to observe the differences between realizations (see Fig. 11). Methods proposed in the paper realize filtering properly offering the expected damping of 28 dB/dec. The numerical method given by (9)-(12) is presenting inferior behavior. Filtration realized in that way introduces unwanted damping for low frequencies, and for frequencies above 3 Hz, the damping is not consistent with expected 28 dB/dec.

Conclusions and Future Work
In this paper, two efficient methods for digital realization of non-integer order filters were compared. Both these methods have their own merits. Discrete time-domain Oustaloup method allows the approximation of the filter for wider spectrum of parameters than LIRA. However, its coefficients are more numerically sensitive to discretization. Also, the necessity of using two integrators leads to complicated matrices as presented in Lemma 1. Also influence of parameters on the eigenvalues of closed-loop systems is not straightforward. These weaknesses do not appear in LIRA method. The structure of state-space representation is simple, and behavior of μ is predictable. On the other hand, verification of the assumptions is difficult. Regarding the quality of approximation, there is no clear 'winner' as the performance strongly depends on the order of approximation and parameters of the system. Impulse response approximation leads to filters with smaller gains for low frequencies. However, for low order of approximation, LIRA is better (or at least comparable) than Oustaloup method. Oustaloup-based method usually works in the desired bandwidth, sometimes with problems on the boundary (see, e.g., Fig. 6). Future work will focus on the following areas. First of them is the improvement of LIRA method in aspects of assumption verification and numerical errors for approximations of order over 20. The second aspect is the verification if the realization of Oustaloup transfer function is optimal. In authors' opinion, it is superior to the currently discussed in the literature (one of them has diagonal state matrix, which should be beneficial, but coefficients in the control matrix are very sensitive to rounding errors and another one is a Frobenius matrix realization, which has the same sensitivity as the transfer function). The third considered area is the real-time implementation of the filters. Authors have experience in implementing simple non-integer order filters on embedded platforms and want to develop it further (see, e.g., [14]).