Frequency Domain Kernel Estimation for 2nd-order Volterra Models Using Random Multi-tone Excitation

We consider the problem of frequency domain kernel estimation using random multi-tone (harmonic) excitation for 2nd-order Volterra models. The basic approach is based on least squares minimization of model output error, and results for the Volterra kernel estimations with random multi-tone inputs and random Gaussian input are compared. We show that kernel estimation with multi-tones are very accurate and efficient compared to the latter. As an illustration, the proposed method is applied to a discrete input–output system obtained from the numerical simulation of a representative hydrodynamic system for modeling semiconductor device transport. We also consider the effect of noise in the kernel estimation.


INTRODUCTION
The Volterra (functional) series approach has been widely used to model the forced response of weakly nonlinear systems [1,2]. In practice, many nonlinear systems can be modeled approximately using only the first few terms of the Volterra series. Identification of such a finite Volterra model requires measurement or estimation of the associated kernel functions from a finite set of representative input -output data (realizations) of the actual system. It is important that the identification be carried out with fewer realizations (to reduce the cost) while yielding maximum information about the kernel functions with a good accuracy. The objective of the present study is to develop an efficient and practical approach, in this sense, for identification of kernel functions of a 2nd-order Volterra model in the frequency domain. The approach is applicable to the analysis of physical systems or numerical simulations, and the focus of the present work is its applications to the analysis of nonlinear mathematical models of physical systems using numerical simulations. In particular, we show that, precise and controlled random multi-tones can be used effectively to obtain realizations.
There is an extensive literature of studies proposing various identification methods using different kinds of inputs for certain systems. In particular, methods for identification of Volterra kernel functions in the frequency domain † (i.e. Volterra transfer functions, VTFs) can be divided into two main groups based on the use of tonal excitation or random excitation. In the latter approach, methods to estimate the linear, quadratic and cubic transfer functions for orthogonal and non-orthogonal models with both non-white Gaussian and non-Gaussian random inputs can be found in the works of Hong et al. [3 -6]. In these studies, kernel estimation with non-Gaussian excitation is based on the least-square minimization (LSM) of the model error. Some drawbacks of using random inputs for estimating VTFs are due to aliasing. First, recall that, it is impossible to generate random input signals with perfect cut-off in its power spectrum. This certainly affects the accuracy of the estimation using random input. Moreover, the number of realizations required for an accurate estimation may be quite large depending on system characteristics. This significantly increases the expense of obtaining realizations. This subject is of general interest in many application areas of nonlinear analysis, including fluid and structural dynamics or circuit and device modeling. For example, in the nonlinear identification of structural systems, the disadvantages of using strictly random or burst-random Inputs have been reported [7]. In practice, the steady-state sine, fast sine, or burst sine are found to be more optimal. The preliminary motivation for our work is to extend LSM techniques for estimation of VTFs with random multi-tone inputs, which eliminates aliasing effects. With this in mind, we, in fact, aim to obtain an estimation method using random multi-tones, which has certain advantages over the conventional techniques to measure the VTFs.
The actual frequency domain formulation of a Volterra model gives the amplitude and phase relationships between sinusoidal tones contained in the inputs and outputs of a system in terms of possible tonal interactions along with the associated VTFs [1,8 -10]. Efficient methods for measuring VTFs by using multi-tone sinusoidal inputs have been reported by several researchers. In these studies, Boyd et al. [11] presented a "quick method" for calculation of the transfer functions of a second-order Volterra system in the frequency domain. Chua and Liao [12] extended the method to third order models, which includes a cubic term. Each multitone input is constructed from the tones at carefully chosen frequencies so that the outputs of different order Volterra operators contain the tones at distinct interaction frequencies, at which the VTF values can be measured explicitly. Thus, the number of frequency coordinates of measured VTFs is maximized for a finite number of multi-tone excitations. On the other hand, with this approach, explicit computation of transfer function values is possible for only some scattered isolated frequency points, not for all discrete spectral points in a given frequency range. (For example, any quadratic Volterra kernel value associated with the DC contribution of the tones cannot be computed explicitly with one multi-tone input.) Furthermore, for multi-tone inputs containing tones at all discrete frequency values within a frequency range, it is not possible to measure the VTFs explicitly since the tonal interactions overlap. The discrete VTF values can be obtained only by solving linear algebraic systems of equations (constructed from the input -output realizations). In principle, the number of realizations can be taken to be equal to the number of unknowns in the algebraic systems. However, in practice, it is desirable to solve over-determined systems (with more realizations), as done in estimations with random inputs. In this way, the effects of imperfections in the Volterra model can be minimized. This is in essence what is done in kernel estimations with random inputs. Our estimation procedure follows this approach and uses multi-tone inputs whose tones are chosen at each consecutive discrete frequency with uniformly distributed random phases and Gaussian distributed random amplitudes. The utilization of random amplitude and phase represents an essential difference from previous multi-tone methods, which, for the most part, used deterministic tones. Sampling the steady-state harmonic outputs to these multi-tone excitations within a record length, which exactly equals one period, prevents aliasing and it also enables a sharp cut-off in the input power spectrum and estimation at each discrete frequency value of the domain. Applying LSM to estimate the VTFs using random realizations makes it possible to check the goodness of the model by using conventional tools such as the power spectrum and coherency. Because of the parallelism in the present approach and conventional techniques using random inputs in kernel identification, we use the term "estimation" rather than "measurement".
The outline of the paper is as follows: We first briefly describe the LSM approach and the harmonic excitation strategy for obtaining the necessary time series of the system inputs and outputs. Then, we compare the frequency domain estimations of the linear and quadratic kernels of the 2nd-order Volterra model using the LSM principle with Gaussian random and random multi-tone inputs. As an illustrative test case, we consider a 2ndorder Volterra model for voltage-current dynamics of a semiconductor device, simulated by a numerical hydrodynamic code. In particular we explore the effects of noise in the system input and output signals, which would be of concern mainly in physical experimental applications as opposed to numerical experimental applications. (However, the effects of approximation error in numerical input can also be interpreted as "noise" and the influence of this on the output and kernel estimation is relevant.) The output of the approximate 2nd-order Volterra model obtained using the present method is compared with the actual output of the numerical model. Finally, we summarize and discuss our findings in the conclusion.

