**3. MATLAB-based Algorithm development**

The systematic process involved in the development of the MATLAB-based SVD-ARMA algorithm for multiexponential transient signal analysis suitable for real-time application is discussed in this section. Apart from performance evaluation and signal generation, the algorithm consists of five major steps: obtaining convolution integral of the exponential signal using modified Gardner transformation; signal interpolation using spline technique; generation of deconvolved data; SVD-ARMA modeling of the deconvolved data; and power spectrum computation to finally estimate the transient signal parameters. A brief summary of the steps involved are as shown in Figure 2, and briefly highlighted as follows:

### **Step 1.** Signal generation

426 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

**Figure 1.** Overview of Techniques for Multiexponential signal analysis

of these techniques can be obtained in (Jibia, 2009; Salam and Sidek, 2000).

Parametric techniques such as Autoregressive (AR), Moving Average (MA) and ARMA models are introduced to the analysis of multiexponential signals to alleviate the drawbacks of the nonparametric techniques. The AR modeling technique requires less computation than the ARMA modeling technique as its model parameters are relatively easy to estimate. However, it is sensitive to noise due to the assumed all-pole model. On the other hand, MA parameters are difficult to estimate and the resultant spectral estimates have poor resolution. Although, the ARMA modeling technique is much better in estimating noisy signal than the AR modeling technique, it requires a lot of computation. A detailed review Generate the required signal, *S*() from MATLAB or from the fluorescence substances. For simulation data, MATLAB inbuilt function is used to generate the white Gaussian noise and the DC offset.

**Figure 2.** Flowchart of the MATLAB-based algorithm for Multicomponent transient signal analysis

#### **Step 2.** Signal preprocessing via Gardner transformation:

This involves the conversion of the measured or generated signal *S*( ) and *p*( ) to y(t) and h(t) respectively using modified Gardner transformation (Salami and Sidek, 2003). This yields a convolution integral as subsequently described.

In general, equation (1) is expressed as

$$S(\tau) = \sum\_{k=1}^{M} A\_k p(\mathcal{A} \kappa \tau) + n(\tau),\tag{2}$$

where the basis function, *p*() *=* exp(*-*). This equation can also be expressed as

$$S(\tau) = \bigcap\_{\alpha}^{\circ} \mathcal{g}(\mathcal{A}) p(\mathcal{A} \tau) d\mathcal{A} + \eta(\tau),\tag{3}$$

where 1 ( ) ( ). *M <sup>k</sup> <sup>k</sup> k g A* 

Both sides of equation (3) are multiplied by in the modified Gardner transformation instead of only in the original Gardner transformation. The value of the modifying parameter, is carefully chosen based on the criteria given by Salami (1995) to avoid poor estimation of the signal parameters because modifies the amplitude of the signal. Salami (1995) suggested to use 0 1 in the analysis of multiexponential signals. According to Salami (1985), the choice of can lead to improved signal parameters estimation as it produces noise reduction effect. The nonlinear transformation  *= e-r* and  *= et* are applied to equation (3), resulting in a convolution integral of the form

$$y(t) = \int\_{-\infty}^{\infty} x(\lambda)h(t - \lambda)d\lambda + v(t),\tag{4}$$

where the output function, *y*(*t*) *=* exp(*t*) *S*{exp(*t*)}, the input function, *x*(*t*) *=* exp{(*-*1)*t*}*g(e-t)*, the impulse response function of the system, *h*(*t*) *=* exp(*t*)*p*(*et* ) and the additive noise, *v*(*t*) *=*  exp(*t*)*n*(*et* ).

