Fractional-Order System Identification and Equalization in Massive MIMO Systems

,


I. INTRODUCTION
Digital communication using massive multiple transmit and multiple receive antennas also known as massive Multiple-Input Multiple-Output (MIMO) has been one of the most important technical developments in modern day communication. In massive MIMO systems, the base station (BS) is equipped with hundreds of antennas to service several mobile user equipment (UE), where the number of base station antennas, is significantly larger than the number of user terminals. Massive MIMO is one of the key enablers for fifth generation (5G) wireless networks [1], and 5G is expected to support large amounts of data, allow the connection of a large number of wireless devices, improve coverage, quality of The associate editor coordinating the review of this manuscript and approving it for publication was Hayder Al-Hraishaw . service (QoS), latency, reliability, security, energy efficiency and spectrum efficiency [2] and [3].
To fully achieve the benefits of massive MIMO, we must overcome several challenges, which include pilot contamination, radio frequency (RF) impairments, channel estimation, etc. Accurate channel state information (CSI) is very important in massive MIMO systems as it has a great impact on the performance of the system, but in reality, channel state information is not readily available. Thus, for improved performance, the massive MIMO wireless channel must be efficiently estimated, and this brings about the idea of channel estimation which in this paper we will refer to as system identification.
System identification is a technique used to develop mathematical models based on the input and output data sets of a system to represent the characteristics or dynamics of VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ that system. It is widely used in many areas such as broadcasting theory, geology, hydrology, communications, control, etc [4]. There is need for exploring novel system identification methods or improving on the already existing system identification techniques especially in the field of wireless communications. The finite impulse response (FIR) and state-space system identification models are models that have been widely used to identify MIMO wireless communication channels and have been proven to provide good model approximations or good channel estimation. State-space identification models provide better model approximation compared to FIR identification models especially when the MIMO channel is frequencyselective [5]. A state-space model is a mathematical representation or modelling of a physical system as a set of input and output data and state variables related by a first-order differential equation [6].
State-space models have given rise to subspace system identification (SSI) algorithms, namely Multiple-Input Multiple-Output Output-Error State Space (MOESP), Numerical Algorithm for Subspace State-Space System Identification (N4SID) and Canonical Variate Analysis (CVA). Subspace system algorithms are algorithms that can estimate a linear time-invariant state-space model, directly from the input-output data [7]. The advantage of subspace system identification (SSI) algorithms is that they can identify a system in a straightforward way using numerically robust computation tools such as linear quadratic (LQ) decomposition (representation of a matrix in a simpler form via orthogonal transformations) and singular value decomposition (SVD). The focus of this paper will be on the MOESP algorithm. The classical MOESP system identification algorithm has given rise to the fractional-order MOESP system identification algorithm. Fractional-order system identification is one technique that has been identified to help improve the identification of any given system. This builds from the integer-order system identification technique and has been mostly applied in the field of control. Fractional-order calculus is very important, especially in explaining many events which traditional mathematics cannot explain [8]. It has also been observed that many real-world physical systems are well characterized by fractional-order differential equations rather than using classical integer-order models [9]. In signal processing, fractional-order operators are used in the design of fractional-order integrators and fractional-order differentiators and for modelling speech signals [10].
The fractional-order time-derivatives of the input-output data are generally not measured meaning that the input-output matrices are not known, as a result the classical subspace methods originally developed for the identification of discrete-time models cannot be directly adapted for the identification of continuous-time fractional-order models. To address this problem the Poisson moment functional (PMF) approach is used when dealing with continuous-time fractional-order system identification, where the input-output data is first filtered using the PMF filter after which the MOESP identification algorithm outlined for discrete-time system identification is then applied to the PMF filtered input-output data. Details of applying the PMF to fractional-order system identification are discussed in Section IV.
Since analytical solutions of fractional-order differentiations and integrals are complicated, the authors in [11] proposed more simplified methods when dealing with fractional-order calculus. The commonly used definitions to solve fractional-order differentiations and integral are the Riemann-Liouville, Grünwald-Letnikov and Caputo definitions, details of which are presented in Appendix A. The Riemann-Liouville and Grünwald-Letnikov definitions are more suitable in describing fractional-order calculus problems with zero-initial conditions whilst the Caputo definition is useful in discussing systems with nonzero initial conditions [12]. In this paper we use the Riemann-Liouville definition when dealing with fractional-order system identification. The geometric interpretation of fractional-order calculus using the Riemann-Liouville definition is presented in Appendix B.
The rest of the paper is organised as follows. In Section II, we present the related work and paper contribution, and the massive MIMO system model is presented in Section III. Section IV discusses on system identification using the fractional-order MOESP algorithm. Channel equalization using the fractional-order model is discussed in Section V. The discussion on symbol error rate is presented in Section VI. Our contribution is evaluated in Section VII, and we conclude our paper in Section VIII.

