## 5 Numerical Aspects  
The plots presented in the above figures point to (24) as the easiest way of computing the solution of our problem. As stated, the coefficients decrease to zero. From a numerical  
Figure 2: Evolution of the coefficients in the series (18) as a function of the order: \(B_{n},n=0,1,\cdots\) in (19).  
point of view, this means that they become comparable to the relative spacing between any two adjacent numbers in the machine's floating point system, \"eps\" in Octave/Matlab. We set \(n_{th}=\lfloor\frac{-\ln(eps)}{ln(q)a}\rfloor\). Therefore, the series in (18) becomes a polynomial.  
We have two different numerically equivalent approaches to compute such summations. The most obvious is the use of (18) together with a suitable algorithm for the computation of the needed MLF [47]. We performed this, but the algorithm has a very severe drawback: the computational burden. For example, taking \(q=1.2\), with a sampling interval equal to \(0.01\), it took around \(17\) minutes to compute \(x(t)\). Therefore, we looked for a different approach. This consisted of approximating a continuous-time function to a discrete-time one [8]. To accomplish this, we worked in frequency domain, using (24) and the bilinear (Tustin) transformation, followed by a (discrete) Fourier-transform inversion. We proceeded according to the following steps.  
1. Use the \(s\to z\) bilinear transformation [8] \[s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}},\] to pass from a continuous-time to a sampled (discrete) function. \(T\) is the sampling interval. We obtain the function \[X_{d}(z)=\frac{T}{2}\frac{1+z^{-1}}{1-z^{-1}}\sum_{n=0}^{n_{th}}C_{n}\frac{ \left(\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}\right)^{\alpha}}{\left(\frac{2}{T }\frac{1-z^{-1}}{1+z^{-1}}\right)^{\alpha}+q^{-n\alpha}}.\] (30)
2. Set \(z=e^{i\omega},\ -\pi<\omega\leq\pi\), to get a discrete-time Fourier transform \(X_{d}(e^{i\omega})\) of \(x(nT)\)[39]. To avoid the effect of the branchcut points at \(z=\pm 1\), we found it better to move them slightly into the unit circle. This is accomplished by setting \[s=\frac{2}{T}\frac{1-0.99e^{-i\omega}}{1+0.99e^{-i\omega}}.\] We can use a simple procedure that consists of the following steps: 1. Set \(N(e^{i\omega})=\mathcal{F}([1,-0.99])\), where \(\mathcal{F}\) represents the discrete Fourier transform that is implemented by the fast Fourier transform algorithm; 2. Similarly, set \(D(e^{i\omega})=\mathcal{F}([1,+0.99])\); 3. Compute \(s(\omega)=\frac{2}{T}\frac{N(e^{i\omega})}{D(e^{i\omega})}\); 4. Obtain the approximation to \(X(s)\) \[X_{d}(s(\omega))=\frac{T}{2}\frac{1}{s(\omega)}\sum_{n=0}^{n_{th}}C_{n}\frac{ s^{\alpha}(\omega)}{s^{\alpha}(\omega)+q^{-n\alpha}}.\] (31)
3. Sample \(X_{d}(e^{i\omega})\) in a uniform grid, making \(\omega_{k}=\frac{2\pi T}{Nk},\ k=0,1,\cdots,N-1\), where \(N\) is in high enough agreement with \(T\) and the interval where we want to compute the function. In our simulations, we used \(T=0.001\) and \(N=2\)[16].
4. We could use the fast Fourier transform to compute the inverse of \(X_{d}(e^{i\omega_{k}})\), \(x(nT)\). However, this is not very good from numerical aspects, due to the singularity at the origin. We found it was better to separate the computation in three steps: 1. Computation of the summation in (32), making a slight transformation \[\bar{X}_{d}(e^{i\omega})=\sum_{n=0}^{n_{th}}C_{n}\frac{1}{1+q^{-n_{th}s-\alpha }(\omega)};\] (32)