Efficient Synthesis of Complex, Driven Systems

The present invention relates generally to synthesis of the behavior of complex physical systems, and in particular to the use of probabilistic modeling techniques to these ends.

Complex physical systems, such as musical instruments, are difficult to physically reproduce or replicate faithfully in terms of their structures, minute variations of which can result in observably different response characteristics. For such systems, it may be preferable to mathematically model the behavior of the system and sythesize the output based on the model; that is, the model becomes (or operates) as a replacement for the physical system itself.

Simple physical systems (or more complex systems treated simply for purposes of approximation) of course can be straightforwardly "modeled" through direct use of the physical parameters that govern the system's behavior; for example, the motion of a simple pendulum can be precisely characterized (ignoring air resistance) in terms of length and angular acceleration, or the vibration of a string in terms of amplitude and angular frequency. Far less tractable is the goal of constructing, for example, an electronic violin capable of emulating a Stradivarius with any degree of fidelity. Complex systems contain tremendous numbers of degrees of freedom; accordingly, "first principles" physical modeling of their responses is highly difficult, both in terms of computational requirements--to keep up with changes in the many degrees of freedom--and sufficiently precise measurements of their values. Not surprisingly, approaches to physical modeling require significant approximation. For example, European Patent Application No. 0583043 ("Tone Generation System") describes the use of digital waveguides to model the behavior of musical instruments. The modeled waveforms do not represent underlying degrees of freedom in any rigorous sense, and are not related in any physically meaningful way to system inputs.

A more computationally tractable approach involves approximation or sampling techniques, whereby, for example, a list of known inputs and observed outputs is constructed, and behavior "modeled" through lookup and interpolation among entries; see, e.g., U.S. Pat. No. 5,521,322. This approach is limited, however, by the coarseness of the entries and the constrained generality of any model unconnected with the underlying physics of the system. Thus, musical instruments have been synthesized by storing output sounds for numerous known inputs--i.e., specific player manipulations of the actual instrument--and the instrument synthesized through interpolation among the sampled sound based on provided input. In practice, while short segments of such recorded sounds can be faithful, the overall response to the player's actions is not.

A related approach, exemplified by the disclosure of U.S. Pat. No. 4,018,121 (to Chowning) is to model the output behavior of a physical system as a mathematical waveform, without reference to the underlying physical degrees of freedom. Although output modeling may reflect greater attention to actual system behavior than a mechanistic interpolation approach, the validity of the output model, once again, is ultimately limited due to the absence of any connection with the true system degrees of freedom.

The "first principles" approach can be made more tractable by recognizing that not all possible physical degrees of freedom participate in the generation of a system output; a faithful model of a Stradivarius, for example, does not require detailed knowledge at the molecular level. Thus, machine learning and state-space reconstruction have been employed as intermediate approaches between a complete but computationally unachievable specification of a physical system, on one hand, and purely numerical techniques without reference to system state on the other. State-space reconstruction by the technique of time-delay embedding permits recovery and modeling of those effective physical degrees of freedom relevant to observed behavior; that is, the part of the "configuration space" (which specifies the values of all potentially accessible physical degrees of freedom) that the system actually explores as its dynamics unfolds. In accordance with the technique of time-delay embedding, a time series is measured from such a system, and the entries used to define a lag vector. For example, if s is the state vector describing the system (in terms of effective degrees of freedom rather than the complete configuration space), ds/dt=.function.(s) denotes the effective governing equations, and the measured quantity observed over time is y=y(s(t)) (where y is a vector or a scalar quantity such as amplitude or temperature), then given a delay time .tau. and a dimension d, a lag vector x may be defined as x=(y.sub.t,y.sub.t-.tau., . . . ,y.sub.t-(d-1).tau.). The central result of time-delay embedding is that the behavior of .function.(s) and x will differ only by a smooth, local, invertible change of coordinates for almost every possible choice of s, y(s) and .tau., so long as d (the "embedding dimension," i.e., the number of time lags) is sufficiently large, y depends on at least some of the components of s, and the remaining components of s are coupled by the governing equations to the ones that influence y; this result can be generalized to linear transformations on lag vectors as well.