II. RELATED WORK AND PAPER CONTRIBUTION
The use of state-space models in system identification in communications has been motivated by the works by Li [13] and Zhang and Bitmead [5]. In [13], Li used both discrete and continuous-time state-space techniques to estimate flat and frequency-selective fading wireless communication channels. In [5], the authors used the discrete-time integerorder MOESP algorithm, a subspace system identification algorithm to identify MIMO frequency-selective wireless channels.
To further improve on the identification of the wireless communication channel, this work proposes the use of the fractional-order identification model. The major advantages of fractional-order calculus are: (i) They make it possible to obtain mathematical models that describe the results closer to the experimental measurements. (ii) Their ability to predict more accurately the dynamics of the system that is being modelled. (iii) They make it possible to obtain simplified system models with just a few physically motivated parameters [14] and [15]. The main contribution of this paper is on the use of the continuous-time fractional-order MOESP system identification algorithm to help identify massive MIMO frequency-selective wireless channels. The use of the fractional-order modelling approach in identifying the massive MIMO frequency-selective wireless channels is motivated by its advantage of smaller residual errors for the identified model compared to the integer-order model, and this further enables the adoption of a fractional-order equalization approach. On the other hand, subspace system identification algorithms have been proven to successfully identify the dynamics of a system and this has motivated their application in this paper [16].
The use of the continuous-time fractional-order MOESP system identification algorithm in the identification of massive MIMO frequency-selective wireless channels provides a novel research perspective in the area of wireless communications.

III. MASSIVE MIMO SYSTEM MODEL
We consider a massive MIMO wireless system as shown in Fig. 1, with a BS equipped with m transmitting antenna elements, and a terminal station equipped with p receiving antenna elements having length-L ISI channel paths, where for massive MIMO, m p, u (m) is transmitting antenna m, y (p) is receiving antenna p, h (m,p) is the channel path between transmit antenna m and receive antenna p.
We assume that the channel is quasi-static, i.e., it is time-invariant within each frame and changes independently from frame to frame. In addition, the antennas are sufficiently spaced such that there is no interference between antennas. Training symbols known to both the transmitter and receiver are inserted at the start of each frame to assist with channel estimation. The received signal of the linear-time invariant length-L ISI channel massive MIMO system is expressed as where y (t) is the p × 1 received signal vector, u (t − l) is the m × 1 transmitted symbols vector at time (t − l), n (t) ∼ CN 0, N 0 I p is the p×1 additive white Gaussian noise vector at the receiver side at time t, with N 0 being the noise power and I p is the p × p identity matrix and H l is the p × ml th path massive MIMO channel matrix coefficients.

