Nonparametric System Identification: Impulse & Frequency Response

Nonparametric identification estimates an impulse response or frequency-response function without first compressing the dynamics into a small pole-zero model. It is the diagnostic layer between raw experiments and parametric fitting, exposing delay, resonances, usable bandwidth, uncertainty and model violations.

A low-order transfer function answers “which compact model will I use?” A nonparametric estimate answers an earlier question: “what input/output dynamics are actually visible in these data?” It spends many degrees of freedom on an impulse sequence or frequency grid, gaining diagnostic resolution and losing the extrapolation, smoothness and compactness supplied by a parametric structure.

Time-domain and frequency-domain routes to nonparametric system identification Measured input and output records feed either a convolution regression producing impulse coefficients or spectral products producing a frequency response. Both estimates carry uncertainty and are compared on independent runs before selecting a compact parametric model. records u[k], y[k] excited · aligned conditioned Toeplitz regression y ≈ U g horizon + regularization impulse estimate g[0 … n] records u[k], y[k] periods · windows averages DFT / spectra Y/U or Sᵧᵤ/Sᵤᵤ smoothing sets resolution frequency response G(eʲω) excited band only uncertainty and checks input power · coherence repeatability · causality independent validation then choose a compact model time-domain route frequency-domain route
Time-domain estimation solves for impulse coefficients; frequency-domain estimation forms or smooths spectral ratios. Both require excitation and uncertainty checks before a parametric model is fitted.

The object being estimated #

For a stable causal discrete-time SISO LTI system,

yk=∑i=0∞giuk−i+vk,y_k=\sum_{i=0}^{\infty}g_i u_{k-i}+v_k,

where gig_i is the sampled impulse response and vkv_k collects output disturbance and measurement noise. The frequency response is its discrete-time Fourier transform,

G(ejω)=∑i=0∞gie−jωi.G(e^{j\omega})=\sum_{i=0}^{\infty}g_i e^{-j\omega i}.

These are two representations of the same linear operator. Finite data force approximations: truncate the impulse response, sample a finite frequency grid, assume stationarity over the record, and quantify which coefficients or frequencies are weakly determined.

Nonparametric does not mean assumption-free. It still assumes a causal channel, a chosen input/output pairing, time invariance over the experiment, sufficient excitation, and a disturbance relationship compatible with the estimator.

Direct impulse-response estimation #

Truncate after ngn_g coefficients:

yk≈g0uk+g1uk−1+⋯+gng−1uk−ng+1+vk.y_k\approx g_0u_k+g_1u_{k-1}+\cdots+g_{n_g-1}u_{k-n_g+1}+v_k.

Stacking samples gives a linear regression

y=Ug+v,\mathbf y=\mathbf U\mathbf g+\mathbf v,

where each row of the Toeplitz matrix U\mathbf U contains delayed input samples. Ordinary least squares gives

g^=(UTU)−1UTy\hat{\mathbf g}=(\mathbf U^T\mathbf U)^{-1}\mathbf U^T\mathbf y

when the matrix is full rank and the noise assumptions justify it. Regularization or smoothness penalties can stabilize a long FIR, but then the prior is part of the estimate and must be reported.

The horizon ngTsn_gT_s should cover the response tail relevant to the use. Too short a horizon forces omitted slow dynamics into the last coefficients and residual. Too long a horizon estimates many noise-dominated coefficients with high variance. Plot coefficient confidence intervals and cumulative step response, not only the mean curve.

Correlation analysis #

For zero-mean stationary signals with input independent of additive output noise,

Ryu[ℓ]=∑igiRuu[ℓ−i].R_{yu}[\ell]=\sum_i g_i R_{uu}[\ell-i].

If the input is ideally white, Ruu[ℓ]=σu2δ[ℓ]R_{uu}[\ell]=\sigma_u^2\delta[\ell], so

gℓ=Ryu[ℓ]σu2.g_\ell=\frac{R_{yu}[\ell]}{\sigma_u^2}.

With colored input, estimating gg requires deconvolving by the input autocorrelation or solving the corresponding normal equations. Finite records replace expectations with sample correlations, creating edge effects and uncertainty.

The independence assumption breaks in feedback: output disturbance changes the measured output, the controller reacts, and the resulting input is correlated with that disturbance. A naive correlation estimate then mixes plant dynamics with feedback-induced correlation. External excitation and closed-loop identification methods are needed; subtracting a mean does not fix causality.

Step response is an integral view #

For a unit step input, the output of a stable discrete LTI system is

sk=∑i=0kgi.s_k=\sum_{i=0}^{k}g_i.

Therefore gk=sk−sk−1g_k=s_k-s_{k-1} in the ideal discrete case. Differencing a noisy measured step response amplifies high-frequency noise, so estimating gg by visually differentiating a plot is fragile. A step is excellent for DC gain, dominant delay and time constants, but its spectrum decays with frequency and often leaves fast modes poorly determined.

Empirical transfer-function estimate #

Take DFTs of a finite input/output record. At an excited frequency bin,

G^ETFE(ωm)=Y(ωm)U(ωm).\hat G_{ETFE}(\omega_m)=\frac{Y(\omega_m)}{U(\omega_m)}.

This is the empirical transfer-function estimate (ETFE). It preserves narrow resonances and is especially effective for periodic inputs observed over complete periods: repeat periods can be averaged, and exact excited lines avoid leakage when synchronization is correct.

For general noisy data, ETFE is jagged because each bin divides by one finite-record input coefficient. Where ∣U∣|U| is small, output noise is magnified. A longer record creates more frequency lines but does not by itself make every unsmoothed line low-variance. Never interpret an ETFE outside the excited band, across spectral notches, or where the input coefficient is near zero.