2ND-ORDER VOLTERRA MODEL IN THE FREQUENCY DOMAIN
The 2nd-order Volterra model v 2 consists of linear and quadratic operators acting on an input function x (t ). That is, where l(t ) and q(t,t 0 ) are the associated linear and quadratic kernel functions in the time domain. In the frequency domain, Eq. (1) transforms to Qðg; f 2 gÞXðgÞXð f 2 gÞ dg ð2Þ where X( f ), V 2 ðXÞð f Þ; L( f ) and Q( f, f 0 ) denote the Fourier transform of x(t ), v 2 ðxÞðtÞ; l(t ) and qðt; t 0 Þ; respectively. Here, the quadratic kernel may be assumed to be symmetric without loss of generality.
If the 2nd-order Volterra model in the frequency domain is written more compactly in terms of the respective linear and quadratic operators in Eq. (2) as then the power spectrum of V 2 (X ) can be expressed as where k·l denotes the expectation operator implying the average value of the spectral quantity applied to a large number of signal realizations. R denotes the real part of the associated expression and symbol * implies the complex conjugate. The first two terms in the right-hand side of Eq. (4) represent the individual contributions to output power from the linear and quadratic operators. The last term is the interference term, and vanishes when L(X ) and Q(X ) are orthogonal. The orthogonality condition is satisfied if the bispectrum kX* ð f ÞX* ð f 0 ÞXð f þ f 0 Þl of the input x(t ) is zero. This corresponds to the fact that the sampling frequency to collect the time series for input and output should be chosen such that, the Nyquist frequency f s /2 is smaller than 2f M to avoid aliasing in the output, i.e. f s $ 4f M : When f M ¼ f s =4; for a given N, the largest spectral range is attained with M ¼ N=4 for the identification domain of the quadratic kernel. Now, assume that a single input -single output (SISO) relationship between an input x and an output y of a nonlinear system can be approximately characterized by a 2nd-order Volterra model in the frequency domain. That is, , corresponding to a set of sampled system input -outputs (realizations). Then, the Volterra kernels given in Eq. (5) can be calculated approximately from the LSM of the power spectrum for the model error. That is we seek the discrete values of the Volterra kernels in Eq. (5) such that the least squares functional is minimized for each k. Here, k·l denotes the expectation operator, which is to be approximated by an ensemble average over N r realizations. Minimizing Eq. (7) requires for each k and m in the identification region. This leads to an algebraic system of the following block-matrix form for each k: where the scalars S XX [k ] and S XY [k ] are the values of discrete auto and cross power spectra defined by The entries in the vectors q (k ) , b (k ) and y (k ) and in the matrix C (k ) are, respectively, Remarks In practice, the symmetry of the quadratic kernel, i.e. Q½m; m 0 ¼ Q½m 0 ; m; is taken into account. This implies we need to consider the kernel values only for m $ m 0 ¼ k 2 m (on the solid part of the line in Fig. 1), so the dimension of the algebraic system is reduced approximately by half.
For Gaussian inputs, the orthogonality condition implies b ðkÞ ¼ 0 so Eq. (9) is decoupled and L[k ] and q (k ) can be obtained independently. Furthermore, it can be shown that the 4th-order spectral moment c i 0 j 0 vanishes for i 0 -j 0 ; making C (k ) diagonal. This permits the values of q (k ) to be calculated explicitly [3]. However, for a finite number of realizations, the conditions c i 0 j 0 ¼ 0 for i 0 -j 0 and b ðkÞ ¼ 0 only approximately hold, so calculations that do not make these assumptions are expected to produce more accurate results for the transfer function estimates [13].