IV. SYSTEM IDENTIFICATION USING FRACTIONAL-ORDER MOESP ALGORITHM
Fractional-order state-space models in controllable, observable and diagonal canonical forms are similar to integer-order state-space models [17]. A fractional-order linear time invariant (LTI) system is mathematically equivalent to an infinite dimensional LTI filter. Thus, a fractional-order system can be approximated using higher-order polynomials having integer-order differ-integration operators.
The following assumptions are necessary when dealing with system identification: i) The system is persistently excited by the training symbols. ii) The system is stable, observable and controllable.
iii) The dimensions of matrix A are known, and rank (A) = n, where n is the order of the system and A is the system matrix. The rank is the number of linearly independent rows. iv) The random white noise is uncorrelated to the input signal. We consider commensurate fractional systems, i.e. where all the differentiation orders are multiple integers of α [18]. The dynamics of the massive multiple-input multipleoutput linear time invariant (MIMO LTI) system can be modelled using fractional-order state-space modelling. Ignoring the effects of the process noise and additive noise, and using continuous-time fractional-order state-space modelling, the linear time invariant length-L ISI channel massive MIMO system, (1), can be expressed as in [17].
is the p × 1 output vector, A is the n × n system matrix and it describes the dynamics of the system, i.e. the eigenvalues of the system, B is the n × m input matrix and it describes the linear transformation by which the inputs influence the next state, C is the p × n output matrix, and it describes how the state is transferred to the output, D is the p × m feed-forward matrix,αis the commensurate fractional-order and D α is the fractional derivative of order α. Equation (2a) is called the fractional-order state equation and (2b) called the output equation.

Y(s) = CX(s) + DU(s)
Rearranging (4), the transfer function of the system can be expressed as where I n is an n × n identity matrix. The general characteristics of the transfer function, (5) are [19] (i) The magnitude curve has a constant slope of −20α dB/decade. (ii) The phase plot is a horizontal line of value −απ/2. (iii) The gain crossover frequency depends on A.
(iv) The Nyquist plot is a straight line which starts from the origin with argument −απ/2. It is important to know if a fractional-order system is stable or not and in the fractional-order LTI case, the stability is different from that of the integer-order one. In the fractionalorder case a stable fractional-order system may have roots on the right half of the complex plane as shown in Fig. 2. A fractional-order system is stable if [12] and [20], as in where 0 < α < 2, eig (A) gives the eigenvalues of matrix A and λ k is the k th eigenvalue of A and −π < arg (λ k ) ≤ π . Since the continuous-time state-space representation of commensurate fractional-order systems is similar to that of integer-order systems [21], the analysis of the fractional-order MOESP model follows the one outlined for the classical integer-order MOESP model as proposed in [7] and [17].

A. CONSTRUCTING MASSIVE MIMO DATA MATRICES FOR CHANNEL IDENTIFICATION
In constructing the data matrices for the continuous-time fractional-order system, (2) and the MOESP algorithm are considered. By computing the successive α-order fractional derivatives of (2) till the (i − 1) α-order derivative, the data equations can be constructed.
Stacking (2b) and all its α-order derivatives into a column vector, the input-output relationship can be written as [22] . . .
where i is a user defined index. In compact form, (8), can be written as where where the notation (·) 0|(i−1)α denotes the zero th order derivative to the (i − 1) α th order derivative, i and i are the observability and Toeplitz matrices, respectively.

