 Software
 Open Access
 Published:
Coding Prony’s method in MATLAB and applying it to biomedical signal filtering
BMC Bioinformatics volume 19, Article number: 451 (2018)
Abstract
Background
The response of many biomedical systems can be modelled using a linear combination of damped exponential functions. The approximation parameters, based on equally spaced samples, can be obtained using Prony’s method and its variants (e.g. the matrix pencil method). This paper provides a tutorial on the main polynomial Prony and matrix pencil methods and their implementation in MATLAB and analyses how they perform with synthetic and multifocal visualevoked potential (mfVEP) signals.
This paper briefly describes the theoretical basis of four polynomial Prony approximation methods: classic, least squares (LS), total least squares (TLS) and matrix pencil method (MPM). In each of these cases, implementation uses general MATLAB functions. The features of the various options are tested by approximating a set of synthetic mathematical functions and evaluating filtering performance in the Prony domain when applied to mfVEP signals to improve diagnosis of patients with multiple sclerosis (MS).
Results
The code implemented does not achieve 100%correct signal approximation and, of the methods tested, LS and MPM perform best. When filtering mfVEP records in the Prony domain, the value of the area under the receiveroperatingcharacteristic (ROC) curve is 0.7055 compared with 0.6538 obtained with the usual filtering method used for this type of signal (discrete Fourier transform lowpass filter with a cutoff frequency of 35 Hz).
Conclusions
This paper reviews Prony’s method in relation to signal filtering and approximation, provides the MATLAB code needed to implement the classic, LS, TLS and MPM methods, and tests their performance in biomedical signal filtering and function approximation. It emphasizes the importance of improving the computational methods used to implement the various methods described above.
Background
Prony’s method
In 1795, Gaspard de Prony [1] proposed a method to explain the expansion of gases as a linear sum of damped complex exponentials of signals that are uniformly sampled. Prony’s method approximates a sequence of N = 2p equally spaced samples to a linear combination of p complex exponential functions with differing amplitudes, damping factors, frequencies and phase angles. The main contribution of this classic method is that it converts a nonlinear approximation of exponential sums by solving a set of linear equations and a rootfinding problem.
The conventional or polynomial Prony method consists of setting out an autoregressive model of order p that assumes that the value of sampled data x[n] depends linearly on the preceding p values in x. Solving this linear system of equations obtains the coefficients of the characteristic or Prony polynomial φ(z). The roots of this polynomial yield two of the parameters of the solution (damping factors and frequency) and provide a second system of equations to calculate the amplitude and phase of the p functions.
Prony’s original method exactly matched the curve of p exponential terms to a dataset of N = 2p elements. When N > 2p, the linear systems of equations are overdetermined and can be approximated by the least squares (LS) method [2]. The conventional leastsquares method considers that in the linear system (A.x ≈ b), only b (observation vector) is contaminated by noise, while A (coefficient matrix) is noisefree. However, generally both matrix A and vector b are noiseperturbed (in Prony’s method, A and b share the same data source, see below) and, in this case, the total leastsquares technique (TLS) [3] can be more advantageous.
In some cases, a problem with the Prony polynomial method is that it can be numerically unstable because of the steps that comprise the algorithm: solving an illconditioned matrix equation and finding the roots of a polynomial. When the number of exponentials is relatively high, the sensitivity of roots of the characteristic polynomial to perturbations of their coefficient is likewise high [4] and Prony’s method may be unstable.
Another alternative is to use the matrix pencil method (MPM). Although similar to Prony’s method, it consists of solving an eigenvalue problem rather than following the conventional twostep Prony method. It has been found through perturbation analysis and simulation that for signals with unknown damping factors the MPM is less sensitive to noise than the polynomial method [5].
In recent years, and due to advances in computing systems, Prony’s method has been successfully applied in various engineering sectors, such as electric power quality analysis [6], materials science [7], antennae [8], etc. In the biomedical field, the classic Prony method is used in [9] to process multifocal visualevoked potentials (mfVEPs) to diagnose the early stages of multiple sclerosis (MS). The LS Prony method is used in [10] to estimate the parameters of the single eventrelated potential; the TLS is used in [11] to discriminate between three cardiac problems, and the MPM is used in [12,13,14].
Various programming languages are widely used in the scientific field. These languages include Python, a free and opensource highlevel programming language [15, 16], and MATLAB®, a proprietary product.
MATLAB® is userfriendly and needs practically no formal programming knowledge [17]. As it implements a wide number and variety of functions (statistics, neural networks, graphics, wavelets, etc.), it is widely accepted as a development platform for numerical software by a significant portion of the computational science and engineering community [18,19,20]. Its open availability ensures reproducibility and knowledge exchange.
Objectives
This paper presents a tutorial on implementation in MATLAB of two families of Prony methods: the polynomial method (classic and extended — LS and TLS) and the matrix pencil method. It presents an overview of the mathematical bases of each method and implements them in MATLAB using the functions directly available. The results produced by the different methods when approximating synthetic signals are obtained and filtering of mfVEP records is implemented in the Prony domain. The Discussion section provides information on possible ways of mitigating the illconditioning problems associated with several of the resolution phases of the Prony methods.
Implementation
Polynomial method
A data sequence x[n] (n = 1,…N) can be represented by the sum of p complex parameters (order p) according to the following expression:
Approximation of signal x[n] occurs in p components, in which A_{k} is the initial amplitude in the same units as x[n], α_{k} is the damping factor in seconds^{−1}, f_{k} is the frequency in Hertz, T_{S} is the sampling period (in seconds) of signal x[n] and θ_{k} is the initial phase in radians. Therefore, signal x[n] is characterized by the parameters A_{k}, α_{k}, f_{k} and θ_{k} (k = 1,…,p). h_{k} is the timeindependent component and z_{k} is an exponential and timedependent component (poles).
Equation 1 is the expression of the general solution of a homogeneous linear difference equation, if the p roots are different [21]. In order to find that equation we have to construct its characteristic equation, which is
where the characteristic roots are the parameters z_{k} in Eq. 1.
Demonstration of the Prony approximation method is found in [22]. Practical implementation requires performance of the following steps:
Step 1: Solve the linear prediction model constructed by the observed dataset and the obtained coefficients (a [1]…a[p]) of the characteristic polynomial. An autoregressive model of order p considers that the value of x[n] depends linearly on the preceding p values in x. Equation 1 can be rewritten as a linear prediction model according to the matrix system T_{pxp}.a_{px1} = − x_{px1}:
Where.
a: Linear prediction coefficients vector.
x: Observation vector.
T: Forward linear prediction matrix (a square Toeplitz matrix).
Solving this linear system (3) reveals that the values of a are the coefficients of the characteristic or Prony polynomial φ(z).
Step 2: Find the roots of the characteristic or Prony polynomial formed from the linear prediction coefficients.
Solving an equation in finite differences is achieved by finding the roots of the characteristic polynomial. As vector a is known from (3), the roots z_{k} of the polynomial φ(z) can be computed to obtain the damping factor (α_{k}) and frequency (f_{k}).
Step 3: Solve the original set of linear equations to yield the estimates of the exponential amplitude and sinusoidal phase.
First, the initial system of equations (Z_{pxp}.h_{px1} = x_{px1}) is solved:
The h_{k} values yield the amplitude (A_{k}) and phase (θ_{k}):
The classic Prony method (N = 2p) obtains an exact fit between the sampled data points and the exponentials if matrices T and Z are nonsingular. However, in many practical cases N > 2p and, in this situation, both systems are overdetermined (more equations than unknowns) and can be approximated using the LS or TLS methods.
Least squares
In general, given the overdetermined linear system: A x ≈ b with A ∈ ℂ^{mxn}, b ∈ ℂ^{mx1}, x ∈ ℂ^{nx1}, m > n; being A the data matrix and b the observation vector, the least squares solution x_{LS} is given by the normal equation:
H represents the Hermitian conjugate of a matrix and A^{+} is the Moore–Penrose pseudoinverse matrix of A. In practice, the normal equation is rarely used, as methods based on QR decomposition or singular value decomposition (SVD), among others, are preferable.
Total least squares
Solution of the system A x ≈ b by the total leastsquares method is a generalization of the LS approximation method when the data matrix A and observation vector b are contaminated with noise. In Prony’s method, eqs. 3 and 6 are constructed from the measured signals. The basic total leastsquares algorithm is [3]:
C ≔ [A : b], matrix A augmented (expansion by columns) by vector b (C ∈ ℂ^{mx(n + 1)}). SVD of C matrix is then performed:
The matrices U_{m × m} (left singular vector matrix) and V_{(n + 1) × (n + 1)} (right singular vector matrix) are orthonormal (U^{H}U = UU^{H} = I_{m}, V^{H}V = VV^{H} = I_{n + 1}) and Σ_{m × (n + 1)} = diag(σ_{1}, σ_{2}, …σ_{min {m, n + 1}})) being σ_{1} ≥ σ_{2}… ≥ σ_{min {m, n + 1}} the singular values of C.
The structure of V is as follows:
The TLS solution exists if v_{(n + 1), ( n + 1)} ≠ 0 [23] and, moreover it is unique if σ_{n} ≠ σ_{n + 1}:
Algorithms in which the solution does not exist or is not unique are considered in detail in [24].
Implementation in MATLAB of the polynomial method
The code presented was developed and tested under MATLAB R2016b. Code 1 presents implementation in MATLAB of a function to perform the Prony approximation using the three polynomial methods mentioned above. The function is defined as follows:
function [Amp,alpha,freq,theta] = polynomial_method (x,p,Ts,method)
The sampled data are given in vector x; p is the number of terms to obtain in the approximation, Ts is the sampling time of the signal and classic, LS or TLS indicates the method used to solve the problem. The function returns the parameters (Amp, alpha, freq, theta) resulting from the approximation.
First, the sample length is obtained (N = length(x)) and consistency between the parameter method, p and the sample data length is checked.
Step 1.
Coding the linear system of Eq. 2 takes into account that the MATLAB function T = toeplitz(c,r) creates nonsymmetrical Toeplitz matrix T (dimensions p × p under the classic method and (N − p) × p under the overdetermined methods), having c as its first column and r as its first row, achieved by the following instruction:
T = toeplitz (x(p:N1), x(p:1:1));
The solution of this system of eqs. (T.a = −x) for the classic and LS methods is obtained in MATLAB using the matrix left division (also known as backslash) operator. If T is square and if it is invertible, the backslash operator solves the linear equations using the QR method. With an overdetermined system, LS should be used. The backslash operator is a collection of algorithms used to solve a linear system [25], selected according to the characteristics of matrix T. Taking into account that vector x is a matrix column:
a = − T \ x(p + 1:N);
In the case of the TLS option, the function a = tls(T,x(p + 1:N)); is called (Code 2).
Step 2.
The p roots of the polynomial are now obtained:
The MATLAB instruction r = roots(c) returns a column vector whose elements are the roots of the polynomial c. Row vector c contains the coefficients of a polynomial, ordered in descending powers. If c has n + 1 components, the polynomial it represents is c_{1}s^{n} + … + c_{n}s + c_{n + 1}.
The input vector for the roots function must be a row vector and must contain the element a[0] = 1, which was not obtained in the previous solution. Its implementation is therefore.
c = transpose([1; a]);
r = roots(c);
Based on r, and having defined the acquisition period Ts, it is possible to find the values of the damping factor (α_{k}) and frequency (f_{k}):
alpha = log(abs(r))/Ts;
freq = atan2(imag(r),real(r))/(2*pi*Ts);
log is the Napierian logarithm and atan2 returns the fourquadrant inverse tangent.
Step 3: Obtain complex parameters h _{ k } from roots z _{ k } .
The number of equations (len_vandermonde) employed for the solution is set (p in classic and N in overdetermined systems) and the data matrix for the system of equations is constructed (6):
Z = zeros(len_vandermonde,p);
for i = 1:length(r).
Z(:,i) = transpose(r(i).^(0:len_vandermonde1));
End
Finally, the following is solved:
h = Z \ x(1:len_vandermonde);
In the case of the TLS option, the function h = tls(Z, x(1: len_vandermonde)); (Code 2) is called. In the TLS algorithm, SVD is used. The infinite values therefore have to be converted into maximum representative values beforehand otherwise the SVD function will yield an error.
The solutions yield the initial amplitude (A_{k}) and initial phase (θ_{k}) values:
Amp = abs(h);
theta = atan2(imag(h),real(h));
The function that solves a linear system using the TLS method (Code 2) receives as arguments matrices A and b, which define the linear system to solve: Function x = tls(A,b). The number of columns in matrix A is obtained ([~,n] = size(A);) and augmented matrix C (C = [A b]) is constructed while matrix V of the SVD decomposition is obtained via the instruction [~,~,V] = svd(C); the TLS solution (if it exists) is obtained by applying the formula (12) to matrix V.
Matrix pencil method
Steps 1 and 2 of the polynomial method yield the roots of the characteristic polynomial that coincide with the signal poles z_{k}. An alternative solution is to use the MPM to find z_{k} directly by solving a generalized eigenvalue problem.
In general, given two matrices Y_{1} ∈ ℂ^{mxn}, Y_{2} ∈ ℂ^{mxn}, the set of matrices of the form Y_{2} − λY_{1} (λ ∈ ℂ) is a matrix pencil [26].
In our case, to implement MPM a rectangular Hankel matrix Y is formed from the signal (x[n], n = 1,…N), where, in this method, p is the pencil parameter:
This matrix is used to create matrices Y_{1} and Y_{2}. Y_{1} is constructed by eliminating the last column of Y while Y_{2} is constructed by eliminating the first column of Y:
Where M is the real number of poles of function x[n], if M ≤ p ≤ (N − M) is fulfilled, the poles z_{k} (k = 1,….p) are the generalized eigenvalues of the matrix pencil Y_{2} − λY_{1} [27]; matrices Y_{1} and Y_{2} are illconditioned and therefore the QZalgorithm is not stable enough to yield the generalized eigenvalues [5]. It is more efficient to obtain the values of z_{k} from the following expression:
Where \( {\mathbf{Y}}_1^{+} \) is the Moore–Penrose pseudoinverse matrix of Y_{1}, defined as:
The values z_{k} yield the parameters α_{k} and frequency f_{k} (Equations 5 and 6); The final step coincides with Step 3 of the Prony polynomial method: solving the system Z_{pxp}.h_{px1} = x_{px1} and obtaining A_{k} and θ_{k} (Equations 8 and 9).
Coding of the MPM in MATLAB is done in Code 3, the function call being.
Function [Amp,alpha,freq,theta] = matrix_pencil (x,p,Ts)
The first step is to obtain the matrix Y then, based on that, matrices Y_{1} and Y_{2}. To achieve this, the following instruction is employed:
Y = hankel (x(1:endp), x(endp:end));
To obtain Y_{1}, the last column is eliminated.
Y1 = Y (:,1:end1);
To obtain Y_{2}, the first column is eliminated.
Y2 = Y (:,2:end);
The eigenvalues are obtained (Equation 16).
l = eig (pinv(Y1)*Y2);
eig (A) is a function that returns the eigenvalues of A while pinv(A) yields the Moore–Penrose pseudoinverse matrix of A which, in this case, corresponds to the expression in Equation 17.
The frequency (f_{k}) and damping factor (α_{k}) values are obtained from the eigenvalues in the same way as the roots are obtained in the polynomial method:
alpha = log(abs(l))/Ts;
freq = atan2(imag(l),real(l))/(2*pi*Ts);
To calculate the initial amplitude and phase values (A_{k} and θ_{k}), the steps followed are exactly the same as in the polynomial method.
Results
The methods described are applied in two situations: a) approximation of synthetic signals and b) filtering of mfVEP signals.
Synthetic functions
1 000 Functions are generated (g_{i}[n]) with N = 1 024 points each (i = 1, …1 000; n = 0, …1 023), according to the following expression
The parameters of the functions have a uniform random distribution at the following intervals: A_{k} ∈ [1, 10]; α_{k} ∈ [0, −4], f_{k} ∈ [1, 31], f_{i} ≠ f_{j}; θ_{k} ∈ [−π, π] and f_{0} = 0.
Due to the possible existence of numerical errors in the computational approximation of the functions it is advisable to evaluate the error between the original function (g_{i}[n]) and its approximation (\( \overset{\sim }{g_i} \)) using Prony’s method. The precision of the approximation obtained from the normalized rootmeansquare error is used:
‖.‖ indicates the 2norm and \( \overline{\ {g}_i} \) is the mean of the reference signal.
If for a certain function G ≥ 0.60 is fulfilled, the approximation is considered correct. Table 1 shows the number of functions correctly approximated by the Prony LS, Prony TLS and MPM methods and for the two different parameters (N, p).
None of the methods implemented works 100% correctly (G ≥ 0.60 for the 1000 g_{i}[n] functions in all the situations tested). If the mean number of functions wellapproximated by each method is considered, the best result is obtained with MPM (\( \overline{MPM}=999.55 \)) and the worst is obtained with TLS (\( \overline{TLS}=677.39 \)). The LS method yields the correct approximation in 60.52% of cases, the TLS method in 2.63% of cases and the MPM method in 92.10% of cases tested in this experiment.
In general, the results obtained using LS and MPM are very similar, as the MATLAB roots(·) function generates the companion matrix of the polynomial and uses the QRalgorithm to obtain its eigenvalues.
Figure 1 shows the roots obtained using the LS and MPM methods for one of the g_{i}[n] signals (N = 256, p = 30). The correct number of roots for signal g_{i}[n] is M = 19; in both examples, p = 30 roots are obtained, though with the MPM method 12 roots are equal to 0. This is because in the LS method the range of the companion matrix is always equal to p and, consequently, p roots are obtained. In the MPM method, the range of matrix \( \left({\mathbf{Y}}_1^{+}{\mathbf{Y}}_2\right) \) is less than or equal to p (\( \mathrm{r}=\operatorname{rank}\left({\mathbf{Y}}_1^{+}{\mathbf{Y}}_2\right)\le p \)) and r roots other than zero and (pr) roots equal to 0 are obtained [5]. In the example shown, r = 18 is fulfilled. The differences in the results between the two methods are evident in Step 3 and are due to computational errors.
mfVEP filtering
The mfVEP technique [28] can be used to obtain the electrophysiological response of the primary visual cortex to stimuli produced in a large number (e.g. 60) of sectors of the visual field. Generation of the visual stimulus is governed by the same pseudorandom sequence [29] used to separate the individual responses of each sector from the continual EEG recording obtained using electrodes. Analysis of mfVEP signals is employed in diagnosis and study of patients with glaucoma, amblyopia, nerve optic drusses, optic neuritis, multiple sclerosis and other pathologies.
The aim of this test is to evaluate whether mfVEP signal filtering in the Prony domain improves the separation between the signals of control subjects and the signals of patients with MS. It uses the signaltonoise ratio (SNR) of the records as the parameter. The discrimination factor is evaluated using the area under the ROC curve (AUC). The results achieved by applying the conventional method to mfVEP records are then compared: signals filtered using the discrete Fourier transform (DFT) between 0 and 35 Hz and the signals filtered in the Prony domain.
This experiment uses a database of mfVEP signals captured from 28 patients (age 34.39 ± 10.09 years, 7 males and 21 females) diagnosed with MS according to the McDonald criteria; the signals were obtained from 44 eyes in 22 control subjects (age 30.20 ± 7.55 years, 10 males and 12 females) with normal ophthalmologic and neurological examination results. The study protocol adhered to the tenets of the Declaration of Helsinki and was approved by the local Institutional Review Board (Comité de Ética en Investigación Clínica del Hospital Universitario Príncipe de Asturias, Alcalá de Henares, Spain). Written informed consent was obtained from all participants.
mfVEP signals were recorded monocularly with VERIS software 5.9 (ElectroDiagnostic Imaging, Inc., Redwood City, CA). The visual stimulus was a scaled dartboard with a diameter of 44.5 degrees, containing 60 sectors, each with 16 alternating checks. The luminance for the white and black checks were 200 and < 3 cd/m^{2}, respectively. The checks in each sector were reversed in contrast using a pseudorandom sequence (msequence) at a frame rate of 75 Hz.
The mfVEP signals were obtained using gold cup electrodes (impedance < 2 KΩ). The reference electrode was positioned on the inion (E_{R}) and the ground electrode on the forehead. The active electrodes were placed 4 cm above the inion (E_{A}) and 1 cm above and 4 cm either side of the inion (E_{B}, E_{C}). The difference between the signals of the active electrodes was used to obtain channels CH_{1} = E_{A}E_{R}, CH_{2} = E_{B}E_{R} and CH_{3} = E_{C}E_{R}. Three additional derived channels were obtained (CH_{4} = CH_{1}CH_{2}, CH_{5} = CH_{1}CH_{3}, CH_{6} = CH_{2}CH_{3}). Therefore, the data from 6 channels were processed. In the analogue phase, the signals were amplified at a gain of 10^{5} at a bandwidth between 3 and 100 Hz. The sampling frequency was 1200 Hz, obtaining 600 samples in each recording (length 500 ms).
The conventional signalprocessing method consists of bandpass filtering between 0 and 35 Hz using the fast Fourier transform; these signals are denominated X_{DFT}.
One method for determining the intensity of the mfVEP records is to use the signaltonoise ratio, defined by the following expression:
In an mfVEP, the physiological response to the stimulus presents in the 45–150 ms interval (signal window) following onset. In the 325–430 ms interval (noise window) only noise is considered to be recorded. RMS_{45–150 ms} and RMS_{325 − 430 ms} are the root mean square (RMS) amplitudes in the signal window and noise window, respectively.
Signal processing using Prony’s method is carried out in the following steps:
1. The Prony approximation is obtained (X_{LS}, X_{TLS}, X_{MPM}, with p = 250, N = 600) for the X_{DFT} signals. The number of MS signals is N_{MS} = 20,160 (28 × 2 × 60 × 6) and the number of control signals is N_{CONTROL} = 15,840 (22 × 2 × 60 × 6).
2. Correct approximation of the X_{DFT} signal is checked against the expression shown in Equation 19 and considering G ≥ 0.45. Figure 2 shows an example of a signal approximated using the LS method.
3. The correctly approximated signals are bandpassfiltered in the Prony domain, selecting the 10 lowestfrequency components. The MATLAB code is shown in Code 4. Figure 3 shows an example of a filtered signal.
4. The SNR value of the X_{DFT} and Pronyfiltered signals (X_{LS_F}, X_{TLS_F}, X_{MPM_F}) is obtained and the discrimination value between the signals of subjects with MS and control subjects is calculated.
Similar to the case of the synthetic signals, the LS method only correctly approximated a low percentage of records (48.79% of the control records and 42.90% of the MS records) (Table 2). The LS and MPM methods yielded the same results, achieving correct approximation in over 99% of cases. The signal intensity value in the control subjects is greater than in the subjects with MS. Filtering the signals using the conventional method yields an AUC value of 0.6538; using the TLS method yields practically the same result (AUC = 0.6472) while using the LS and MPM methods yields a value of 0.7055. This example shows that filtering in the Prony domain can increase the capacity to discriminate between the signals of control subjects and those of patients with MS.
Discussion
In this paper we have used general MATLAB functions to implement the principal methods of function approximation based on the linear combination of exponentials: the polynomial Prony method (classic, LS and TLS) and the matrix pencil method. In the polynomial method, signal poles (frequencies and damping factors) are found as roots of a polynomial while the MPM obtains the poles by finding the eigenvalues of a matrix pencil.
Currently, the most common method is Fourier analysis, which represents a signal as a summation of continuous undamped sinusoidal functions with frequency and integer times the fundamental frequency (harmonics). In contrast, the p components of a Prony series may be complex exponentials. In general, the Prony spectrum will be nonuniformly spaced in the frequency scale (as it is one of the estimated parameters), depending on the observed data [30].
Prony modelling produces higher frequency resolution than DFT methods due to its reliance on autoregressive modelling [31]. Another advantage is that it is a natural transformation for impulse responses since it uses damped sinusoids as a basis and therefore representation is efficient in terms of the number of coefficients required [32].
Not all mathematical signals can be approximated using Prony’s method [33] and computational finite arithmetic also generates errors. Consequently, the results of computational implementation of the Prony methods depend on the characteristics and number of points of the signal to be interpolated, on the p number of functions and on the use of computational procedures not susceptible to illconditioning issues. Furthermore, these potentially illconditioned operations are concatenated, thereby increasing the instability issues. For example, since the second step of Prony’s method is an illconditioned problem and roundoff errors must exist for the linear prediction parameters to be computed in the first step, the estimation of z_{k} in the second step of Prony’s method can contain significant error [34].
In our experimental approximation of synthetic functions, the best result was obtained using the MPM and LS methods, while the effectiveness of the TLS method was shown to be highly dependent on the number of points and on the p number of functions (Table 1). In some cases, when the number of exponentials is relatively high, the sensitivity of roots of the characteristic polynomial to perturbations of their coefficient is likewise high [4] and Prony’s method may be unstable.
In a second experiment, we lowpassfiltered mfVEP signals in the Prony domain in order to evaluate the improvement in the capacity to discriminate between signals of control subjects and those of MS patients. Selecting the first 10 components of each record reveals that the AUC value between the signals of healthy subjects and those of MS subjects increases by between 0.3% and 4.7% depending on the method compared. The smallest improvement was obtained with the TLS method and the greatest improvement with the LS and MPM methods.
Coding in MATLAB used the functions directly available in this programming language. However, these evidently have their computational limitations and could be replaced with better alternatives. Various aspects that could improve the code presented in this paper are discussed below.
Solution of linear systems
Solution of the linear systems using the classic and LS methods was implemented with the MATLAB mldivide (\) operator. Although the mldivide operator is valid for most cases (it selects between the LU, Cholesky, LDLT or QRfactorization methods, among others, depending on the characteristics of matrix A [35]), it may be more efficient to implement other algorithms.
The numerical stability of the solution in linear algebra may be evaluated by the condition number and the numerical rank of matrix A. The condition number is defined as: \( {k}_2\left(\mathbf{A}\right)=\frac{\sigma_{max}}{\sigma_{min}} \); a low condition number usually means that the system is wellconditioned. The rank (r) of a matrix is the number of linear independent rows (or columns) (r ≤ min {m, n}) and is equal to the number of singular values (σ_{i}) in the matrix other than zero. When r = min {m, n} the matrix has full range, otherwise it is rankdeficient. If A is nearly rankdeficient (σ_{min} is small), then the solution x is illconditioned and possibly very large. A more robust solution to obtain the effective rank may be to evaluate the number of singular values of AA^{H} or A^{H}A above a specified tolerance. Analysing the condition number and the rank of a matrix may make it possible to select the best method for system solution.
Least squares
In general, although the normal equation is the fastest method it is not used to solve systems by LS as it yields worse numerical results than other methods. In the normal equation, accuracy depends on \( {k}_2\left({\mathbf{AA}}^{\mathbf{H}}\right)={k}_2^2\left(\mathbf{A}\right) \), although this method may be used if A is wellconditioned [36]. If A is rankdeficient, then x = A\B is not necessarily the minimum norm solution. The more computationally expensive x = pinv(A)*B computes the minimum norm leastsquares solution. Specifically, the function pinv(A, tol) returns the Moore–Penrose pseudoinverse, obtained by SVD decomposition where the values above tolerance (tol) are set to zero; this may be adapted to an illconditioned problem (A is not of full rank). Another option to obtain the Moore–Penrose pseudoinverse is proposed in [37], which makes use of QRfactorization and an algorithm based on a reverse order law for generalized inverse matrices; this method was later refined in [38]. An iterative solution to obtain the Moore–Penrose pseudoinverse was published in [39].
TLS
The TLS method implemented is the one that yielded the worst computational performance. This method performs SVD of the augmented matrix C ≔ [A : b]; If C is rankdeficient or nearly rankdeficient (its singular values decay gradually), it may be advisable to truncate its small singular values [40]. [41] presents basic information, references and applications for the TLS method. In [42], an interactive method is proposed which combines orthogonal projections to a sequence of generalized Krylov subspaces of increasing dimensions and Newton’s method. The introduction to [43] presents various alternatives to obtaining the solution using the TLS method and the authors present a solution based on randomized algorithms.
Roots
Numerical solution of a polynomial is a classic problem in mathematical research [44]. Methods available with which to obtain the roots of a polynomial include Laguerre [45], Bairstow, Græffe and Müller, Horner, Jenkins and Traub, and Newton [46], etc., with differing performances in terms of accuracy, convergence and speed. The code presented uses the roots() function used by the QRalgorithm on the balanced Frobenius companion matrix to compute its eigenvalues.
Eigenvalues and SVD
The eigenvalues of a square matrix A are the roots of its characteristic polynomial det(A − λI) = 0. However, singular values of A are nonnegative square roots of eigenvalues of (A^{T}A), meaning that both methods are related. The general idea is to diagonalize the target matrix as the values of the diagonal are the eigenvalues. All methods to solve the eigenvalue problem are of an iterative nature [47]. The builtin MATLAB function eig(A) uses the generalized Schur decomposition method (implemented via the QRalgorithm or its variants), which consists of interactively obtaining an upper triangular matrix U, in which the values of its diagonal are the eigenvalues of A. The QRalgorithm can be adapted to small or moderately large nonsymmetrical eigenvalue problems [48]. For large matrices, [49] provides possible alternatives.
Pronylike methods
Other modifications have been made to the Prony method, generally with the intention of improving its numerical stability. If any of the parameters of equation (1) are known, the Prony method makes it easier to find a robust solution. In the polynomial method, small variations in the coefficients of equation (2) due to signal noise can result in large variations in its zeros and, consequently, the frequencies of the approximation will vary greatly. Parametric spectral estimation techniques, such as MUSIC (MUltiple SIgnal Classification) [50], ESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques) [51] or fast ESPRIT [52] offer an alternative that in many cases make it possible to obtain more robust solutions. [53] presents an algorithm for the factorization of a matrix pencil based on QR decomposition of a rectangular Hankel matrix, which simplifies the ESPRIT method.
The NAFASS (Nonorthogonal Amplitude Frequency Analysis of the Smoothed Signals) approach [54] makes it possible to obtain the set of frequencies that make up strongly correlated random sequences with N < 1500. [55] presents the physical interpretation of the Prony decomposition as a linear recording of memory effects that can exist in random sequences in which the Fourier decomposition is a partial case. [56] improves the NAFASS method, presenting a linear recurrence expression that obtains the set of frequencies.
Another way to obtain more robust results is to act on the signals before obtaining their decomposition in a Prony series by using prefiltering [57]. In the modified instantaneous Prony method [58] the input data used in an application of speech parameter extraction are those derived from the signal x[n] instead of adjacent samples.
Applying the Prony method to a time window that can be moved along the x[n] signal makes it possible to perform time–frequency analysis. One such example could be the shorttime matrix pencil method (STMPS) successfully used to obtain antennae responses [59]. The Piecewise Prony Method [60] essentially consists of dividing the signal to be interpolated into windows of variable length and sampling rate.
Conclusions
Decomposition of a signal using Prony’s method can be considered a generalization of the Fourier decomposition. Although the method has been known since 1795, its application in engineering has only increased since about the 1970s as computer use has grown. This paper has presented the theoretical bases and their piecebypiece implementation in MATLAB. It has also shown some of their limitations and the benefit of improving the quality of the mfVEP signals. With the information provided, readers can begin practical implementation of the most common Prony methods, test the reliability of the results and, if applicable, research other methods more appropriate to their areas of research.
Abbreviations
 AUC:

Area under the ROC curve
 DFT:

Discrete Fourier transform
 LS:

Least squares
 mfVEP:

Multifocal visualevoked potential
 MPM:

Matrix pencil method
 MS:

Multiple sclerosis
 RMS:

Root mean square
 ROC:

Receiveroperatingcharacteristic
 SNR:

Signaltonoise ratio
 SVD:

Singular value decomposition
 TLS:

Total least squares
References
 1.
Prony R. Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques, et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alkool, á différentes temperatures. J L’école Polytech. 1795;1:24–76.
 2.
Householder, AS. On Prony’s method of fitting exponential decay curves and multiplehit survival curves. Oak Ridge National Laborator, 1950.
 3.
Markovsky I, Van Huffel S. Overview of total leastsquares methods. Signal Process. 2007;87:2283–302. https://doi.org/10.1016/j.sigpro.2007.04.004.
 4.
Guillaume P, Schoukens J, Pintelon R. Sensitivity of roots to errors in the coefficient of polynomials obtained by frequencydomain estimation methods. IEEE Trans Instrum Meas. 1989;38:1050–6. https://doi.org/10.1109/19.46399.
 5.
Hua Y, Sarkar TK. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise. IEEE Trans Acoust. 1990;38:814–24. https://doi.org/10.1109/29.56027.
 6.