In other words, if an experimentally observed quantity arises from deterministic governing equations, it is possible to use time-delay embedding to recover a representation of the relevant internal degrees of freedom of the system from the observable; and because the relevant degrees of freedom typically are a relatively small subset of the configuration space, the solution, while highly accurate, is also computationally tractable. System behavior can thus be modeled based on the mapping between the lag vector, whose time-varying components are measurable, and internal (and therefore generally inaccessible) system states.

In a driven system, some user input u is imposed on the system and affects its dynamic behavior; for example, the system might be a violin and the input the player's drawing of the bow. In this case, while it would be desirable to predict the behavior of the system from the input rather than the observable or the inaccessible internal degrees of freedom, this is generally not possible. Put differently, the evolution of the system cannot be described simply as y=.function.(u), since the system's behavior depends on its prior history as well as the input; the system is said to have "memory." If x is the embedding vector--which, again, maps smoothly to internal system behavior and therefore provides a complete specification of the internal degrees of freedom relevant to the observable--then the behavior of the system can be modeled given knowledge of the embedding vector x, which includes time lags on both y and u. Accordingly, it is possible to use time-delay embedding to predict y from time lags on y and u, i.e., the history of the output and the input. For example, in the above-noted case of a violin, y is amplitude (a scalar quantity) and u is the time-varying bow position and/or fingering. With x characterized, a function y=.function.(x,u) can be derived to model and predict system behavior in terms of the unobservable degrees of freedom (modeled by x) and the observable input u.