B. ESTIMATING THE EXTENDED OBSERVABILITY MATRIX USING THE POISSON MOMENT FUNCTIONAL FILTERED DATA
With reference to (10) and (11), the data matrix, W 0|(i−1)α,N can be constructed as Equation (14) contains successive α-order fractional derivatives of the input-output data which in most practical cases are not measured. To address this problem, we propose the use of the Poisson Moment Functional (PMF) approach, details of which are presented in Appendix C. Applying the PMF to the data matrix, W 0|(i−1)α,N , the PMF filtered data can be expressed as Having filtered the input-output data, the MOESP identification algorithm can then be applied. Applying linear quadratic (LQ) decomposition to matrix (15) yields Using (16) compute the SVD of L 22 as where S 1 is a diagonal matrix containing the singular values different from zero and the order of the system, n can be estimated from S 1 as the number of singular values different from zero. The extended observability matrix, i , can then be estimated from (17) The estimates of A, B, C and D can then be obtained directly fromˆ i . Matrix C can be obtained directly from the first p rows ofˆ i asĈ Matrix A can be obtained using the shift property as where i isˆ i with the last p rows removed and¯ i isˆ i with the first p rows removed. The estimate of A can also be expressed aŝ Once the estimates of A and C are known, the estimates of B and D can be obtained using the following input-output equation where ˆ i is the estimated extended observability matrix written aŝ whereˆ ⊥ i is a full row rank matrix (i.e. all the rows and columns are linearly independent) satisfyingˆ (25) can be rewritten as Equation (26) can then be written as which can be solved using least squares.

C. ANALYTICAL SOLUTION OF THE FRACTIONAL-ORDER MASSIVE MIMO SYSTEM
Considering (2a), and applying the fractional-order integration of order α on both sides results in where I α is the fractional-order integral operator and x (0) is the initial condition. Assuming that the general solution of (28) can be written as [24] x where and where x k−1 (t) is the previous state vector. Using (30) and (31) the state vector can be recursively written as follows At k = 1 Substituting (30) into (32) yields At k = 2 Substituting (33) into (34) yields At k = 3 Substituting (35) into (36) yields and so on. With reference to the above equations, the general expression for x k (t) can be written as Using the property [24] equations (33), (35), (37), etc. can be written as respectively, where γ is taken to be γ = 0 since we do not have the t term. Using the above property, the general expression for x k (t) can be rewritten as Using (43), equation (29) can now be rewritten as In terms of the Mittag-Leffler function in [25] and [26], equation (44) can be written as Substituting (45) into (2b), the general solution of the massive MIMO system having used fractional-order system identification can be written as