Lobos T, Leonowicz Z, Rezmer J, Schegner P. Highresolution spectrumestimation methods for signal analysis in power systems. IEEE Trans Instrum Meas. 2006;55:219–25. https://doi.org/10.1109/TIM.2005.862015.
 7.
Park SW, Schapery RA. Methods of interconversion between linear viscoelastic material functions. Part I—a numerical method based on Prony series. Int J Solids Struct. 1999;36:1653–75. https://doi.org/10.1016/S00207683(98)000559.
 8.
Sarkar T, Weiner D, Jain V. Some mathematical considerations in dealing with the inverse problem. IEEE Trans Antennas Propag. 1981;29:373–9. https://doi.org/10.1109/TAP.1981.1142573.
 9.
Fernández A, de Santiago L, Blanco R, PérezRico C, RodríguezAscariz JM, Barea R, MiguelJiménez JM, GarcíaLuque JR, Ortiz del Castillo M, SánchezMorla EM, Boquete L. Filtering multifocal VEP signals using Prony’s method. Comput Biol Med. 2015;56:13–9. https://doi.org/10.1016/j.compbiomed.2014.10.023.
 10.
Hansson M, Gänsler T, Salomonsson G. Estimation of single eventrelated potentials utilizing the Prony method. IEEE Trans Biomed Eng. 1996;43:973–81. https://doi.org/10.1109/10.536898.
 11.