The variables y and x are related by an unknown joint probability density p(y,x) (which for a deterministic system reduces to a prediction surface y(x), while for a stochastic system samples can be drawn from p(y,x) to emulate the system's behavior). Given a set of experimental measurements {y.sub.n,x.sub.n }.sub.n=1.sup.N, the goal of data analysis is to infer a model that can predict the output or its distribution from a measurement of a new input, and to characterize the relationship between y and x. Since there is rarely enough data to estimate the unconditional density p(y,x), it is more common to seek conditional quantities such as the likelihood y.vertline.xand the error .sigma..sub.y.sup.2 .vertline.x. Traditionally these have been written in the form ##EQU1## where the .beta..sub.m quantities are unknown linear coefficients and the f.sub.m quantities are known basis functions (such as polynomial expansions). More powerful techniques, such as neural networks, utilize coefficients inside nonlinear basis functions of the form ##EQU2##

Any phenomenological model must balance underfitting (in which the model is not flexible enough to describe the data) and overfitting (in which the model is so flexible that it describes noise in the data that does not generalize). While this has been done by varying the number M of basis functions f.sub.m, such an approach is rarely justified because globally simple behavior might require a large number of terms to be represented in a given basis.

The present invention approaches the problem of synthesis of driven systems by modeling the effective underlying degrees of freedom, preferably using a probabilistic framework. In accordance with the invention, a replica of the system's behavior in response to external stimulus is developed computationally, and the model used to replace (or facilitate replacement) of the system with, for instance, a physical representation programmed to behave in accordance with the model. For example, the invention may be applied to develop a model capturing the behavior of a complex musical instrument such as a violin; the model then may be embodied in any physically appealing format (e.g., as a plastic replica of the original violin that would, absent the implemented model, produce no sound if bowed; or a keyboard or other musical instrument whose response to being "played" is to generate the sounds of the original violin). It is important to stress that once the model of system behavior (based on effective degrees of freedom) is complete, its utilization is entirely within the user's discretion.

The preferred embodiment of the invention utilizes a flexible, probabilistic model of behavior that includes information based on prior beliefs. The simplest way to approach such an analysis is to constrain the probabilities with extra terms (called regularizers) that express prior beliefs about the data. One may use a Bayesian framework to identify the best model m as the one that maximizes the likelihood that it was generated from the measured data d: ##EQU3## The quantity p(d.vertline.m) measures the model mismatch (i.e., the fitting error, so that for a Gaussian error model it leads to the familiar least-squares statistic), p(m) is the "regularizer" and represents prior beliefs characterizing a good model, and p(d) represents prior beliefs about the likelihood of a data set (and becomes relevant only when analyzing multiple data sets).

Common regularizers include maximum entropy (which handles discontinuities well but cannot capture local smoothness) and integrated curvature (which enforces smoothness but rounds out discontinuities). A more general problem with regularization is the global nature of its operation; the statements concerning prior probabilities may not be valid locally. Regularized models also do not simplify in a natural fashion; a linear system might be described, for example, by many nonlinear basis functions plus a regularizer. Though such a nonlinear model may accurately describe a data set, it fails to offer substantial insight into the behavior of the modeled system. Finally, maximization of the overall cost function required by regularization necessitates a nonlinear search that may be extremely computationally intensive in a high-dimensional space.

In accordance with the preferred embodiment of the present invention, the joint probability density p(y,x) is expanded in clusters, each of which is associated with a local model (which may be linear or nonlinear). As a result, the system's behavior is represented locally, allowing for separate modeling of different aspects of the system's dynamics, and ultimately leading to a more accurate, simplified overall representation that avoids both underfitting and overfitting. The modeling technique of the present invention is computationally efficient, exhibits good convergence and stability properties, is capable of expressing prior knowledge about the system, and is equally suited to description of simple and complex behavior.

In a first aspect, the method of the invention is used to emulate the behavior of a nonlinear physical system that generates an output in response to a stimulus such as physical manipulation. The system may be stochastic or fully determinstic, and the term "physical system" as used herein connotes either type. In accordance with the method of the invention, an input-output data set is first recorded by imposing a plurality of input manipulations on the system and measuring the resulting outputs; for example, in the context of a musical instrument, sensors are applied to the instrument and associated with the player, and both the audio outputs and player inputs are recorded as the instrument is played arbitrarily to explore the range of possible responses. These data are used to recover the internal system degrees of freedom, e.g., using state-space reconstruction embedding techniques. The system internal degrees of freedom are related to system outputs as a prediction surface (for a fully determinstic system) or a probability density (for a stochastic system) in the joint input-output space, and a predicted system output is thereby obtained for a given manipulation by (i) estimating values for the internal degrees of freedom that the manipulation would cause if applied to the system, and (ii) based on the modeled relationship between internal degrees of freedom and system output, computing an expected system output. For a determinstic system, the expected output is single-valued; for a stochastic system, the expected output is represented as a distribution that may be randomly sampled.

The approach of the invention is distinct from prior-art systems, such as those that merely sample system output or calculate a complete representation of system physics; in effect, the present invention samples the physics.

In a second aspect, the invention comprises a method of modeling one or more output characteristics of a system based on a set of input parameters. An emulation process relating each output characteristic to the input parameters is modeled according to a joint probability density therebetween by expanding the joint probability density in a plurality of clusters, each cluster (i) being associated with a local model relating the ouput characteristic to the input parameters and (ii) valid over a range of values of the input parameters according to a probability distribution, the joint probability density for an input set of values for the input parameters being a weighted sum of the clusters in accordance with the probability distributions thereof. Using this model, a set of input values is developed for the input parameters, and the emulation process used to produce the output characteristic.

In a third aspect, the invention comprises hardware apparatus for implementing the foregoing.

The foregoing discussion will be understood more readily from the following detailed description of the invention, when taken in conjunction with the accompanying drawings, in which:

FIG. 1 graphically depicts the behavior of a driven, dynamic system and its representation in accordance with the present invention;

FIG. 2 graphically depicts the relationship between the embedding dimension and expected error;

FIG. 3 graphically depicts the relationship between the number of clusters and error;

FIG. 4 is a flow chart illustrating operation of the method of the invention;

FIG. 5 schematically depicts a representative hardware implementation of the invention; and

FIGS. 6A and 6B graphically illustrate application of the invention to the behavior of a bowed string.

1. Mathematical Framework

A physical system responsive to external manipulation can be modeled as a space of vectors s representing the set of effective degrees of freedom available to the system (i.e., those relevant to the observed output); an embedding vector x at any time is representative of (that is, can be mapped to through a smooth change of coordinates to) the effective internal degrees of freedom s. A time-varying input stimulus u(t) applied to the system produces a time-varying output state y, which may be a scalar (e.g., temperature or amplitude) or vector quantity. The present approach allows for characterization of system behavior in terms of a plurality of local models, each of which is valid over a range of relevant degrees of freedom according to a probability distribution, so that each applies with greatest accuracy within a probabilistically defined subset of the degrees of freedom.

This is shown in FIG. 1. Suppose that the true physical behavior of the system (expressed in terms of the embedding vector) is nonlinear and somewhat stochastic, as described by the indistinct parabola 120; in other words, various output states of the system are associated, in a nonlinear fashion, with different values of x that are determined or affected by an input stimulus. The true behavior 120 can be approximated by a pair of local linear models 122, 124, each of which is valid over--and in accordance with--a respective probability distribution 126, 128 associated therewith and defining a domain of influence. Distributions 126, 128 can also be viewed as defining an "input width," i.e., a probabilistic expression of the relevant input and internal degrees of freedom. The "output" (y value) produced by the local model can similarly be expressed probabilistically (as indicated, in the case of local model 122, by an output probability distribution 130), so that the output value has a stochastic width rather than a precise value (concording with the true behavior 120).

In accordance with the invention, the joint probability density p(y,x) is expanded into a series of clusters c.sub.m, each of which is associated with a local model, a probabilistically defined domain of influence, and an output distribution, all as defined above. The domain of influence can be defined in the input space, and can optionally be defined in time, e.g., according to the theory of recursive estimation--that is, based on a history of inputs and outputs, and cluster-transitioii probabilities.

Thus, while the complete joint density is ultimately used to predict values of y given an input stimulus and its expected effect on x, the joint density is represented by a plurality of local models, thereby facilitating the use of simple (or at least locally accurate) models whose respective contributions to the joint density do not extend beyond their true relevance to system behavior.

The joint density is separated into a conditional probability given a cluster, multiplied by a cluster weight: ##EQU4## The conditional probability may then be separated into input and output parts as follows: ##EQU5## for M clusters (and, consequently, M local models), each local model being represented by a different instance of the function .function.(x,.beta..sub.m) with unknown parameters .beta..sub.m. The first factor, corresponding to p(y.vertline.x,c.sub.m), represents the output probability distribution, while the second factor, corresponding to p(x.vertline.c.sub.m), corresponds to the input probability distribution over the local model's domain of influence. The parameter p(c.sub.m) represents the weighting associated with the cluster. The quantity .sigma..sub.m,y corresponds to the width of the output distribution of the local model m, while .sigma..sub.m,d denotes the width of the input distribution. The quantity .mu..sub.m,d specifies the mean of the input distribution, i.e., the point corresponding to maximum probability (and maximum accuracy of the model m). The vector x is of the form x=(x.sub.1, . . . ,x.sub.d, . . . ,x.sub.D), where D is the embedding dimension; for example, in modeling a musical instrument, each of the x terms may be a time-lagged value of the audio output and the sensor readings that represent a player's manipulations of the instrument.

The foregoing framework assumes a scalar output and separable Gaussians; it can be straightforwardly generalized to vector output (i.e., where y has more than one component), and using non-separable input Gaussians with covariances. Specifically, the input term ##EQU6## uses separable Gaussians and just the variances; non-separable Gaussians can be represented by replacing this term with ##EQU7## where [C.sub.m ].sub.ij ={(x.sub.i -.mu..sub.m,i)(x.sub.j -.mu..sub.m,j }.sub.m is the cluster-weighted covariance matrix. The first input term (Eq. 2c) requires D parameters (the variances), while Eq. 2d requires D.sup.2 parameters (the covariances). The former is better suited to high dimensions, but forces the clusters to line up with the axes. The latter is more flexible, but requires excessive storage in high dimensions.

Similarly, the output term ##EQU8## can be generalized to ##EQU9## where [C.sub.m ].sub.ij ={(y.sub.i -.function..sub.i (x, .beta..sub.m))(y.sub.j -.function..sub.j (x, .beta..sub.m))}.sub.m.

In practice, x is generally determined by embedding. In particular, time-delay embedding of a set of N known input stimuli u.sub.1 . . . u.sub.n and the corresponding system outputs y.sub.1 . . . y.sub.N (i.e., a set of experimental or "training" measurements {y.sub.n,u.sub.n }.sub.n=1.sup.N) characterizes y in terms of the relevant system degrees of freedom comprehended in x (i.e., as a set of corresponding parameters of the form {y.sub.n,x.sub.n }.sub.n=1.sup.N). The greater the embedding dimension, the more accurate the system characterization will be; an excessively large dimension, however, wastes resources and may result in some inaccuracy (as discussed further below).

To use the approach of the invention, the user need specify at the outset only the number M of clusters and a form of the local model f reasonable for the system under study. M in effect varies the number of system "features" to be explained and is used to control underfitting and overfitting; that is, if the model captures too many output features, it will be "overfitted" and inappropriately include noise in its representation of system behavior. A technique for optimizing M is discussed below.

It is possible to use clusters with different types of local models if more than one type of local relationship is expected. The unknown model parameters (i.e., all terms of Eq. 2b having subscripts that include m) and the optimal embedding dimension D are both determined by separate iterative procedures, also as discussed below.

The output distribution p(y,x.vertline.c.sub.m) is given in terms of the likelikhood of the output data given the model; it is necessary to invert this probability in order to maximize the likelihood of the model given the data: ##EQU10##

The posterior probabilities can be used to find the cluster weights as follows: ##EQU11##

Eq. 4 follows because an integral over a density can be approximated by an average over variables drawn from the density. The practical result is the ability to utilize a sum of experimentally derived measurements rather than an integral over a continuous function.

As stated previously, the approach of the present invention does not require preliminary specification of model parameters. Using the foregoing framework, the model parameters may be found from the data itself using an iterative technique such as the Expectation-Maximization ("E-M") procedure, alternately calculating the likelihood of the data p(y,x.vertline.c.sub.m) given the current parameters, then choosing new parameters that maximize the posterior probability p(c.sub.m .vertline.y,x) given known data. Initially, a trial set of cluster parameters is chosen at random, and the process continued until the overall likelihood reaches a maximum.

More specifically, following initial calculation of p(c.sub.m .vertline.y,x), the expected input mean .mu. of each cluster is computed as follows: ##EQU12## where the last quantity defines the cluster-weighted expectation of the mean (that is, with reference to FIG. 1, the mean value of x within either input probability distribution 126, 128 associated with a respective local model 122, 124). The input variances (again with reference to FIG. 1, the variances of the probability distributions 126, 128) are computed from the updated mean, current values being iteratively updated by taking the cluster-weighted expectation of the input distribution

Starting with randomly chosen cluster parameters, the iterated parameters rapidly converge on the local maximum of the probability of the data set given the model, and through the sum over clusters in the denominator of the posterior probability, the clusters interact so that they each specialize in data poorly explained by other clusters.

Once the cluster input parameters are found from these expectations, the model parameters .beta..sub.m are found by maximizing the likelihood of the data as a function of the cluster parameters: ##EQU13## Since the logarithm function is monotonic, it is possible to maximize the expected value of the logarithm instead: ##EQU14##

The .beta..sub.m values are used to find the output variances:

The object of the analysis, of course, is to predict y from an instance of the feature vector x (that is, the inaccessible relevant system degrees of freedom represented by accessible time lags on y and u)--for example, to obtain the expected value of y given a particular feature vector x, or {y.vertline.x}. This quantity is given as ##EQU15##

Thus, with the model parameters and cluster weights characterized, the expected value of y is straightforwardly obtained from a given feature vector x. Optionally, it is possible to store the clusters hierarchically, such that it is unnecessary to include contributions from distant clusters.

The total variance in y--that is, the total output error--is the weighted sum of the variances .sigma..sub.m,y.sup.2 associated with each local model. It is also possible, however, to characterize the expected error for y given a particular feature vector x, or {.sigma..sub.y.sup.2 .vertline.x}. This quantity is given as ##EQU16##

For a stochastic system, it is also possible to generate samples drawn directly from the distribution to emulate the behavior of the modeled system.

If the output model is written as a linear sum of basis functions ##EQU17## Defining a.sub.j ={y.function..sub.j (x)}.sub.m and B.sub.ji ={.function..sub.j (x).function..sub.i (x)}.sub.m, the new value of .beta. can be found from .beta..sub.m =B.sup.-1 .multidot.a.

For example, the local models 122, 124 in FIG. 1 are linear functions in .beta.. The model parameters B.sup.-1 and a are calculated by iterative analysis of the joint probability distribution according to Eqs. 2a-2g. The inverse of the matrix B is preferably a pseudo-inverse.

The foregoing discussion describes use of the invention to model the physical behavior of a driven, dynamic system based on a training set of inputs and outputs. The analysis assumes, however, an optimal number of embedding dimensions as well as M, the number of local models (and clusters). The variance quantities discussed above can be used to identify optimal values for these parameters. Refer to FIG. 2, which illustrates the manner in which varying the embedding dimension (given the training set {y.sub.n,u.sub.n }.sub.n=1.sup.N) affects the total error {.sigma..sub.y.sup.2 }. The total error stops decreasing--that is, the curve reaches a plateau and no further improvements in error are obtained--when the embedding dimension d is sufficiently high. Values of d above this optimal level are superfluous, wasting computational resources and possibly introducing some unwanted modeling of noise.

For purposes of determining optimal numbers of clusters, one can withhold a subset of the measured data for testing the model; then by the technique of cross-validation, the model is developed using the training data, and its performance evaluated against both the training and testing data. FIG. 3 illustrates the relationship between M and the total variance .sigma..sub.y,train.sup.2 within the training set {y.sub.n, u.sub.n }.sub.n=1.sup.N (i.e., the data not withheld for testing), and the total variance .sigma..sub.y,test.sup.2 within a set of test performances of the model (that is, predicted values of y from sensed input stimuli u). Although increasing values of M naturally reduces training error on a continuous basis, excessive levels of M actually increase testing error. As noted earlier, this results from "overfitting," where the model is excessively precise and captures noise as well as legitimately generalizable behavior. On the other hand, too few clusters "underfit" and fail to provide an adequately complete behavior model. Accordingly, goodness of fit can be modeled as an error function, and an optimal value of M chosen by minimizing testing error. If the testing error cannot be brought sufficiently low or if the behavior shown in FIG. 3 is not observed, the number of sensors used to record input stimuli may be insufficient.

2. Operation of the Invention

As shown in FIG. 4, operation of the invention can be expressed in a series of steps, with iterative loops among steps followed as necessary for optimization. The illustrated steps assume, for purposes of illustration, modeling of a musical instrument, it being understood that the invention is not limited in its capability to this domain of complex, driven physical systems. In a first step 400, a series of sensors for measuring input and output is arranged about the instrument; for example, the output sensors might be electronic transducers that convert sound or vibration into an electronic signal, while the input sensors might detect the instantaneous position of a user's hand or finger with respect to the instrument. A sensor arrangement suitable for, e.g., bowed musical instruments is disclosed in U.S. Pat. No. 5,247,261, the entire disclosure of which is hereby incorporated by reference.

In step 405, the user plays the instrument with the sensors active and providing signals to, for example, a digital computer (by means of suitable analog-to-digital converters, if necessary) programmed to accumulate and assign values to these signals. In step 410, the data are arranged in timewise sets, i.e., segregated so that data obtained at identical times remain grouped. These training data form the basis of a time-lag embedding step, the resulting lag vector including time lags on sensed values for both input and output.

Based on the lag vector and a selected number of local models of system behavior (whose form, but not parameter values, is known), the behavior of the system is modeled in accordance with the cluster-weighted approach set forth above. In particular, the joint probability density is characterized from the clusters, enabling prediction of future output values based on current and lagged states and a given input stimulus (i.e., {y.vertline.x}). In step 425, the model is employed to synthesize the original physical system. For example, the user might "play" a noiseless replica of the original instrument, the input sensors generating predicted output (e.g., driving loudspeakers) based on the model.

Steps may be repeated in accordance with any of four processing loops 450, 455, 460, 465, iteration along which is used to optimize various aspects of performance. Step 415, for example, requires multiple time-lagged measurements in accordance with the embedding dimension D, necessitating repeated collections of data to form the training set. The first execution of the embedding step 415 uses a small embedding dimension D, which is increased until the approximation is adequate; however, as expected output variance (i.e., {.sigma..sub.y.sup.2 }) is characterized through implementation of the model, D can be adjusted (loop 455) to an optimal value. Similarly, modeling step 420 can reveal inadequacies in the number and/or configuration of sensors if the testing error cannot be brought sufficiently low or if the behavior shown in FIG. 3 is not observed; in this case, step 400 is repeated (i.e., the sensor array is altered and/or augmented) until the testing error observed as a result of step 420 is abated (loop 460).

Finally, actual use of the model accomplishes two important objectives. First, a testing set of inputs and outputs is used to establish the optimal number of clusters and local models (loop 465) in accordance with observed output errors; and also to ensure that the set of training inputs is adequate to generate a model robust enough for its intended use. For example, if the training set is too limited in terms of exploring the system's dynamics, the model will not include relevant internal degrees of freedom, and will not produce accurate output through the full range of possible inputs. In this sense, however, "too limited" is meaningful only in terms of intended use; so long as the training set is consistent with the range of possible inputs, the model will produce accurate output.

3. Application of the Invention

Refer now to FIG. 5, which illustrates, in block-diagram form, a representative hardware embodiment of the present invention. The apparatus, indicated generally at 500, is driven by a central-processing unit ("CPU") 510, which is typically a high-speed microprocessor or microprocessor array capable of processing, in real-time, the model derived as described above to generate predicted output values as the user provides input. Apparatus 500 includes a main computer memory 515, which contains a group of modules that control the operation of CPU 510 and its interaction with the other hardware components. An operating system (not shown) directs the execution of low-level, basic system functions such as memory allocation, file management and operation of mass storage devices. At a higher level, an analysis module 520, implemented as a series of stored instructions, directs execution of the primary functions performed by the invention, as discussed above. In particular, analysis module 520 executes steps 410, 415, 420 shown in FIG. 4. Associated with analysis module 520 is a database 525, which is generally a memory partition, and which accumulates data associated with training, testing, and use of the model generated by analysis module 520, as well as data representing the parameters of the model itself.

A set of output sensors 530 (as in, for example, the '261 patent) gather data representative of the output of the system under study in response to physical manipulations, while a series of input sensors 540 (e.g., in the context of a musical instrument, a similar set of sensors arrayed on a replica of the instrument, or a keyboard each of whose keys represent a predetermined input mode) gather data representative of those manipulations. The outputs of all sensors are provided as digital signals to CPU 510, which stores them in database 525. As indicated previously, data corresponding to the training set is grouped in a timewise fashion to facilitate embedding synthesis.

Following generation of the model, the user provides input to the apparatus 500, for example by "playing" a replica of the system under study so that the user's gestures, movements or other operations relevant as system input are monitored and provided to the model, which computes an expected output (e.g., an audio time series) therefrom. It should be stressed that the arrangement used to "play" (i.e., provide input to) the model is a matter of design choice, and may be quite different from the sensor arrangement initially used to gather data to build the model. All that is necessary is an arrangement capable of capturing or otherwise representing manipulations representative of (or which can be mapped to) the manner in which the original system is manipulated.

That output is presented in real-time by means of an output stage 550, which ideally is capable of producing responses characteristic of the physical system itself. In the case of a musical instrument, for example, output stage 550 includes suitable digital-to-analog, amplifier and audio-processing electronics, as well as a loudspeaker, for producing an audio response.

Operation of the invention, and the dimensional reduction obtainable therewith, are depicted in FIGS. 6A, 6B. The physical system illustrated by the figures is a sticky bow drawn across a damped string. With reference to FIG. 6A, a numerical model is used to represent the action of the bow along the string; the slowly descending, discretely sampled curve 610 shows the changing position of the bow as it is drawn, and the curve 620 is the resulting output audio time series. (In other words, the y-axis with respect to curve 610 represents position, while with respect to curve 620 represents amplitude.) The highly accurate but computationally intensive numerical model used to generate FIG. 6A requires 54 degrees of freedom.

Using the techniques set forth above, a predictive model was generated using an embedding vector consisting of five lags of audio and two samples of bow position, for a total of seven degrees of freedom. An audio time series 630 was generated from the predictive model by specifying the illustrated bowing action. As can be seen from comparison of the two figures, the predictive model accurately represents the far more complex "true" behavior of the numerical model with far fewer degrees of freedom.

It should again be emphasized that the example of a musical instrument is for purposes of presentation and not limitation. The present invention is capable of modeling the behavior of a wide range of driven, dynamic systems. The invention may be used, for example, to model the flight behavior of an aircraft in response to various manipulations of aircraft controls (stabilizer, rudder, flaps, thrust, etc.) or environmental inputs (e.g., to model the effects of turbulence), the model being implemented as a flight simulator. In another implementation, the invention is used to model non-linear circuit elements (e.g., microwave devices the outputs of which depend nonlinearly on the input); the device is placed in a circuit, and the above-described techniques used to define a model of device behavior.

It will therefore be seen that the foregoing represents a computationally tractable and highly flexible approach to synthesis of complex physical systems. The terms and expressions employed herein are used as terms of description and not of limitation, and there is no intention, in the use of such terms and expressions, of excluding any equivalents of the features shown and described or portions thereof, but it is recognized that various modifications are possible within the scope of the invention claimed.