IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 51, NO. 2, FEBRUARY 2004 237 Identification of Hammerstein Models With Cubic Spline Nonlinearities Erika J. Dempsey, Student Member, IEEE and David T. Westwick, Member, IEEE Abstract—This paper considers the use of cubic splines, instead of polynomials, to represent the static nonlinearities in block structured models. It introduces a system identification algorithm for the Hammerstein structure, a static nonlinearity followed by a linear filter, where cubic splines represent the static nonlinearity and the linear dynamics are modeled using a finite impulse response filter. The algorithm uses a separable least squares Levenberg-Marquardt optimization to identify Hammerstein cascades whose nonlinearities are modeled by either cubic splines or polynomials. These algorithms are compared in simulation, where the effects of variations in the input spectrum and distribu- tion, and those of the measurement noise are examined. The two algorithms are used to fit Hammerstein models to stretch reflex electromyogram (EMG) data recorded from a spinal cord injured patient. The model with the cubic spline nonlinearity provides more accurate predictions of the reflex EMG than the polynomial based model, even in novel data. Index Terms—Cubic splines, mean square optimization, non- linear system identification, separable least squares, stretch reflex dynamics. I. INTRODUCTION S YSTEM identification is the process of building mathemat- ical models based on measurements of a system’s inputs and outputs [1]. These models are often used in biomedical en- gineering to gain insight into the behavior of physiological sys- tems. Some biological processes, such as the stretch reflexes in postural muscles, are highly nonlinear. Studying them requires both nonlinear models and methods to identify them. The Hammerstein cascade [2], shown in Fig. 1, consists of a memoryless nonlinearity followed by a linear dynamic system. It is less general than functional expansion models, such as Volterra and Wiener series [3], but when used appropriately it can reduce the number of parameters necessary to accurately represent a system. A Hammerstein cascade can be a useful model, especially when representing systems containing hard nonlinearities, such as stretch reflex dynamics [4], where func- tional expansions would be impractical. It has also been used to identify lung tissue dynamics [5], [6], a weakly nonlinear system whose long memory length renders functional expan- sion models impractical. Manuscript received November 28, 2002; revised May 17, 2003. This work was supported in part by the Natural Sciences and Engineering Research Council (Canada) and in part by the University of Calgary Research Grants Office. Asterisk indicates corresponding author. E. J. Dempsey is with the Department of Electrical and Computer Engi- neering, University of Calgary, Calgary, AB T2N 1N4, Canada. *D. T. Westwick is with the Department of Electrical and Computer Engi- neering, University of Calgary, 2500 University Drive NW, Calgary, AB T2N 1N4, Canada (e-mail: westwick@enel.ucalgary.ca). Digital Object Identifier 10.1109/TBME.2003.820384 Fig. 1. Block diagram of a Hammerstein structure, consisting of a static nonlinearity, followed by a dynamic linear system. The Hammerstein model is nonlinear in its parameters, so it cannot be identified with a linear regression. Several identifi- cation algorithms have been developed, including an iterative, cross-correlation based algorithm by Hunter and Korenberg [7], and a method proposed by Chang and Luus [8] that uses linear regression to identify an auxiliary model that is linear in its pa- rameters, and then reconstructs the Hammerstein cascade from it. Westwick and Kearney [9] have proposed an optimization algorithm that uses the principle of “separable least squares” (SLS) to limit the iterative search to finding the optimal param- eters for the static nonlinearity. These methods all used polynomials to represent the static nonlinearities. However, block structured models are primarily used to represent highly nonlinear systems, thus requiring the use of high-order polynomials. These often include oscillations, and their accuracy tends to degenerate at the edges of the data set, so that values obtained from extrapolating beyond the limits, and even those near the limits, are of questionable accuracy. Finally, and most significantly, the coefficients describing the polynomial do not act locally; changing one coefficient of the polynomial changes the shape of the entire nonlinearity. Fitting a polynomial to a sharp corner often results in ripples, especially in areas where the data is sparse. Because of these characteris- tics, it is unlikely that polynomials are the ideal representation for hard nonlinearities. The SLS optimization used by Westwick and Kearney [9] is not limited to polynomial nonlinearities. Indeed, any represen- tation that is differentiable with respect to its parameters can be used. Rational polynomials, Fourier series, orthogonal wavelets and artificial neural networks were suggested as possible alter- natives. This paper considers a cubic spline [10] as the nonlin- earity in a Hammerstein cascade. A cubic spline is a piecewise cubic function defined by a series of knot points. Between each pair of knots the spline is defined by a third degree polynomial. The polynomials are chosen such that the spline is a continuous, twice differentiable function. Splines have been used instead of polynomials in curve fitting applications for several reasons. Because they are defined by the positions of their knots, the parameters defining cubic splines act primarily locally. Moving one knot changes the third-order polynomials between it and the two adjacent knots, only slightly affecting the functions beyond. Unlike 0018-9294/04$20.00 © 2004 IEEE