Chen SW. A twostage discrimination of cardiac arrhythmias using a total least squaresbased Prony modeling algorithm. IEEE Trans Biomed Eng. 2000;47:1317–27. https://doi.org/10.1109/10.871404.
 12.
Bhuiyan M, Malyarenko EV, Pantea MA, Seviaryn FM, Maev RG. Advantages and limitations of using matrix pencil method for the modal analysis of medical percussion signals. IEEE Trans Biomed Eng. 2013;60:417–26. https://doi.org/10.1109/TBME.2012.2227318.
 13.
Bauman G, Bieri O. Matrix pencil decomposition of timeresolved proton MRI for robust and improved assessment of pulmonary ventilation and perfusion. Magn Reson Med. 2017;77:336–42. https://doi.org/10.1002/mrm.26096.
 14.
Gopalakrishnan R, Machado AG, Burgess RC, Mosher JC. The use of contact heat evoked potential stimulator (CHEPS) in magnetoencephalography for pain research. J Neurosci Methods. 2013;220:55–63. https://doi.org/10.1016/j.jneumeth.2013.08.015.
 15.
Goodman D. Brian: a simulator for spiking neural networks in Python. Front Neuroinform. 2008;2:5. https://doi.org/10.3389/neuro.11.005.2008.
 16.
