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

Multiexponential transient signals are particularly important due to their occurrences in many natural phenomena and human applications. For instance, it is important in the study of nuclear magnetic resonance (NMR) in medical diagnosis (Cohn-Sfetcu et al., 1975)), relaxation kinetics of cooperative conformational changes in biopolymers (Provencher, 1976), solving system identification problems in control and communication engineering (Prost and Guotte, 1982), fluorescence decay of proteins (Karrakchou et al., 1992), fluorescence decay analysis (Lakowicz, 1999). Several research work have been reported on the analysis of multicomponent transient signals following the pioneer work of Prony in 1795 (Prony, 1975) and Gardner et al. in 1959 (Gardner, 1979). Detailed review of several techniques for multicomponent transient signals’ analysis was recently reported in (Jibia, 2010).


Introduction
Multiexponential transient signals are particularly important due to their occurrences in many natural phenomena and human applications. For instance, it is important in the study of nuclear magnetic resonance (NMR) in medical diagnosis (Cohn-Sfetcu et al., 1975)), relaxation kinetics of cooperative conformational changes in biopolymers , solving system identification problems in control and communication engineering (Prost and Guotte, 1982), fluorescence decay of proteins (Karrakchou et al., 1992), fluorescence decay analysis (Lakowicz, 1999). Several research work have been reported on the analysis of multicomponent transient signals following the pioneer work of Prony in 1795 (Prony, 1975) and Gardner et al. in 1959(Gardner, 1979. Detailed review of several techniques for multicomponent transient signals' analysis was recently reported in (Jibia, 2010).
Generally, a multiexponential transient signal is represented by a linear combination of exponentials of the form where M is the number of components, Ak and k respectively represent the amplitude and real-valued decay rate constants of the kth component and n() is the additive white Gaussian noise with variance n 2 . The exponentials in equation (1) are assumed to be separable and unrelated. That is, none of the components is produced from the decay of another component. Therefore, in determination of the signal parameters, M, Ak and k from equation (1), it is not sufficient that equation (1) approximates data accurately; it is also important that these parameters are accurately estimated.
There are many problems associated with the analysis of transient signals of the form given in equation (1) due to the nonorthogonal nature of the exponential function. These problems include incorrect detection of the peaks, poor resolution of the estimated decay and inaccurate results for contaminated or closely-related decay rate data as reported in (Salami et.al., 1985). These problems become increasingly difficult when the level of noise is high. Although Gardner transform eliminated the nonorthogonality problem, it introduced error ripples due to short data record and nonstationarity of the preprocessed data. Apart from these problems, analysis of multiexponential signal is computationally intensive and requires efficient tools for its development and implementation in real-time.
To overcome these problems, modification of Gardner transform has been proposed recently with Multiple Signal Classification (MUSIC) superposition modeling technique (Jibia et al., 2008); with minimum norm modeling technique (Jibia and Salami, 2007), with homomorphic deconvolution, with eigenvalues decomposition techniques, and the Singular Value Decomposition (SVD) based-Autoregressive Moving Average (ARMA) modeling techniques (Salami and Sidek, 2000;Jibia, 2009). As reported in (Jibia, 2009), performance comparison of these four modeling techniques has been investigated. Though, all the four techniques were able to provide satisfactory performances at medium and high signal-to-noise ratio (SNR), the SVD-ARMA was reported to have the highest resolution, especially at low SNR.
Hence, the development of SVD-ARMA based algorithm for multiexponential signal analysis using MATLAB software package is examined in this chapter.
MATLAB provides computational efficient platform for the analysis and simulation of complex models and algorithms. In addition, with the aid of inbuilt embedded MATLAB Simulink block, it offers a tool for the integration of developed algorithm/model in an embedded application with little programming efforts as compared to the use of other programming languages (Mathworks, 2008). This functionality is explored in integrating the developed MATLAB-based algorithm into National Instrument (NI) Labview embedded programming tool. Hence, an integrated MATLAB-Labview software interface is proposed for real-time deployment of the algorithm. To this end, the analytical strength of MATLAB together with simplicity and user-friendly benefits of the National Instrument (NI), Labview design platforms are explored in developing an efficient, user-friendly algorithm for the real-time analysis of multiexponential transient signal.
The rest of the chapter is organized as follows. Section 2 provides brief review of techniques for multiexponential signal analysis. The MATLAB algorithm development for the signal analysis is presented in section 3. The development of an integrated MATLAB-Labview realtime software interface is then examined in section 4. Section 5 presents sample real-time data collection together with results and analysis. The chapter is concluded in section 6 with recommendation for future study.

Techniques of multiexponential transient signal analysis
Several techniques have been reported for the analysis of transient multiexponential signal. They are classified as: (i) time domain or frequency domain, and (ii) parametric or nonparametric techniques. The main objective of these techniques in analyzing the multiexponentially decaying signals is to estimate the signal parameters as accurately as possible and to get better display of the signal spectra. Time-domain techniques constitute the oldest methods of multiexponential signal analysis prior to the advent of Gardner transformation technique (Gardner et al, 1959;Salami, 1985). Gardner transformation is one of the most important methods of the transient signal analysis based on spectral analysis. Generally, spectral analysis involves transformation of a time-domain signal to a frequency domain so that certain features of the signals that characterized them are easily discerned such as its decay constant, k and amplitude, Ak. In other words, it is the process of obtaining the frequency content (spectrum) of a signal (Proakis and Manolakis 1996). The spectral analysis approach is further categorized into nonparametric and parametric techniques. Nonparametric technique is a frequency-domain technique that obtains the signal spectra directly from the deconvolved data, while the parametric technique obtains the signal spectra indirectly by determining a finite set of parameters that defines a closed form mathematical model for the deconvolved data. Therefore the techniques of multiexponential transient signal analysis are sub-divided into time-domain, nonparametric frequency domain and parametric frequency domain techniques as shown in Figure 1 with their associated methods.
Among the earliest time-domain technique is the peeling technique. However, this technique produces poor results when ( ) S  contains more than two components. Other time-domain techniques such as Prony's method and its variants produce better performance than the peeling technique, however they are very sensitive to noise (Smyth, 2002;Salami, 1995). The nonlinear least squares technique (Smyth, 2002) is computationally inefficient and the solution sometimes fails to converge, which means the estimate of the signal parameters cannot be accurately obtained.
The nonparametric techniques of spectral analysis are introduced to overcome some limitations of the time-domain techniques. The Gardner transformation technique produces error ripples which obscure the real peaks of the spectrum due to the cutoff points. This technique is good in analyzing signal with high SNR. The fast Fourier transform (FFT) technique, which is an improvement over the original Gardner transformation produces improved resolution. However, the problem of error ripples still exists. The extension of this technique involving the use of digital signal processing and Gaussian filtering is sensitive to noise and its data range has to be limited to get accurate estimates of the signal parameters. Whilst the modified FFT technique is better than the previous methods, it often fails to estimate k correctly especially when the peaks are closely related. The differential technique (Swingler, 1977) provides some improvements over the existing techniques such as better resolution display but it is not suitable for analyzing noisy signal. Furthermore, modifying the FFT technique by incorporating integration procedure (Balcou,1981) does not produce better results as compared to the previous digital technique and the modified FFT technique. Moreover, this technique is still affected by error ripples. 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 of these techniques can be obtained in (Jibia, 2009;Salam and Sidek, 2000).

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 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. 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 where the basis function, p() = exp(-). This equation can also be expressed as 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  = e t are applied to equation (3), resulting in a convolution integral of the form 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(e t ) and the additive noise, v(t) = exp(t)n(e 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 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 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: 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: where  denotes the complex conjugate. Substituting equation (8) into equation (7) yields 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 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 where N is the number of samples. Using both the second and fourth order backward difference operations in equation (10) yield where 4 ( ) 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 where 1  k  2Nd+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) 1 1 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: 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) 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) Eventually, M and ln k 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: where N is the number of data points, the CRB(k) is the CRLB for estimating k 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).

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.