Spectral ratio estimates #

For an open-loop experiment with input independent of additive output noise,

Syu(ω)=G(ejω)Suu(ω),S_{yu}(\omega)=G(e^{j\omega})S_{uu}(\omega),

so a common estimate is

G^(ω)=S^yu(ω)S^uu(ω).\hat G(\omega)=\frac{\hat S_{yu}(\omega)}{\hat S_{uu}(\omega)}.

Spectral densities are estimated by averaging or smoothing finite-record Fourier products. Blackman–Tukey methods window estimated covariance sequences before transformation. Welch methods split the record into overlapping windowed segments and average their periodograms. Both trade frequency resolution for reduced variance.

The tuning parameter—segment length, covariance window or frequency-dependent smoothing width—is part of the model:

  • wider smoothing averages more independent fluctuation and lowers variance;
  • wider smoothing also merges nearby resonances and biases sharp dynamics;
  • short windows broaden spectral features through leakage;
  • long windows preserve resolution but leave fewer averages.

Use logarithmic or frequency-dependent resolution only with an explicit reason. A constant fractional resolution can be useful across decades; a narrow mechanical resonance may require locally finer resolution than a slow process pole.

Leakage, transients and periodic experiments #

The DFT treats each record as one period of a periodic signal. If record endpoints do not join, energy leaks across bins. A taper reduces the discontinuity but also changes amplitude, bandwidth and the number of effectively independent samples.

Periodic multisine and PRBS experiments provide a stronger alternative when feasible:

  1. drive the plant to periodic steady state;
  2. record an integer number of periods with synchronized sampling;
  3. average DFT coefficients across periods;
  4. estimate noise variance from period-to-period deviations;
  5. reserve unexcited detection lines to reveal nonlinear distortion or leakage.

Discarding initial periods removes transient contamination only when the plant is stable and the discarded interval is long enough. For an unstable plant or closed-loop experiment, initial conditions and controller action must be modeled rather than wished away.

Coherence is a diagnostic, not a verdict #

Magnitude-squared coherence is

γyu2(ω)=∣Syu(ω)∣2Suu(ω)Syy(ω),0≤γ2≤1.\gamma_{yu}^2(\omega)=\frac{|S_{yu}(\omega)|^2}{S_{uu}(\omega)S_{yy}(\omega)},\qquad 0\le\gamma^2\le1.

Under the estimator's assumptions, low coherence can indicate poor input signal-to-noise ratio, unmeasured disturbances, nonlinear response, time variation, leakage or timing error. It does not identify which cause is responsible. High coherence says a repeatable linear relationship dominates at that frequency; it does not prove the input caused the output, that the loop is open, or that the model extrapolates to another operating point.

For MIMO data, pairwise coherence can be misleading when inputs are correlated. Use multivariable spectral matrices and examine their conditioning; two nearly collinear inputs cannot separately identify two channels even if each appears coherent with the output.

Reading uncertainty #

A useful nonparametric result carries uncertainty with each time lag or frequency. Report at least:

  • confidence bands or covariance of gkg_k and G(ejω)G(e^{j\omega});
  • input spectrum and excited frequencies;
  • coherence or an equivalent residual/noise diagnostic;
  • number and length of averages, window and overlap;
  • sample period, frequency resolution and retained band;
  • repeated-run and operating-point variability.

Pointwise confidence intervals do not automatically form a simultaneous band over all frequencies. Inspecting hundreds of bins guarantees some apparent excursions by chance. For control design, uncertainty should be converted into a frequency-dependent plant-error description over the intended bandwidth, then checked against sensitivity and robustness requirements.

What the shape tells you #

Observation Plausible mechanism Next check
Approximately linear phase with frequency Pure delay Recheck channel timestamps and command/sensor latency
Narrow magnitude peak with phase rotation Lightly damped mode Repeat with finer resolution and lower amplitude; inspect nonlinear distortion
High-frequency gain floor Measurement noise divided by weak input Plot SuuS_{uu} and uncertainty; do not fit poles to the floor
Different estimates across periods Time variation, nonlinear behavior or nonstationary disturbance Separate periods/operating states and inspect amplitude dependence
Low coherence only near one band Weak excitation, local disturbance or nonlinearity Redesign input energy and use detection lines/repeated runs
Noncausal-looking impulse energy Timing error, feedback correlation or estimator/window artifact Audit alignment and causal architecture before adding model delay

These are hypotheses, not automatic classifications. Change one experimental condition and repeat the estimate before turning a visual feature into a pole, zero or delay.

From nonparametric evidence to a model #

  1. Establish the usable frequency band from excitation and uncertainty.
  2. Estimate delay, DC gain, resonances, roll-off and likely order without forcing a rational form.
  3. Compare time-domain impulse and frequency-domain estimates; they should tell a compatible story within their windows and uncertainty.
  4. Select the smallest parametric structure able to represent the decision-relevant features.
  5. Fit on estimation data and compare its response against the nonparametric estimate.
  6. Validate simulation, prediction and residuals on independent runs.

Nonparametric identification is not merely a preliminary plot. It is the empirical constraint that prevents a convenient parametric model from inventing certainty outside the information carried by the experiment.

Sources

  1. paperPrediction Error Estimation Methodspdfnotes
  2. docsMathWorks: Selecting the Method for Computing Spectral Models
  3. docsMathWorks: Nonparametric Impulse Response Estimation
  4. paperWelch — The Use of Fast Fourier Transform for the Estimation of Power Spectra

Searches titles, aliases, headings, summaries, and full text.↑↓ move ↵ open esc close