Meyer R, Obermayer K. Pypet: a Python toolkit for data management of parameter explorations. Front. Neuroinform. 2016;10. https://doi.org/10.3389/fninf.2016.00038.
 17.
Sen SK, Shaykhian GA. MatLab tutorial for scientific and engineering computations. Nonlinear Anal Theory, Methods Appl. 2009:e1005–20. https://doi.org/10.1016/j.na.2009.01.069.
 18.
Ihlen EAF. Introduction to multifractal detrended fluctuation analysis in Matlab. Front Physiol. 2012;3. https://doi.org/10.3389/fphys.2012.00141.
 19.
MiguelJiménez JM, Boquete L, Ortega S, Cordero CA, Barea R, Blanco R. mfERG_LAB: software for processing multifocal electroretinography signals. Comput Methods Prog Biomed. 2012;108:377–87. https://doi.org/10.1016/j.cmpb.2012.02.013.
 20.
Perakakis P, Joffily M, Taylor M, Guerra P, Vila J. KARDIA: a Matlab software for the analysis of cardiac interbeat intervals. Comput Methods Prog Biomed. 2010;98:83–9. https://doi.org/10.1016/j.cmpb.2009.10.002.
 21.
Saber E. An introduction to difference equations. New York: SpringerVerlag; 2005. https://doi.org/10.1007/0387276025.
 22.