V. CHANNEL EQUALIZATION USING THE FRACTIONAL-ORDER MODEL
The equalization technique we employ here follows the works by Zhang and Bitmead [27], Zhang and Bitmead [28] and Al-Dhahir and Sayed [29] where they proposed the use of the minimum mean square error -decision feedback equalizer (MMSE-DFE) to reduce the effects of ISI in a multiple-input multiple-output (MIMO) system. The application of the MMSE-DFE in the equalization of the fractional-order massive MIMO channel model is a novelty in wireless communications. The choice of the MMSE-DFE is that it benefits from both the advantages of the minimum mean square error (MMSE) and decision feedback (DF) equalizers.
Since the continuous-time state-space representation of commensurate fractional-order systems is the same as that of integer-order systems [21] and the analysis of the fractionalorder MOESP model follows that of the classical MOESP model as proposed in [7] and [17], we conclude that the MMSE-DFE for the commensurate fractional-order system will also follow the procedure outlined in [27] for the MMSE-DFE for the state-space model. Based on this, the equalization of the fractional-order channel model is presented in the following section.
Having obtained matrices A, B, C and D using the fractional-order system identification algorithm, these matrices can then be used for massive MIMO channel equalization in the context of fractional-order channel modelling. Following the equalization procedure outline in [27] which uses the MMSE-DFE to recover the transmitted signal, a length N f block of the input-output data, i.e. the blockwise data model for the commensurate fractional-order model (2) can be expressed as (47), which is shown on the bottom of the next page, where N f is the length of the feedforward filter matrix or the number of feedforward filter taps and is the equalizer's decision delay. The decision delay helps in determining which symbol is detected at the current time, t. A more compact representation of (47) can be written as where Y FO (t) is the effective fractional-order output vector, U FO (t) is the effective fractional-order input vector, x t + − N f + 1 is the state vector, N FO (t) is the effective fractional-order noise vector, is the observability matrix and is the Toeplitz matrix.
With reference to [27] the blockwise data model (47) can be decomposed in terms of the impulse response as (49), which is shown on the bottom of the 9th page. Since the term x t + − N f + 1 is a contribution from distant past symbols, it can be subtracted from the effective fractionalorder output vector, Y FO (t) resulting in This signal is then passed to the feedforward filter which helps in minimising the effects of ISI and noise from future symbols. The signal at the output of the feedforward filter is expressed as F H Y FO (t), where F represents the feedforward filter coefficients. Signal F H Y FO (t) is then passed through the detector to obtain an estimate of the originally transmitted signal, which is expressed asû (t) FO−MMSE−DFE . After detection, the detected symbols,û (t) FO−MMSE−DFE , are fed back to the DFE through the feedback filter, B. The fed back symbols are then subtracted from the incoming symbols to help minimise the effects of the previously detected symbols from the incoming signal, Y FO (t).
Next, we outline the procedure to follow when calculating the feedforward and feedback filter coefficients. We start by defining matrix W which is a combination of the feedforward and feedback filter coefficient matrices expressed as Expressing W as a combination of the feedforward and feedback filter coefficients allows for the joint calculation of these matrix coefficients, thus circumventing the need to first calculate the feedforward filter coefficients and then using these to calculate the feedback filter coefficients and viceversa. A term Z FO (t) which is a combination of the effective fractional-order output vector, Y FO (t) and the previously detected symbols, D (t) can also be defined as . . .
is a vector of previously detected symbols. Assuming correct detection, these are the same as the transmitted symbols, and N b is the length of the feedback filter. Matrix Z FO (t) plays a great role in the joint calculation of the feedforward and feedback filter coefficients. Using (54), the auto-correlation matrix, R Z FO Z FO of Z FO (t) can be defined as Substituting (54) into (56) results in Using (52) in (57) yields where P s is the total transmit power which is assumed to be evenly distributed between all the m transmit antenna elements. The cross-correlation matrix, R Z FO u of Z FO (t) and u (t) is defined as Using (54) and then (52) in (59) yields According to [27] and using (58) and (60), the MMSE solution of the MMSE-DFE can be rewritten as Using (54) and (61), the estimate of the transmitted symbols at the output of the MMSE-DFE,û (t) FO−MMSE−DFE having considered fractional-order channel modelling can be expressed aŝ

VI. SYMBOL ERROR RATE PERFORMANCE
According to [30] and [31], the output SNR of the unbiased MMSE-DFE can be expressed as and the output SNR of the biased MMSE-DFE can be written as where R ee,min is the covariance matrix of the error vector which can be expressed as With reference to [32] and using the Moment Generation Function (MGF) approach, the average symbol error 86488 VOLUME 8, 2020 rate (SER) having considered M-ary quadrature amplitude modulation (M-QAM) for the frequency-selective quasi-static Rayleigh fading channel can be expressed as where M is the number of constellation points. and is the MGF of the SNR, γ for Rayleigh fading. The MGF approach is considered on the bases that it is a bit easier to work with as it avoids the need to find the probability density function (PDF) of the output SNR which can be a bit tedious for high-order systems [33]. Substituting (68) into (66) and after some manipulations, equation (66) can then be expressed as [34] where a = Replacing γ with SNR MMSE−DFE in (69), the average SER for M-QAM having considered the MMSE-DFE signals can then be expressed as