Labview real-time module and target
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.

MATLAB-labview software integration
The developed algorithm with MATLAB scripts was integrated inside Labview for real-time embedded set-up as shown in Figure 5. Labview front-panel and block diagrams were developed with inbuilt Labview math Scripts to run the MATLAB algorithm for the analysis of multicomponent signals. The use of Labview allows for ease of programming, and realtime deployment using the Labview real-time module described in section 4.1. Also, it provides a user friendly software interface for real-time processing of the fluorescence signals.
As shown in Figure 5, the sampled data produced by the spectrofluorometer system are preprocessed by the NI-DAQ Cards. These signals are then read by the embedded real-time software, analyze the signals and display the results in a user friendly manner. The user is prompted to enter the number of samples to be analyzed from the front-panel using the developed SVD-ARMA algorithm. this requires re-structuring of the codes to be compatible for embedded Simulink implementation, and hence deployment inside Labview MathScript. Figure 7 shows the MATLAB Simulink blocks configuration with embedded MATLAB function together with cross-section of the algorithm. The DATA with time vector is prepared inside the workspace and linked to the model input. Figure 10 (a-f) shows the simulation results with experimental data described in section 5. iii. The integrated software interface is evaluated with the real-time fluorescence data collected from a spectrofluorometer.

Sample collection, results and discussion
Due to unavailability of spectrofluorometer system that can be directly linked with the setup, sampled data collected from fluorescence decay experiment conducted using Spectramax Germini XS system were used to test the performance of the integrated system. The schematic diagram of the spectrofluorometer operation is shown in Figure 9, and itemized as follows: Step 1. The excitation light source is the xenon flash lamp.
Step 2. The light passes through the excitation cutoff wheel. This wheel reduces the amount of stray light into the movable grating.
Step 3. The movable grating selects the desired excitation wavelength. Then, this excitation light enters a 1.0 mm diameter fiber.
Step 4. This 1.0 mm diameter fiber focuses the excitation light before entering the sample in the micro-plate well. This focusing prevents part of the light from striking adjacent wells.
Step 5. The light enters the wheel and if fluorescent molecules are present, the two mirrors focus the light from the well into a 4.0 mm optical bundles.
Step 6. The movable, focusing grating allows light of chosen emission wavelength to enter the emission cutoff wheel.
Step 7. This emission cutoff wheel will further filter the light before the light enters the photomultiplier tube.
Step 8. The photomultiplier tube detects the emitted light and passes a signal to the instrument's electronics which then send the signal to the data acquisition system inside the spectrofluorometer.
Three intrinsic fluorophores (Acridine Orange; Fluorescein Sodium and Quinine) were used in the experiment. The details of the substances are given in Table 1.  Figure 10 (a-f) respectively. Figure 11 to Figure 13 show the sample of results obtained from the integrated MATLAB-Labview real-time software which has been developed. Both results of simulation and realtime software interface yield accurate estimates of the fluorescence data as shown in Figure 10- Figure 13, and presented in Table 2. The singular values for each of the samples combination are given in Table 3. The results indicate accurate determination of the constituent samples

Conclusion/Future study
The development of MATLAB-based algorithm for real-time analysis of multicomponent transient signal analysis based on SVD-ARMA modeling technique has been presented in this chapter. To enhance real-time interface and rapid prototyping on target hardware, complementary benefits of MATLAB and Labview were explored to develop real-time software interface downloadable into single board computer by NI Labview. In the absence of the spectrofluorometer system, the developed user friendly software for real-time deployment was validated with the collected real-time data. The obtained results indicate the effectiveness of the proposed integrated software for practical application of the proposed algorithm.
Future direction of this research will be directed towards development of customized spectrofluorometer sub-system that can be directly integrated to the overall system. This will eventually facilitate direct application of the developed algorithm in practical applications involving transient signal analysis. Also, other algorithms based on homomorphic and eigenvalues decomposition developed by the authors in similar study are to be made available as option on the user interface.  Table 3. Singular values for SVD-ARMA using experimental data