Reddy DC. Biomedical signal processing: principles and techniques. New York: McGrawHill; 2005.
 23.
Van Huffer S, Zha H. The total least squares problem. Handb. Stat. 1993;9:377–408.
 24.
Van Huffel S, Vandewalle J. The total least squars problem: computational aspects and analysis. Philadelphia: Siam; 1991.
 25.
Yang WY, Cao W, Chung TS, Morris J. Applied numerical methods using MATLAB. Hoboken: Wiley; 2005.
 26.
Gantmacher FR. The theory of matrices, vol. 2. Moscow: GITTL; 1953.
 27.
Sarkar TK, Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials. IEEE Antennas Propag Mag. 1995;37:48–55. https://doi.org/10.1109/74.370583.
 28.
Baseler HA, Sutter EE, Klein SAA, Carney T. The topography of visual evoked response properties across the visual field. Electroencephalogr Clin Neurophysiol. 1994;90:65–81. https://doi.org/10.1016/00134694(94)901147.
 29.
Müller PL, Meigen T. Msequences in ophthalmic electrophysiology. J Vis. 2016;16(1):15.
 30.
Mitrofanov G, Priimenko V. Prony filtering of seismic data. Acta Geophys. 2015;63:652–78. https://doi.org/10.1515/acgeo20150012.
 31.