VII. SIMULATION RESULTS
For our simulation, we used MATLAB having considered the input signal, u (t), to be a chirp signal with frequency ranging from 20Hz to 20 kHz with a sampling frequency twice the highest frequency range, which is a persistently excitation signal. This frequency range was chosen to cover the range of audible frequencies for humans. The fractionalorder, α, was chosen to be α = 0.1. This is because for fractional-orders greater than α = 0.1 the system proved to be unstable, that is, the conditions in (6) and (7) fail. Given that we assumed zero initial conditions for our fractional-order model, the Riemann-Liouville definition was chosen for our simulations as it is more suitable in describing fractionalorder calculus problems with zero-initial conditions [12]. The fractional-order MOESP system identification algorithm was based on the parameters in Table 1.
In system identification, the performance of the model depends on how much data is available to describe the system and whether the sampling time is good enough to capture the system dynamics. In our simulations, the sampling time, 1s was identified as the sampling time that could reliably capture the system dynamics. The user defined index, was selected to be i = 3 as for lower values of i the dynamics of the massive MIMO system could not be well captured. The selection of the PMF filter order, r, Poisson filter gain, β, and Poisson filter constant, λ, was based on the criteria defined in the Appendix C, where for simplicity we assumed that, λ = β = 1. Fig. 3(a) shows the fractional-order identification results compared to the integer-order results, and (b) shows the zoomed part of the results to give a clear picture. These results show that the fractional-order MOESP algorithm is able to reliably capture the dynamics of the massive MIMO system when the appropriate fractional-order, α value is selected. Comparing the output of the fractional-order MOESP identified system and the measured output/ actual output, gives a mean square error (MSE) of 0.0455. This is better than the MSE of 0.4237 for the same massive MIMO system identified using the integer-order MOESP algorithm of order n = 7. Order n = 7 being the best performing system order when applying the integer-order MOESP algorithm to our massive MIMO system.
To study the performance of the MMSE-DFE for the fractional-order massive MIMO channel model, simulations were run in MATLAB based on the parameters in Table 2.
The choice of these parameters was based on [27], owing to the fact that the continuous-time state-space representation of commensurate fractional-order systems is the same as that  of integer-order systems [21]. In [27], they stated that given the number of channel paths, L, the decision delay of the equalizer is given as = L+2, the optimal feedforward filter length is, N f = + 1 and the optimal feedback filter length is, N b = L. The fractional-order MOESP identified system with fractional-order of order α = 0.1 was considered. Fig. 4 shows the performance of the MMSE-DFE for the fractional-order massive MIMO channel model for the optimal feedforward filter length of N f = 6 and the optimal feedback filter length of N b = 3 compared to the originally transmitted signal in each transmitting antenna for fractional-order α = 0.1.
We can observe that the MMSE-DFE for the fractionalorder massive MIMO channel model is capable of carrying out the desired equalization to recover the transmitted signal.  The mean square error (MSE) for this equalization model is 0.15273387. This compares to the MMSE-DFE for the integer-order massive MIMO channel model of order n = 7. Whilst the MSE for the MMSE-DFE for the FIR equalization model is 0.15273753. These MSE results show that a high order integer-order model is needed to achieve similar results as the fractional-order model.
Next, we present the simulation results for the symbol error rate performance for the fractional-order massive MIMO channel model in comparison to the integer-order and FIR channel models. For our simulations we considered BPSK, QPSK and 256-QAM modulated signals. Fig. 5 to Fig.7 show the SER performances for the FIR, integer-order and fractional-order massive MIMO channel models having applied MMSE-DFE plotted against the signal-to-noise ratio for BPSK, QPSK and 256-QAM modulated signals.
We can observe from Fig.5 to Fig. 7 that the MMSE-DFE for the fractional-order massive MIMO channel model has lower symbol error rate results compared to the MMSE-DFE for the FIR massive MIMO channel model. This implies that the fractional-order MMSE-DFE model can be a more robust choice than the FIR MMSE-DFE model when it comes to massive MIMO channel equalization. It can also be observed that the SER performance of the MMSE-DFE integer-order state-space channel model of order n = 7 compares to that of the fractional-order channel model of fractional-order α = 0.1. This proves the viability of the fractional-order state space models in channel equalization.