Coherency
Coherency is defined as the ratio of the model output power to actual output power.
It can be shown that the coherency function, ‡ calculated using the estimated kernels, is always bounded by zero and unity. Based on the model output as given in Eq. (3), the coherency function may be also written as a sum of partial coherency functions, i.e. g 2 ¼ g 2 l þ g 2 q þ g lq . When orthogonality holds the interference term g lq vanishes, and the individual linear and quadratic coherency terms are less than unity. Partial coherencies are useful in understanding the "strength" of the linear and quadratic mechanisms in terms of their contribution to the total power as a function of frequency [3]. In this study, we are concerned with the total coherency g 2 and the coherency of the first order model g 2 l for comparing the 1st and 2nd-order Volterra model outputs.

Identification With Multi-tone Excitation
In harmonic excitation, multi-tone inputs may be prescribed in the time domain in the form where M is the total number of (tones) harmonics. The selection of amplitudes {A m } determines the class of input signals for the identification. In particular, here, we consider the amplitudes {A m }, and the phases {f m }, taken to be random numbers with uniform distributions in the ranges ½2A; A and ½0; 2p; respectively, where A is the reference input amplitude. Here, we stress the fact that, in conventional methods using multi-tone inputs, constructing the tones with deterministic amplitudes and phases can result in some useful simplifications in measuring the Volterra kernels explicitly. In the present method, however, choosing the amplitudes and phases randomly as described above actually makes the matrix C (k ) in Eq. (9) much more diagonally dominant, ‡ Note that, with this definition, the coherency function reduces to the well known form g 2 ¼ jS XY j 2 =S XX S YY for linear systems.
and therefore, yields more robust computations for the kernel estimations, as confirmed in our numerical experiments.
In kernel identification using random multi-tone inputs, each distinct system realization corresponds to a pair of a multi-tone inputs with a different set of {A m } and {f m } and the steady-state harmonic response of the system to this input. For accurate identification (or estimation) of Volterra kernels, it is essential that the output data from the simulations be collected after the transient dynamics of the system die out.
The complete kernel identification procedure is outlined below:  9) and solve for the L and Q values at the discrete points in the associated range. 7. Having obtained the discrete kernels L and Q in the entire identification domain, calculate the model output from Eq. (5) for the same inputs. 8. Compute the coherency from Eq. (15) to check the "goodness" of the model.