Lander P, Jones D, Berbari E, Lazzara R. Timefrequency structure of the highresolution ECG. J Electrocardiol. 1994;27:207–12. https://doi.org/10.1016/S00220736(94)800936.
 32.
Demiralp T, Ademoglu A, Istefanopulos Y, Gülçür HO. Analysis of eventrelated potentials (ERP) by damped sinusoids. Biol Cybern. 1998;78:487–93. https://doi.org/10.1007/s004220050.
 33.
Slivinskas V, Šimonyte V. On the foundation of Prony’s method. In: Stoch. Control, Elsevier; 1987. p. 121–6. https://doi.org/10.1016/B9780080334523.500259.
 34.
James Hu SL, Yang WL, Li HJ. Signal decomposition and reconstruction using complex exponential models. Mech Syst Signal Process. 2013;40:421–38. https://doi.org/10.1016/j.ymssp.2013.06.037.
 35.
Davis TA, Duff IS. An unsymmetricpattern multifrontal method for sparse LU factorization. SIAM J Matrix Anal Appl. 1997;18:140–58. https://doi.org/10.1137/S0895479894246905.
 36.
Demmel JW. Applied numerical linear algebra, SIAM; 1997. https://doi.org/10.1137/1.9781611971446.
 37.
Katsikis VN, Pappas D, Petralias A. An improved method for the computation of the Moore–Penrose inverse matrix. Appl Math Comput. 2011;217:9828–34. https://doi.org/10.1016/j.amc.2011.04.080.
 38.
Ataei A. Improved Qrginv algorithm for computing MoorePenrose inverse matrices. ISRN Appl Math. 2014;2014:1–5. https://doi.org/10.1155/2014/641706.
 39.
Petković MD, Stanimirović PS. Two improvements of the iterative method for computing Moore–Penrose inverse based on Penrose equations. J Comput Appl Math. 2014;267:61–71. https://doi.org/10.1016/j.cam.2014.01.034.
 40.
Fierro RD, Golub GH, Hansen PC, O’Leary DP. Regularization by truncated total least squares. SIAM J Sci Comput. 1997;18:1223–41. https://doi.org/10.1137/S1064827594263837.
 41.
Markovsky I. Bibliography on total least squares and related methods. Stat Interface. 2010;3:329–34.
 42.
Lampe J, Voss H. Largescale Tikhonov regularization of total least squares. J Comput Appl Math. 2013;238:95–108. https://doi.org/10.1016/j.cam.2012.08.023.
 43.