VIII. CONCLUSION
We presented the continuous-time MOESP system identification and equalization algorithms for the massive multipleinput multiple-output (MIMO) frequency-selective wireless channels having considered the fractional-order channel model. The results showed that the proposed continuous-time fractional-order MOESP system identification algorithm can actually be used to identify the dynamics of the massive MIMO system when the appropriate fractional-order, α value is selected.
In addition, MSE results showed that the MMSE-DFE for the fractional-order massive MIMO channel model with α = 0.1 and the MMSE-DFE for the integer-order massive MIMO channel model with n = 7 performed better than the MMSE-DFE for the FIR massive MIMO channel model.
The fractional-order massive MIMO channel model and integer-order massive MIMO channel model also showed an improved symbol error rate performance when compared to the FIR massive MIMO channel model. All these because of the reduced residual errors compared to the FIR model. This shows that the fractional-order massive MIMO channel model can be a more robust choice than the FIR massive MIMO channel model when it comes to channel estimation and equalization. It also shows that for the integer-order massive MIMO channel model to match the performance of the fractional-order massive MIMO channel model a higher system order should be selected. These results continue to show the viability of the fractional-order models in system identification, i.e., channel estimation and in channel equalization in wireless communications. They also show that fractional-order models can play a major role in improving performance in 5G networks.

APPENDIX A DEFINITIONS OF FRACTIONAL-ORDER CALCULUS A. THE GRÜNWALD-LETNIKOV DEFINITION
The α th order Grünwald-Letnikov (G-L) derivative of a function f (t) is written as [9] and [35]  where α j = (α+1) (j+1) (α−j+1) , the notation GL indicates ''Grünwald-Letnikov definition'', the operator D α defines fractional differentiation or integration depending on the sign of α. If α is negative then D −α will define fractional integration, and if α is positive then D α will define fractional differentiation, (·) is the Euler's Gamma function, [·] means to round off to the nearest integer, t 0 is the initial time instant and h is the finite sampling interval.

APPENDIX C THE POISSON MOMENT FUNCTIONAL FILTERING FOR CONTINUOUS-TIME FRACTIONAL-ORDER SYSTEM IDENTIFICATION
To handle the estimation of U 0|(i−1)α,N and Y 0|(i−1)α,N the PMF transform of order ris applied to (13), resulting in where P r t [D α g (t)]is the r th order PMF transform of signal D α g (t) at time instant t and is given by the following convolution product [41] P r t D α g (t) = D α g (t) * p r (t) = t 0 D α g (τ ) p r (t − τ ) dτ where p r (t) is the r th order Poisson pulse function expressed as p r (t) = β r+1 t r e −λt r! (C3) where r ∈ N, r ≥ i, whereby i is the user-defined index and λ ∈ R + and β ∈ R + are the Poisson filter constant and gain respectively. The r th order PMF of D α g (t)can be measured as the output at time instant t of a cascaded low-pass filter chain of (r + 1) identical stages, each with a transfer function expressed as [41] P r (s) = β s + λ (C4) According to [23], if we consider two α th differentiable functions g (t) and f (t), then In (C5) it can be seen that the α th order derivative of f (t) can be transmitted to the α th order derivative of g (t) and viceversa. Applying the property of (C5) to the (i − 1) α-order fractional derivatives of the input data, the following PMF based filtered input data is obtained P r t U 0|(i−1)α,N =      g (t) * u (t) D α g (t) * u (t) . . .
where f (t) in (C5) is replaced by u (t). In (C6) the α th order derivatives of g (t) can be easily obtained using the Riemann-Liouville, Grünwald-Letnikov or Caputo definitions. Applying the property of (C5) to the (i − 1) α-order fractional derivatives of the output data and the noise, the following PMF based filtered output data and noise is obtained and P r t N 0|(i−1)α,N =      g (t) * n (t) D α g (t) * n (t) . . .
respectively. According to [42] the PMF filter order and the PMF filter parameters are chosen as follows: The PMF filter order, r has to respect the following condition and for the PMF filter parameters, it is generally assumed that: This assumption helps to reduce the number of design parameters.