The discrete impulse response function, *h*[*n*] is obtained by sampling *<sup>α</sup>p*(*λk*) at 1/*Δt* Hz. Later, *H*(*k*) is obtained from the discrete Fourier transform (DFT) of *h*[*n*]. Next, equation (4) is converted into a discrete-time convolution. This is done by sampling *y*(*t*) at a rate of 1/*Δt* Hz to obtain

$$\text{lg}[n] = \sum\_{m=-\text{N min}}^{\text{N max}} \text{x}[m]h[n-m] + v[n] \tag{5}$$

where the total number of samples, *N* equals *nmax -nmin+*1, both *nmax* and *nmin* represent respectively the upper and lower data cut-off points. The criteria for the selection of these sampling conditions have been thoroughly discussed by Salami (1987) and Sen (1995) and these are not discussed further here. Equation (5) forms the basis of estimating the signal parameters since ideally *x*(*n*) can be recovered from the observed data by deconvolution. It is necessary to interpolate the discrete time convolution of *y*[*n*] since the log samples, *n =*  exp(*n∆t*) are not equally spaced.

#### **Step 3.** Cubic spline Interpolation

428 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

h(t) respectively using modified Gardner transformation (Salami and Sidek, 2003). This

( ) ( ) ( ),

*S gp d n* ( ) ( ) ( ) ( ),

 

*<sup>k</sup> <sup>k</sup>*

in the original Gardner transformation. The value of the modifying

in the analysis of multiexponential signals. According to

*t*) *S*{exp(*t*)}, the input function, *x*(*t*) *=* exp{(

*t*)*p*(*et*

can lead to improved signal parameters estimation as it

 *= e-r* and

(4)

(5)

is carefully chosen based on the criteria given by Salami (1995) to avoid poor

*y*( ) ( ) ( ) ( ), *t x ht d vt* 

Later, *H*(*k*) is obtained from the discrete Fourier transform (DFT) of *h*[*n*]. Next, equation (4) is converted into a discrete-time convolution. This is done by sampling *y*(*t*) at a rate of 1/*Δt*

[ ] [ ] [ ] [ ],

where the total number of samples, *N* equals *nmax -nmin+*1, both *nmax* and *nmin* represent respectively the upper and lower data cut-off points. The criteria for the selection of these sampling conditions have been thoroughly discussed by Salami (1987) and Sen (1995) and

*y n xmhn m vn*

max

*n*

*m n*

min

The discrete impulse response function, *h*[*n*] is obtained by sampling

   

). This equation can also be expressed as

1

*k S Ap n*

0

 *M*

(2)

(3)

in the modified Gardner transformation

modifies the amplitude of the signal. Salami

 *= et*

are applied to

) and the additive noise, *v*(*t*) *=* 

*<sup>α</sup>p*(*λk* *-*1)*t*}*g(e-t)*,

) at 1/*Δt* Hz.

 and *p*( ) 

to y(t) and

This involves the conversion of the measured or generated signal *S*( )

**Step 2.** Signal preprocessing via Gardner transformation:

yields a convolution integral as subsequently described.

) *=* exp(*-*

In general, equation (1) is expressed as

where the basis function, *p*(

1 ( ) ( ). *M*

*k g A* 

instead of only

parameter,

exp(*t*)*n*(*et* ).

Hz to obtain

*<sup>k</sup> <sup>k</sup>*

Both sides of equation (3) are multiplied by

estimation of the signal parameters because

equation (3), resulting in a convolution integral of the form

the impulse response function of the system, *h*(*t*) *=* exp(

produces noise reduction effect. The nonlinear transformation

 

(1995) suggested to use 0 1

where the output function, *y*(*t*) *=* exp(

Salami (1985), the choice of

where

The discrete-time signal y[n] obtained in (5) consists of non-equally spaced samples which would be difficult to digitally process. A cubic interpolation is therefore applied to y[n] using MATLAB function 'spline' to obtain equally spaced samples of y[n].

#### **Step 4.** Generation of deconvolved data

In this stage, deconvolved data is generated from y[n] using optimally compensated inverse filtering due to its ability to handle noisy signal when compared to conventional inverse filtering (Salami,1995).

Conventionally, this is done by taking the DFT of equation (5) to produce:

$$Y(k) = X(k)H(k) + V(k)\tag{6}$$

The deconvolved data, *X k*( ) can be obtained by computing *Y*(*k*)/*H*(*k*), that is

$$
\hat{X}(k) = \frac{Y(k)}{H(k)} = X(k) + \frac{V(k)}{H(k)}\quad\text{for}\ 0 \le k \le N - 1\tag{7}
$$

where *Y*(*k*), *X*(*k*), *H*(*k*) and *V*(*k*) represent the discrete Fourier transform of *y*(*n*), *x*(*n*), *h*(*n*) and *v*(*n*) respectively. This inverse filtering operation is called the conventional inverse filtering. It yields deconvolved data with decreasing SNR for increasing values of *k*. Therefore, the accuracy of *X k*( ) deteriorates when the noise variance level is high.

To overcome this problem, an optimally compensated inverse filtering is introduced. In this approach, *H*(*k*) is modified by adding an optimally selected value, into it. This procedure is done according to Riad and Stafford (1980) by initially designing a transfer function, *HT* (*k*) that yields a better *X* (*k*) in equation (7), where *HT* (*k*) is given by:

$$H\_T(k) = \frac{H^\*(k)}{\left[\lceil H(k) \rceil^2 + \mu\right]},\tag{8}$$

where denotes the complex conjugate. Substituting equation (8) into equation (7) yields

$$\overset{\circ}{X}(k) = \frac{Y(k)H^\*(k)}{\left[\|H(k)\|^2 + \mu\right]},\tag{9}$$

which is referred to as the optimally compensated inverse filtering. It is noted that a small value of *μ* has a very little effect in the range of frequency when |*H*(*k*)|2 is significantly larger than . However, if|*H*(*k*)|2 is very small, the effect of *μ* on the deconvolved data is quite substantial, that is tends to make ( ) *X k* less noisy. Therefore, *μ* puts limit on the noise amplification because the denominator becomes lower bounded according to Dabóczi and Kollár (1996). The parameter *μ* is carefully selected according to the SNR of the data to obtain good results. The choice of the optimum value of *μ* according to Salami and Sidek (2001) is best determined by experimental testing.

Equation (9) shows one-parameter compensation procedure, however, multi-parameter compensation is considered in this study. Thus, a regularization operator L(k) is introduced into (8), that is

$$F(k) = \frac{H^\*(k)}{\left[\left|H(k)\right|^2 + \alpha \left|L(k)\right|^2\right]}\tag{10}$$

where is the controlling parameter and the regularization operator, *L*(*k*) is the discrete Fourier transform of the second order backward difference sequence. |*L*(*k*)|2 is given as

$$\|L(k)\|^2 = 16\sin^4\left(\frac{\pi k}{N}\right). \tag{11}$$

where *N* is the number of samples. Using both the second and fourth order backward difference operations in equation (10) yield

$$F(k) = \frac{H^\*(k)}{\left[\left|H(k)\right|^2 + \mu + \alpha \left|L(k)\right|^2 + \beta \left|L(k)\right|^4\right]}\tag{12}$$

where <sup>4</sup> *L k*( ) denotes the fourth order backward difference operator and , and are the varied compensation parameters to improve the SNR of the deconvolved data.

Unwanted high frequency noise can still be introduced by this optimally compensated inverse filtering which can make some portion of *X k*( ) unusable. Therefore, a good portion of ˆ *X k*( ) denoted as *f* (*k*), is given as

$$f(k) = \sum\_{i=1}^{M} B\_i \exp(j2\pi k \Delta f \ln \lambda\_i) + V(k),\tag{13}$$

where 1 *k* 2*Nd+*1, �� � ��/ � �, *Nd* ( / 2) 1 *N* , *Nd* is the number of useful deconvolved data points, *N* is the number of data samples and *V*(*k*) is the noise samples of the deconvolved data. Equation (13) is interpreted as a sum of complex exponential signals. The number of deconvolved data points, *Nd* is carefully selected to produce good results from *f* (*k*).

#### **Step 5.** Signal parameters estimation using SVD-ARMA modeling

SVD-ARMA algorithm is applied to f(k) to estimate the signal parameters *M* and *k.* as it provides consistent and accurate estimates of AR parameters with minimal numerical problem which is necessary for real-time application. In addition, it is a powerful computational procedure for matrix analysis especially for solving over determined system of equations. The detailed mathematical analysis of this algorithm is reported in (Salami, 1985).

Generally, the ARMA model assumes that *f* (*k*) satisfies the linear difference equation (Salami 1985)

$$f(k) = -\sum\_{i=1}^{p} a(i)f(k-i) + \sum\_{i=1}^{q} b(i)V(k-i) \tag{14}$$

where V(k) is the input driving sequence, *f(k)* is the output sequence, *a*[i] and *b*[i] are the model coefficients with AR and MA model order of *p* and *q* respectively. Usually, the white Gaussian noise becomes the input driving sequence in the analysis of exponentially decaying transient signals.

One of the most effective procedures for estimating these model parameters is by solving a modified Yule-Walker equation (Kay and Marple 1981). This procedure is subsequently discussed.

Equation (14) is multiplied by *f \**(*k-m*) and taking the expectation yields:

430 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

tends to make ( ) *X k*

(2001) is best determined by experimental testing.

difference operations in equation (10) yield

*X k*( ) denoted as *f* (*k*), is given as

where 1 *k* 2*Nd+*1, �� � ��/ 

inverse filtering which can make some portion of *X k*( )

�

. However, if|*H*(*k*)|2 is very small, the effect of *μ* on the deconvolved data is

2 2

() ()

is the controlling parameter and the regularization operator, *L*(*k*) is the discrete

*N* 

2 24

 

(13)

�, *Nd* ( / 2) 1 *N* , *Nd* is the number of useful deconvolved data

 

() () ()

*Hk Lk Lk* 

Unwanted high frequency noise can still be introduced by this optimally compensated

( ) exp( 2 ln ) ( ),

points, *N* is the number of data samples and *V*(*k*) is the noise samples of the deconvolved data. Equation (13) is interpreted as a sum of complex exponential signals. The number of

provides consistent and accurate estimates of AR parameters with minimal numerical

*fk B j kf Vk* 

*i i*

 

*Hk Lk* 

 

less noisy. Therefore, *μ* puts limit on the

(10)

(11)

, and are the

unusable. Therefore, a good portion

(12)

*k.* as it

noise amplification because the denominator becomes lower bounded according to Dabóczi and Kollár (1996). The parameter *μ* is carefully selected according to the SNR of the data to obtain good results. The choice of the optimum value of *μ* according to Salami and Sidek

Equation (9) shows one-parameter compensation procedure, however, multi-parameter compensation is considered in this study. Thus, a regularization operator L(k) is introduced

( ) ( )

Fourier transform of the second order backward difference sequence. |*L*(*k*)|2 is given as

2 4 | ( )| 16sin , *<sup>k</sup> L k*

where *N* is the number of samples. Using both the second and fourth order backward

( ) ( )

*H k F k*

varied compensation parameters to improve the SNR of the deconvolved data.

where <sup>4</sup> *L k*( ) denotes the fourth order backward difference operator and

1

deconvolved data points, *Nd* is carefully selected to produce good results from *f* (*k*).

SVD-ARMA algorithm is applied to f(k) to estimate the signal parameters *M* and

*i*

**Step 5.** Signal parameters estimation using SVD-ARMA modeling

*M*

*H k F k*

larger than

into (8), that is

where 

of ˆ

quite substantial, that is

$$R\_{ff}(k) = -\sum\_{n=1}^{p} a[n]R\_{ff}(k-n) + \sum\_{n=0}^{q} b[n]h(k-m),\tag{15}$$

where *Rff* (*k*) is the autocorrelation function of *f k*( ) and *h*(*k*) is the impulse response function of the ARMA model. Next, considering the AR portion of equation (15) leads to the modified Yule-Walker equation

$$\mathcal{R}\_{\mathcal{f}\mathcal{f}}(k) + \sum\_{n=1}^{p} a[n]\mathcal{R}\_{\mathcal{f}}(k-n) = 0; k \ge q+1. \tag{16}$$

Equation (16) may not hold exactly in practice because both *p* and *q* are unknown prior to analysis and *Rff* (*k*) has to be estimated from noisy data. This problem is solved by using an SVD algorithm. This algorithm is used by first expressing equation (16) in matrix form as **Ra**  = **e** with **R** having elements ( , ) ( 1 ), *ff e ril R q i l* where 1 *i r* ; 1 *l pe* + 1. Note that both *pe* and *qe* are the guess values of the AR and MA model order respectively, **a** is a *pe* 1 and **e** is a *r*  1 error vector with *r* > *pe*. The SVD algorithm is applied to **R** to produce

$$\mathbf{R} = \mathbf{U}\mathbf{Z}\mathbf{V}^{\mathrm{T}=} \sum\_{n=1}^{p\_c} \sigma\_n \boldsymbol{\mu}\_n \boldsymbol{\upsilon}\_n^H \tag{17}$$

where the *<sup>r</sup>* (*pe*+1) unitary matrix **U** = [*u*<sup>1</sup> *<sup>u</sup>*2 … <sup>1</sup> *<sup>e</sup> <sup>p</sup> <sup>u</sup>* ], (*pe+*1) (*pe+*1) unitary matrix **V =** [*v*<sup>1</sup> *<sup>v</sup>*<sup>2</sup> … <sup>1</sup> *<sup>e</sup> <sup>p</sup> <sup>v</sup>* ] and is a diagonal matrix with diagonal elements ( 1, 2 … <sup>1</sup> *<sup>e</sup> <sup>p</sup>* ). These diagonal elements are called singular values and are arranged so that 1 > 2 > … > <sup>1</sup> *<sup>e</sup> <sup>p</sup>* > 0. Only the

first *M* singular values will be nonzero so that *M+*1 = *M+*2 = … = <sup>1</sup> *<sup>e</sup> <sup>p</sup>* = 0. However, *M+*<sup>1</sup> *M+*<sup>2</sup> … <sup>1</sup> *<sup>e</sup> <sup>p</sup>* 0 due to noise contamination. The problem is solved by constructing a lower rank matrix **RL** from **R** using the first M singular values, that is

$$\mathbf{RL} = \mathbf{U}\mathbf{M}\mathbf{\Delta}\mathbf{V}\mathbf{u}^{H} = \sum\_{n=1}^{M} \sigma\_{n}\boldsymbol{u}\_{n}\boldsymbol{v}\_{n}^{\ H} \tag{18}$$

where **UM**, **M** and **VM** are the truncated version of **U**, and **V** respectively. The AR coefficients are then estimated from the relation **a** = -**RL**#**r**, where **r** corresponds to the first column of **RL** and **RL**# is given as

$$\mathbf{RL}^{\mathfrak{s}} = \sum\_{n=1}^{M} \sigma\_n^{-1} u\_n^H v\_n. \tag{19}$$

The estimated AR coefficients are then used to generate the residual error sequences:

$$\beta(k) = \sum\_{l=0}^{\frac{p\_c}{c}} \sum\_{m=0}^{p\_c} a[l] a^\* \left[m\right] \mathbb{R}\_{\frac{p}{\ell'}}(k+m-l) \tag{20}$$

from which the actual MA parameters are obtained directly from equation (20). However, MA spectra can be obtained from the DFT of the error sequences, (*k*) . An exponential window is applied to (*k*) to ensure that the MA spectra derived from the error sequences are positive definite. Next, the ARMA spectrum is computed from

$$S\_f(\mathbf{z}) = \frac{\sum\_{k=-p\_\varepsilon}^{p\_\varepsilon} \beta(k)\mathbf{z}^{-k}}{\left|A(\mathbf{z})\right|^2} \tag{21}$$

and the desired power distribution of *x*(*t*), denoted as *Px*(*t*) is obtained by evaluating *Sf* (*z*) on the unit circle exp 2 *<sup>t</sup> z j N t* , that is:

$$P\_x(t) = S\_f(z)\bigg|\_{z=\exp\left\{\frac{j2\pi t}{N\Lambda t}\right\}} = \sum\_{k=1}^{M} B\_k^2 \delta(t - \ln \mathcal{A}\_k). \tag{22}$$

Eventually, *M* and ln *<sup>k</sup>* are obtained from *Px*(*t*).

#### **Step 6.** Graphical presentation of output

Power distribution graph has been used to display the results of multiexponential signal analysis. This is computed from the power spectrum of the resulting output signal from SVD-ARMA modeling method as shown in equation (21).

**Step 7.** Performance evaluation

The efficiency of the algorithm with SVD-ARMA modeling technique in estimating *k* correctly is determined by the Cramer-Rao lower bound (CRLB) expressed as:

Matlab-Based Algorithm for Real Time Analysis of Multiexponential Transient Signals 433

$$\text{CRLB ( $\mathcal{A}\_k$ )}=\frac{6(1+0.7(\mathcal{A}\_k N)^{3/2})^2}{N^3 \text{SNR}},\tag{23}$$

where *N* is the number of data points, the CRB(*k*) is the CRLB for estimating *<sup>k</sup>* and SNR equals to ܣ ଶ divided by the variance of the white Gaussian noise.

The Cramer-Rao lower bound will determine whether the estimator is efficient by comparing the variance of the estimator, *var*( *k* ) with the Cramer-Rao lower bound. Variance that approaches the CRLB is said to be optimal according to Kay and Marple (1981) and Sha'ameri (2000). Consequently, the closer the variance of the estimator is to the CRLB, the better is the estimator.

A MATLAB-mfile code has been written to implement the steps 1 to 6 described above, details of this have been thoroughly discussed in (Jibia, 2009).

## **4. Integrated MATLAB-labview for real-time implementation**

This section discusses the development of proposed integrated Labview-MATLAB software interface for real-time (RT) implementation of the algorithm for multicomponent signal analysis as described in section 3. Real-time signal analysis is required for most practical applications of multicomponent signal analysis. In this study, reference is made to the application to fluorescence signal analysis. A typical multicomponent signal analyzer comprises optical sensor which is part of the spectrofluorometer system, signal conditioning/data acquisition system, embedded processor that runs the algorithm in realtime, and display/storage devices as shown in Figure 3. National Instrument (NI) real-time hardware and software are considered for the development of this system due to its ease of implementation.

**Figure 3.** Block diagram of the multicomponent signal analyzer prototype.

#### **4.1. Labview real-time module and target**

432 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

lower rank matrix **RL** from **R** using the first M singular values, that is

 *M+*1 = 

0 due to noise contamination. The problem is solved by constructing a

1

*H n nn*

.

*ff*

2

( )

*A z*

( )

*k z*

*k*

2

 (22)

*M*

*M*

*n u v*

where **UM**, **M** and **VM** are the truncated version of **U**, and **V** respectively. The AR coefficients are then estimated from the relation **a** = -**RL**#**r**, where **r** corresponds to the first

1

( ) [ ] \*[ ] ( ) *e e <sup>p</sup> <sup>p</sup>*

*k ala mR k m l*

from which the actual MA parameters are obtained directly from equation (20). However,

*e*

*p*

*k p*

and the desired power distribution of *x*(*t*), denoted as *Px*(*t*) is obtained by evaluating *Sf* (*z*) on

2 exp 1 () ( ) ( ln ).

*x f j t k k <sup>z</sup> <sup>k</sup> N t Pt S z B t* 

Power distribution graph has been used to display the results of multiexponential signal analysis. This is computed from the power spectrum of the resulting output signal from

The efficiency of the algorithm with SVD-ARMA modeling technique in estimating

correctly is determined by the Cramer-Rao lower bound (CRLB) expressed as:

 

*e*

*M*

*n u v*

The estimated AR coefficients are then used to generate the residual error sequences:

0 0

( )

*f*

, that is:

*<sup>k</sup>* are obtained from *Px*(*t*).

SVD-ARMA modeling method as shown in equation (21).

*S z*

*l m*

MA spectra can be obtained from the DFT of the error sequences,

are positive definite. Next, the ARMA spectrum is computed from

 *M+*2 = … = <sup>1</sup> *<sup>e</sup> <sup>p</sup>* 

,

*H nnn*

= 0. However,

(18)

(19)

(21)

(*k*) . An exponential

*k*

(20)

(*k*) to ensure that the MA spectra derived from the error sequences

 *M+*<sup>1</sup> 

first *M* singular values will be nonzero so that

 **RL** = **UMMVM***<sup>H</sup>*

 **RL**# <sup>1</sup>

**Step 6.** Graphical presentation of output

**Step 7.** Performance evaluation

column of **RL** and **RL**# is given as

window is applied to

the unit circle exp 2 *<sup>t</sup> z j N t*

Eventually, *M* and ln

*M+*<sup>2</sup> … <sup>1</sup> *<sup>e</sup> <sup>p</sup>* 

> The Labview real-time module (RT-software) together with NI sbRIO-9642 (RT-hardware) has been adopted in this study. The Labview Real-Time software module allows for the creation of reliable, real-time applications which are easily downloaded onto the target hardware from Labview GUI programming tool.

The NI sbRIO-9642 is identical in architecture to CompactRIO system, only in a single circuit board. Single-Board RIO hardware features a real-time processor and programmable FPGA just as with CompactRIO, and has several inputs and outputs (I/O) modules as shown in Figure 4.

System development involves graphical programming with Labview on the host Window computer, which is then downloaded and run on real-time hardware target. Since the algorithm has been developed with MATLAB scripts, an integrated approach was adopted in the programme development as subsequently described.

**Figure 4.** NI sbRIO-9642 for real-time hardware target