Xie, P., Wei, Y., Xiang, H. Perturbation analysis and randomized algorithms for largescale total least squares problems. arXiv preprint arXiv:1401.6832, 2014.
 44.
Pan VY. Solving a polynomial equation: some history and recent progress. SIAM Rev. 1997;39:187–220. https://doi.org/10.1137/S0036144595288554.
 45.
Hansen E, Patrick M, Rusnak J. Some modifications of Laguerre’s method. BIT. 1997;17:409–17. https://doi.org/10.1007/BF01933450.
 46.
Madsen K. A rootfinding algorithm based on Newton’s method. BIT. 1973;13:71–5. https://doi.org/10.1007/BF01933524.
 47.
Golub GH, van der Vorst HA. Eigenvalue computation in the 20th century. J Comput Appl Math. 2000;123:35–65.
 48.
Chandrasekaran S, Gu M, Xia J, Zhu JA. Fast QR algorithm for companion matrices. In: Recent Adv. Matrix Oper. Theory, Birkhäuser Basel, Basel; 2007. p. 111–43. https://doi.org/10.1007/9783764385392_7.
 49.
Saad Y. Numerical methods for large eigenvalue problems, Society for Industrial and Applied Mathematics; 2011. https://doi.org/10.1137/1.9781611970739.
 50.
Schmidt R. Multiple emitter location and signal parameter estimation. IEEE Trans Antennas Propag. 1986;34:276–80. https://doi.org/10.1109/TAP.1986.1143830.
 51.
Roy R, Paulraj A, Kailath T. ESPRIT—A subspace rotation approach to estimation of parameters of cisoids in noise. IEEE Trans Acoust Speech Signal Process. 1986;34:1340–2. https://doi.org/10.1109/TASSP.1986.1164935.
 52.
Potts D, Tasche M. Fast ESPRIT algorithms based on partial singular value decompositions. Appl Numer Math. 2015;88:31–45. https://doi.org/10.1016/j.apnum.2014.10.003.
 53.
Potts D, Tasche M. Parameter estimation for nonincreasing exponential sums by Pronylike methods. Linear Algebra Appl. 2013;439:1024–39. https://doi.org/10.1016/j.laa.2012.10.036.
 54.
Nigmatullin RR, Osokin SI, Toboev VA. NAFASS: discrete spectroscopy of random signals. Chaos, Solitons Fractals. 2011;44:226–40. https://doi.org/10.1016/j.chaos.2011.02.003.
 55.
Nigmatullin RR, Khamzin AA, Machado JT. Detection of quasiperiodic processes in complex systems: how do we quantitatively describe their properties? Phys Scr. 2013;89:015201. https://doi.org/10.1088/00318949/89/01/015201.
 56.
Nigmatullin RR, Gubaidullin IA. NAFASS: fluctuation spectroscopy and the Prony spectrum for description of multifrequency signals in complex systems. Commun Nonlinear Sci Numer Simulat. 2018;56:252–69. https://doi.org/10.1016/j.cnsns.2017.08.009.
 57.
Kumaresan R, Feng Y. FIR prefiltering improves Prony’s method. IEEE Trans Signal Process. 1991;39:736–41. https://doi.org/10.1109/78.80860.
 58.
Azarov E, Vashkevich M, Petrovsky A. Instantaneous harmonic representation of speech using multicomponent sinusoidal excitation. Proc Annu Conf Int Speech Commun Assoc Interspeech. 2013:1697–701.
 59.
Rezaiesarlak R, Manteghi M. Shorttime matrix pencil method for chipless RFID detection applications. IEEE Trans Antennas Propag. 2013;61:2801–6. https://doi.org/10.1109/TAP.2013.2238497.
 60.
Garoosi V, Jansen BH. Development and evaluation of the piecewise Prony method for evoked potential analysis. IEEE Trans Biomed Eng. 2000;47:1549–54. https://doi.org/10.1109/10.887935.
Acknowledgements
Thanks to Dr. Román Blanco (Universidad de Alcalá, Spain), for ceding the mfVEP recordings
Availability and requirements
Project name: MAT_PRONY.
Project home page: Not applicable.
Operating system(s): Platform independent.
Programming language: Matlab.
Other requirements: Not applicable.
License: Matlab licences.
Any restrictions to use by nonacademics: Not applicable.
Funding
Universidad de Alcalá grant “Diagnóstico precoz de neuropatías ópticas mediante análisis de registros de potenciales evocados visuales multifocales” (UAH GC2016–004) and Agencia Estatal de Investigación (AEI) – European Regional Development Fund (ERDF) grant “Investigación de la técnica de potenciales evocados visuales multifocales. Aplicación en estudios de evolución de esclerosis múltiple y evaluación de medicamentos” reference: DPI2017–88438R (AEI/FEDER,UE). The funding organizations did not play any role in the design of the study, data collection and analysis, or preparation of the manuscript.
Availability of data and materials
The Matlab code is included in this published article.
Author information
Affiliations
Contributions
AFR and LB conceived and designed the experiments; LdSR and JMRA performed the experiments; ELG and JMMJ revised the mathematics concepts; LB and AFR. wrote the paper; all the authors reviewed and approved the manuscript.
.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
The mfVEP signals adquisition study protocol was approved by the Comité de Ética en Investigación Clínica del Hospital Universitario Príncipe de Asturias, Spain and adhered to the tenets of the Declaration of Helsinki. All participants provided written informed consent.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Fernández Rodríguez, A., de Santiago Rodrigo, L., López Guillén, E. et al. Coding Prony’s method in MATLAB and applying it to biomedical signal filtering. BMC Bioinformatics 19, 451 (2018). https://doi.org/10.1186/s128590182473y
Received:
Accepted:
Published:
Keywords
 Prony’s method
 Matrix pencil
 Least squares
 Total least squares
 Multifocal evoked visual potentials
 Multiple sclerosis