COMPARISONS OF ESTIMATIONS WITH RANDOM INPUTS AND RANDOM MULTI-TONE INPUTS
In this section, we are interested in the performance of the LSM procedure with random and random multi-tone inputs. For simplicity and to isolate the effects due to higher order nonlinearities, we consider the following pure quadratic discrete SISO system which is stable when jaj , 1: This difference equation has a delay effect on the input -output relationship between x and y. If x[n ] and y[n ] are assumed to be sampled with sampling frequency f s ¼ 1; in the frequency domain the input -output system given by Eq. (17) can be easily shown to be represented by a quadratic Volterra operator, whose kernel function is known analytically as Thus, the accuracy of the estimations made with Gaussian random and multi-tone input functions can be compared directly. In our "numerical" experiments, Eq. (17) was sampled with N ¼ 64 points for each realization.
In the random excitation case, the entire set of realizations was obtained from the solution to the difference equation in Eq. (17) with the numerically generated Gaussian random input x[n ]. The input power spectrum is band-limited to f =f s ¼ 0:25: Figure 2 shows the power spectra of the input and output signals obtained for the discrete system Eq. (17) with a ¼ 0:85: Here, the number of realizations was taken to be 100 where each system realization corresponds to the non-overlapping data records (of N ¼ 64 points) for the entire time series of x[n ] and y[n ].
In the random multi-tone (harmonic) excitation case, 16 multi-tones with randomized amplitudes and phases were used in the input signal ðM ¼ 16Þ with the same number of sampled points ðN ¼ 64Þ: Figure 3 gives the power spectra of the input and output signals obtained with the number of realizations N r ¼ 100 for the same discrete system ða ¼ 0:85Þ: Note that the case of M ¼ 16 The magnitude of the quadratic kernels estimated from Gaussian and random multi-tone inputs with N r ¼ 100 realization are compared in Fig. 4. With multi-tone input, the exact quadratic kernel is recovered within the machine accuracy, being independent of number of realizations used for the estimation. The error in the quadratic kernel estimated using random input can be reduced by increasing the number of realizations; however, the convergence to the exact kernel is rather slow. Figure 5 shows the errors in the quadratic kernels with estimated random input of 100 and 3200 realizations. In this figure, it is clear that, even though, the error in the estimation generally reduces with more realizations, the decrease in the maximum local error is small for a ¼ 0:85 with even large N r . One should notice that, the maximum error in the estimated kernel occurs on the frequency coordinate ð0:25; 20:25Þ: This implies that the lack of a sharp cut-off in the power spectrum of the input near f ¼ 0:25 causes this error in the estimation. This situation is more apparent when the discrete system is influenced more and more by the time delay as the parameter a approaches unity. In Fig.  6, we compare the normalized L 2 norms of the errors in the quadratic kernels estimated using Gaussian random input for different number of realizations with a ¼ 0:75; 0.85 and 0.95. As can be seen in this figure, it becomes very difficult to obtain accurate estimations as the parameter a approaches unity using a practical number of realizations with Gaussian random inputs. (See also the test studies presented in Ref. [13]. On the other hand, using random multi-tone inputs recovers the exact kernels within machine accuracy (1.0E 2 6) for any jaj , 1: As expected, additive input noise in the multi-tone realizations result in inaccuracies in the kernel estimation. By using more (multi-tone) realizations, one can anticipate decreasing the error in the kernel estimation. This aspect will be discussed in more detail in "Noise effects in kernel estimation".

APPLICATION TO SEMICONDUCTOR DEVICE MODELING
In this section, we demonstrate an application of the proposed random multi-tone input method for the identification of a 2nd-order Volterra model in the frequency domain to approximately represent the transient FIGURE 4 Quadratic kernels esimated using Gaussian random (left) and multi-tone (right) inputs with 100 realizations. voltage-current dynamics of a simple n þ 2 n 2 2 n þ diode about an operating (steady state voltage-current) point. The identification is carried out using computer simulations of the system. There are several models and numerical techniques in use and under investigation for the steady-state device problem [14 -16]. Our interest is in the nonlinear system identification for the transient behavior. We believe that this topic has not previously been investigated. Here we consider a representative model for transient simulation of a very simple but standard test device [17]. We emphasize that the focus is the nonlinear identification problem so the numerical scheme and test problem is of less importance.
The computer simulation of the device dynamics is based on the numerical solution of the one-dimensional hydrodynamic transport equations [18] using the Mac-Cormack scheme [19]. In these hydrodynamic transport equations, the heat transfer effect is not included while the full convection model for drift kinetic energy is adopted [20] (the approach can be applied to other device models [14]). In the case study, we consider a 0.3 micron device with lattice temperature 290 K. The permitivity of the device is taken to be e ¼ 1:0536 £ 10 210 C 2 /J m. The electron concentration at the gates of the device is 5 £ 10 5 /m 3 while the doping is 2 £ 10 3 =m 3 : Here, a moderate doping profile has been chosen in order to decrease the time-step in the simulations and to ensure weak nonlinearity in the system (more demanding device scales and doping could be similarly studied but we have not elected to do this because out objective is specifically to demonstrate the identification problem). The steadystate solution to the device equations obtained for the nominal voltage V o ¼ 1:6 V gives the corresponding nominal current value I o . Then, a time-dependent voltage perturbation vðtÞ ¼ VðtÞ 2 V o is applied as an input function. The current perturbation from the nominal current, i.e. iðtÞ ¼ IðtÞ 2 I o ; is computed from the numerical solution of the hydrodynamic transport equations, and taken as a response to this voltage perturbation. Thus, we define our "system" as the SISO relationship between temporal perturbations from the nominal voltage v(t ) and current i(t ), as observed in the numerical simulations, and obtain a 2nd-order Volterra model approximating this SISO system. Since the characteristic time scale of the system is in the order of 10 212 s, we present and discuss the results of this analysis in terms of scaled frequency ( f 0 ) and time (t 0 ) which are related to the physical frequency and time by f ¼ f 0 £ 10 12 Hz and t ¼ t 0 £ 10 212 s; respectively. The kernels of the 2nd-order Volterra model were estimated using 200 input -output realizations obtained from the numerical solutions of the hydrodynamic transport equations. The input signals included 16 tones, and the frequency of the maximum tone was 2.5 £ 10 12 Hz. The reference input amplitude for the randomly selected multi-tones is A ¼ 0:05 V: The power spectra of the multi-tone inputs and outputs (averaged over the 200 realizations) used for the estimation are given in Fig. 7 in the scaled frequency range f 0 [ ½0; 5: The resulting estimated linear and quadratic kernels are shown in Figs. 8 and 9, respectively. The high coherency of the model output as seen in Fig. 10 indicates the Volterra model estimated at this amplitude level represents the dynamics of the numerically simulated device quite closely. This is mainly because the estimation has been carried out at small amplitude levels and the input does not strongly excite the higher order nonlinearities.
The estimated 1st-and 2nd-order Volterra models (based on A ¼ 0:05 V) were tested with new sets of higher amplitude random multi-tone inputs. The 1st-order model includes only the linear operator associated with the estimated linear transfer function while the 2nd-order model also includes both the linear and quadratic operator associated with the estimated quadratic transfer function. The frequency range of the test inputs is half that used for the kernel estimation; that is the number of tones in the test inputs is M ¼ 8: Comparisons of power spectra and coherency were made for the outputs in response to the inputs whose reference amplitudes are A ¼ 0:1 and 0.3 V. Note that these amplitudes are two and six times as large as the reference amplitude ðA ¼ 0:05 VÞ used to determine FIGURE 7 Power spectra of the input and output semiconductor device. ð f ¼ f 0 £ 10 12 HzÞ:   the 2nd-order Volterra model. Figure 11 shows the power spectra of the outputs of the actual system (which is based on the numerical solutions of hydrodynamic equations) and the estimated Volterra model for these two amplitude levels of the input. In this figure, we see that the cut-off frequency of the input signals is f 0 c ¼ 1:25 while the system outputs have power in the entire frequency range. The 1st-order model gives outputs in the frequency range f 0 [ ½0; 1:25 whereas the 2nd-order model gives outputs in the frequency range f 0 [ ½0; 2:5 since it contains the contribution of the quadratic operator, too. As one can see in this figure, the system output power in the range f 0 [ ½2:5; 1Þ is due to the cubic and higher order nonlinearities in the system, which are not included in the 2nd-order model. As the amplitude level of the inputs are increased from A ¼ 0:1 to 0.3 V, the system output power in the range f 0 [ ½2:5; 5 increases significantly. Of course, the higher order nonlinearities also affect the system output in the range f 0 [ ½0; 2:5: Coherency functions of the 1st-and 2nd-order Volterra models are given in Fig. 12 for the two amplitude levels and indicate this effect. For A ¼ 0:1 V; the coherency function is close to unity for the two models in the associated frequency ranges. The 1st-order model is quite accurate in the range f 0 [ ð0; 1:25; whereas in the 2ndorder model, which contains the quadratic operator, the total system power is recovered in the entire range f 0 [ [0,2.5]. For, the higher amplitude case A ¼ 0:3 V; the coherency is degraded in ½0; 1:25; and with more significant departures in ½1:25; 2:5: This shows that, at this high amplitude, the total system power cannot be recovered by the linear and quadratic operators in these truncated models. One should note that the estimated 2ndorder Volterra model gives the best coherency (as shown in Fig. 10) for inputs at the same level of amplitude as the inputs that are used in the estimations; that is (A ¼ 0.05 V). The Volterra models are expected to become less accurate at the higher amplitude levels of inputs, and results in Fig. 12 demonstrates this.
The comparisons of the system output and model output in time, were made for an arbitrary input -output realization after the inverse DFT of the model output in the frequency domain were carried out. Figure 13 shows the actual and model outputs to a typical multi-tone input in the time domain. The reference amplitude of the input is A ¼ 0:3 V: The differences between the model outputs and actual output are plotted in Fig. 14 for the 1st-and 2nd-order Volterra models. The superiority of FIGURE 11 Power spectra of the input and outputs of the actual system and Volterrs models at the amplitude levels A ¼ 0:1 (left) and A ¼ 0:3 (right). the 2nd-order Volterra model to the 1st-order (linear) model in approximating the actual output locally is evident in this figure. As Fig. 11 indicates, the power of the actual output signal is low by a factor . 10 23 at high frequencies ð f 0 . 2:5Þ: Therefore, the output from the 2nd-order Volterra model approximates the actual output to a considerable accuracy in the time domain, even for this high amplitude input case. Although the 2nd-order Volterra model approximates the actual output waveform quite well, two singular exceptions occur at t 0 < 4:5 and 10.5 (see Fig. 14). Note that at these times the input amplitude is at its maxima, thereby possibly "exciting" higher-order nonlinearities. In such a case, the 2nd-order model cannot capture these higher-order effects.
Finally, in Fig. 15, the L 2 norms of the errors in 1st-and 2nd-order Volterra model outputs in the time domain are plotted for increasing reference amplitude of the excitation. When the reference input amplitude is increased, the effect of higher-order nonlinearities in the system, which are neglected in the truncated 2nd-order Volterra model, becomes stronger. Since the 2nd-order model takes into account the quadratic nonlinearity as well as the linear part in the system response, it gives a more accurate approximation to the system output compared to the 1st-order model as can be seen clearly in Fig. 15. However, for the higher amplitude levels, limitations of the 2nd-order Volterra model are also evident.

NOISE EFFECTS IN KERNEL ESTIMATION
Although noise is not a primary issue in system identification using numerically generated data, it may be present to a significant extent in experimental studies and affect the kernel estimation adversely. Here we consider a SISO system for which a 2nd-order Volterra model representation is sought by estimating the Volterra kernels using random multi-tone inputs. In this process, even though pure multi-tone inputs are to be applied to the system directly, the actual inputs on the system are expected to contain some noise due to imperfections in the exciter. Furthermore, the measured output signals also contain some noise. Figure 16 is a schematic for the kernel estimation process applied to a physical system with independent noises n 1 (t ) and n 2 (t ) added to the input and output signals, respectively. We now examine the (independent) effects of additive noise in the input and output signals on the kernel estimation of a 2nd-order Volterra system. Our 2nd-order test system to be identified is where the coefficient a controls the relative linearity. In the frequency domain, the discrete system given by Eq.  Table I gives the calculated errors in the estimated kernels of the 2nd-order Volterra series representation of Eq. (19) with the above a values for different N r and SNR when n 1 or n 2 is present in the system (without being applied simultaneously). Here, the errors are normalized as where E L and E Q are the (complex) differences between the estimated and the exact kernel functions, L and Q, and the domain of integration is the frequency domain in which the kernels are estimated. The four sections of the table correspond to the following cases, respectively: ðaÞ a ¼ 0:1; n 1 -0; n 2 ¼ 0; ðbÞ a ¼ 0:1; n 1 ¼ 0; n 2 -0; ðcÞ a ¼ 0:001; n 1 -0; n 2 ¼ 0; ðdÞ a ¼ 0:001; n 1 ¼ 0; n 2 -0: In Table I, we see that, for these four cases, as the SNR and the number of realizations increase, the errors in both kernels decrease regularly although the effect of increasing N r is quite weak. It is interesting to observe that the number of realizations needed for accurate estimation of kernels is quite high when the noise is present. This is obviously similar to the case of kernel estimation using random noise excitation described previously.
Here we can summarize our qualitative observations as follows: . The number of realizations needed for accurate estimation of kernels is quite high when the noise level is high. This is similar to the case of kernel estimation using random noise excitation described previously. . When the linear operator of the system is dominant, the linear kernel is estimated more accurately compared with the quadratic kernel. In this case, the accuracy of the estimated kernels are less for input noise relative to output noise at the same signal-to-noise-ratio (SNR) and N r . . When the system has a dominant quadratic nonlinearity, the estimation of the quadratic kernel is more accurate compared to the linear kernel. . Compared to output noise, input noise has a slightly more adverse effect on the accuracy of the estimated linear kernel. . For the same number of realizations, as compared to output noise, input noise decreases the accuracy in the quadratic kernel less when SNR is low, and more when SNR is high. . The accuracy of the quadratic kernel is not influenced by the relative strength of the linear operator when noise exists in the output.
Note that the above observations are valid when the identified system is actually a 2nd-order Volterra system.

FREQUENCY DOMAIN KERNEL ESTIMATION
If a system has higher-order nonlinearities, noise effects in the estimation of the kernels of a 2nd-order Volterra model will generally be different from those observed here, depending on the strength of higher order nonlinearities in the system input -output.
The methodology presented in this paper is well suited for numerical experiments. The fact that controlling inputs and observing output with high precision is essential in the present method raises questions about its applicability to physical experimental studies where noise can be an important issue. Naturally, it can be expected that the present method would give accurate and efficient estimations when the noise level is low. On the other hand, our study suggests that, when large amplitude noise is present in the input -output, considerably more realizations are needed in order to reduce the noise effect on the kernel estimation by ensemble averaging. In this case, the advantage of multi-tone inputs over random inputs diminishes drastically since both are expected to require a large number of realizations.

CONCLUSION AND DISCUSSION
In this work, we have presented an extension of the LSM approach using random multi-tone inputs for identification of the frequency domain kernels of a 2ndorder Volterra model. The present method employs a collection of independent multi-tones with random amplitudes and phases at the consecutive frequencies within a frequency range. Comparisons of the kernel estimations using random and random (amplitude and phase) multi-tone inputs for a simple test case (in the absence of noise) show that the latter gives considerably more accurate results. This is due to the fact the present method uses the steady-state responses to the multi-tone inputs and sampling the signals within their period.
The method has been successfully illustrated for a discrete SISO system for the voltage-current dynamics of a semi-conductor device (about an operating point) based on numerical solution of the hydrodynamic transport TABLE I Errors in the kernels of the 2nd-order Volterra model while noise is in input ðn 1 -0; n 2 ¼ 0Þ or noise in output ðn 1 ¼ 0; n 2 -0Þ for a ¼ 1 and a ¼ 0:001 cases (a) Noise in input signal (n 1 -0; Err. in Err. in L 0.97E þ 0 0.68E þ 0 0.81E-1 0.68E-2 Err. in Q 0.10E þ 0 0.10E-1 0.10E-2 0.11E-3 system of partial differential equations. It has given estimations for VTFs for regularly distributed frequency points on a discrete domain with a desired frequency resolution. Based on this experience, we found the present method very efficient, and therefore, promising for estimation of VTFs using realizations obtained from computer simulations of more challenging systems. The present method can also be extended to estimate cubic transfer functions for more accurate Volterra models. Note that the approach can be easily utilized with multiple realizations on parallel PC cluster systems or applied across a GRID composed of heterogeneous clusters.