Chapter 03

Identification in the Frequency Domain

The other black-box route. Instead of one experiment and a Hankel matrix, run H single-sinusoid experiments and turn each into one point of the frequency response, then fit a transfer function to those points. The dataset is manufactured rather than collected — and the payoff is a model whose error is small exactly where a control design will look.

Reading: ~26 min Interactive: 1 widgets Source: Polimi MIDA2 2025/26 — Lecture 17 (Ch. 2.2 — frequency-domain identification) · Polimi MIDA2 2025/26 — Lecture 21, slides 1–7 (Ch. 2.2 — closing remarks)

01 · Motivation

The other black-box route

Chapter 02 identified a model from a single experiment. This chapter identifies one from $H$ experiments, gets a transfer function instead of a state-space model, and costs far more test time. It exists because of where the resulting model is accurate.

Both are black-box methods; they differ in almost everything else:

4SID — time domain

One experiment. Non-parametric and constructive: no model family, no cost function, no iteration. Returns $\{F,G,H,D\}$ — a state-space model, which is what a Kalman filter needs.

Frequency domain

$H$ experiments, one per frequency. A classic parametric method: choose a family $M(\theta)$ , define an index $J(\theta)$ , optimise. Returns $W(z;\hat\theta)$ — a transfer function, which is what a frequency-domain control design needs.

The method follows the standard four-step parametric recipe:

Experiment design, data collection, pre-processing
Choice of a parametric family $M(\theta)$
Definition of a performance index $J(\theta)$
Optimisation, $\hat\theta=\arg\min_\theta J(\theta)$

Steps ② and ④ are shared with every other parametric method in the course. Only ① and ③ are special here, and they are what this chapter is about: step ① manufactures a dataset of complex numbers instead of collecting a time series, and step ③ measures error in the complex plane.

why

Why bother, if the model is the same?

In principle the transfer function you obtain is the same one a time-domain method would give. The reason to prefer this route is that identification, like any fitting, distributes its error — and this method distributes it across the frequency grid you chose. Many control techniques are designed in the frequency domain, and a controller is only as good as the model’s accuracy in the band where the loop will operate. Choosing the grid is choosing where to be right.

tip

The running example — a car's yaw dynamics

The lectures hang the method on vehicle dynamics: the input is the steering angle $\delta$ , the output is the yaw rate $\omega_z$ (the car’s rotational speed about its vertical axis). A model of $\delta\to\omega_z$ is what Electronic Stability Control and the lateral-dynamics controllers of autonomous cars are designed against — and those designs are done in the frequency domain, which is why the model is identified there too.

02 · The core skill

Estimating amplitude and phase from a noisy sinusoid

The most repeated question in the bank

“Describe the procedure to estimate a sinusoidal signal — its amplitude and phase — from a noisy one, assuming the frequency is known a priori.” That exact stem has been set in four separate sittings (08/01/2024, 04/09/2024, 10/02/2025 and 26/06/2026), more than any other question in the archive. It is a pure theory question worth 3 marks and it wants the reasoning below, not a number.

Note that the question bank deduplicates those four into one entry, which is why the automatic frequency analysis under-reads this chapter’s weight — the reason it carries a manual override.

Strip away the context and the problem is self-contained. You measure

The measurement

$y(t)=B\sin(\Omega t+\varphi)+\text{noise},\qquad t=1,\dots,N,$

with the frequency $\Omega$ known. Recover $B$ and $\varphi$ .

The one idea that matters

The obvious move is to fit $B\sin(\Omega t+\varphi)$ directly by least squares. Do not. The phase sits inside the sine, so the model is nonlinear in $\varphi$ : the cost surface is non-convex, the solution needs iteration, and it has local minima.

Instead, use the trigonometric identity in reverse. Since $B\sin(\Omega t+\varphi)=B\cos\varphi\,\sin(\Omega t)+B\sin\varphi\,\cos(\Omega t)$ , the very same signal can be written

The linear parameterisation

$\hat y(t)=a\sin(\Omega t)+b\cos(\Omega t),\qquad a=B\cos\varphi,\quad b=B\sin\varphi.$

key

Same signal, different unknowns

$\{B,\varphi\}$ and $\{a,b\}$ describe exactly the same set of sinusoids — nothing has been approximated. But $\{a,b\}$ enter linearly, because $\sin(\Omega t)$ and $\cos(\Omega t)$ are fixed known signals once $\Omega$ is known. A model linear in its parameters has a quadratic cost, and a quadratic cost has a closed-form minimum.

Re-parameterising to make a fit linear is a move worth recognising in general. Here it converts an iterative optimisation with local minima into a $2\times2$ linear solve.

Worked derivation The closed-form estimate

Write the quadratic cost

$J_N(a,b)=\frac1N\sum_{t=1}^{N}\Big(y(t)-\big(a\sin(\Omega t)+b\cos(\Omega t)\big)\Big)^{2}$

This is the sample variance of the estimation error — the same performance index shape used everywhere in parametric identification.

Set both partial derivatives to zero

$\frac{\partial J_N}{\partial a}=0,\qquad \frac{\partial J_N}{\partial b}=0$

The factor $-\tfrac2N$ cancels from both, leaving

$\sum_{t=1}^{N}\sin(\Omega t)\Big(y(t)-a\sin(\Omega t)-b\cos(\Omega t)\Big)=0,$ $\sum_{t=1}^{N}\cos(\Omega t)\Big(y(t)-a\sin(\Omega t)-b\cos(\Omega t)\Big)=0.$

Collect into the normal equations

$\begin{bmatrix}\sum\sin^{2}(\Omega t) & \sum\sin(\Omega t)\cos(\Omega t)\\[4pt] \sum\sin(\Omega t)\cos(\Omega t) & \sum\cos^{2}(\Omega t)\end{bmatrix}\begin{bmatrix}a\\ b\end{bmatrix}=\begin{bmatrix}\sum y(t)\sin(\Omega t)\\[4pt] \sum y(t)\cos(\Omega t)\end{bmatrix}$

A $2\times2$ linear system whose matrix depends only on $\Omega$ and $N$ — not on the data. Invert it once:

$\begin{bmatrix}\hat a\\ \hat b\end{bmatrix}=\big(\text{matrix}\big)^{-1}\big(\text{vector}\big).$

Convert back to amplitude and phase

From $\hat a=\hat B\cos\hat\varphi$ and $\hat b=\hat B\sin\hat\varphi$ , dividing the second by the first gives

$\hat\varphi=\arctan\!\left(\frac{\hat b}{\hat a}\right),\qquad \hat B=\sqrt{\hat a^{2}+\hat b^{2}}.$

✓

Two forms of the amplitude formula

The lectures write the amplitude as the average of the two individually recovered values,

$\hat B=\frac12\left(\frac{\hat a}{\cos\hat\varphi}+\frac{\hat b}{\sin\hat\varphi}\right),$

whereas most textbooks write $\hat B=\sqrt{\hat a^{2}+\hat b^{2}}$ . They are the same number. Both bracket equations hold exactly by construction, so each term in the average already equals $\hat B$ , and squaring and adding $\hat a=\hat B\cos\hat\varphi$ , $\hat b=\hat B\sin\hat\varphi$ gives $\hat a^{2}+\hat b^{2}=\hat B^{2}$ . Either is acceptable; the square root is safer, since the averaged form divides by $\cos\hat\varphi$ or $\sin\hat\varphi$ , which can be zero.

Three ways this goes wrong

Fitting $\{B,\varphi\}$ directly. The whole point of the $\{a,b\}$ form is to avoid a nonlinear fit. An answer that sets up an iterative optimisation has missed the question.
Forgetting that $\Omega$ known is load-bearing. If the frequency were also unknown, the regressors $\sin(\Omega t)$ and $\cos(\Omega t)$ would themselves depend on a parameter and the problem would be nonlinear again. State this explicitly — it is worth marks.
Not resolving the quadrant of $\arctan$ . The ratio $\hat b/\hat a$ is unchanged when both signs flip, so the arctangent alone leaves a $\pi$ ambiguity. The individual signs of $\hat a$ and $\hat b$ pick the half-plane.

why

Why this rejects noise

Averaging $N$ samples against a known-frequency reference is correlation, and broadband noise is asymptotically uncorrelated with a pure sinusoid. The noise contributions to the right-hand side sums partially cancel while the signal contributions accumulate, so the estimate improves like $\sqrt N$ . Running the experiment longer buys accuracy — which is precisely why the method can afford to work one frequency at a time.

2024-01-08-q42024Q04Sinusoid amplitude and phase estimationmedium3 pts

Describe the procedure to estimate a sinusoidal signal (estimation of its amplitude and phase) from a noisy one, assuming that the frequency $\Omega$ of the sinusoid is a-priori known.

03 · Step ①

Designing the experiment

Step ① is asked as a question in its own right

“Describe in detail step 1 (experiment design and pre-processing) of a frequency-domain identification method” — 3 marks, set in January 2024. Together with the sinusoid question above, this chapter’s two stems span 5 of the 12 sittings in the archive. A full-credit answer covers the grid, the amplitudes, the linearity check, and the pre-processing fit, and ends by saying what step ① produces.

Choosing the frequency grid

Select $H$ excitation frequencies $\omega_1,\dots,\omega_H$ , usually uniformly spaced by a constant $\Delta\omega$ . The upper edge is bounded above by the Nyquist frequency $\omega_N$ , and in practice you want real headroom, not the bare bound: $\omega_H \ll \omega_N$ .

The lectures give a complete design chain, anchored on the expected closed-loop bandwidth $\omega_c$ of the controller you intend to build:

The design chain

$\omega_H\approx 3\,\omega_c,\qquad \omega_N\approx 10\,\omega_c,\qquad \omega_s=2\,\omega_N.$

Identify three times past the control bandwidth, sample ten times past it, then double for the sampling frequency.

tip

Worked: stability control in a car

The expected control bandwidth is $f_c=4$ Hz. The chain gives

$f_H\approx 12\ \text{Hz},\qquad f_N=40\ \text{Hz},\qquad f_s=80\ \text{Hz}.$

So: excite up to 12 Hz, sample at 80 Hz. The 12-to-40 Hz margin is what stops measurement noise and unmodelled high-frequency dynamics from folding back into the identified band.

Choosing the amplitudes

The amplitudes $A_1,\dots,A_H$ may all be equal, but usually they decrease with frequency, and the reason is actuator power rather than anything statistical.

Derivation Where the amplitude schedule comes from

Model the actuator effort

For the steering example, torque is proportional to angle, $T(t)=K\delta(t)$ , and mechanical power is torque times rotational speed:

$\text{power}=K\,\delta(t)\,\dot\delta(t).$

Substitute the sinusoidal input

With $\delta(t)=A_i\sin(\omega_i t)$ , so $\dot\delta(t)=A_i\omega_i\cos(\omega_i t)$ :

$\text{power}=\underbrace{K A_i^{2}\,\omega_i}_{\text{peak demand}}\sin(\omega_i t)\cos(\omega_i t)$

The demand scales with the square of amplitude and linearly with frequency.

Impose a constant power budget

If the actuator has a power limit, hold $K A_i^{2}\omega_i$ constant across the grid:

$A_i=\sqrt{\frac{c}{K\,\omega_i}}$

The excitation amplitude must fall like $1/\sqrt{\omega_i}$ — large and slow at the bottom of the grid, small and fast at the top.

A constant amplitude fails at both ends

Pick one amplitude for the whole grid and you must size it for the worst case, which is the highest frequency. That leaves the low-frequency experiments excited far below what the actuator could deliver — wasting signal-to-noise exactly where the control bandwidth is. Size it for the low end instead and the high-frequency experiments saturate the actuator, which makes the plant nonlinear and voids the method.

Running them

Perform $H$ independent experiments, one sinusoid each:

$\text{experiment } i:\quad u_i(t)=A_i\sin(\omega_i t)\ \longrightarrow\ \text{system}\ \longrightarrow\ y_i(t)$

each yielding a dataset $\{u_i(1),\dots,u_i(N)\}$ , $\{y_i(1),\dots,y_i(N)\}$ . Wait for the transient to die before recording: the frequency response is a steady-state property.

Do not excite all frequencies at once

A broadband input would collect all the data in one run and seems obviously better. It is not what this method does, and two things break if you try. The per-frequency fit stops being a two-parameter closed-form problem, and the linearity test of §04 becomes impossible — with energy already at every frequency, you can no longer tell whether output energy at $2\omega_i$ came from the plant or from your own input. One sinusoid at a time is what buys both.

The honest cost is test time: $H$ experiments means $H$ separate runs of the plant. That is the real counterweight to 4SID’s single dataset.

2024-01-29-q42024Q04Frequency-domain experiment designmedium3 pts

Describe in detail the "step 1" (experiment design and pre-processing) of a frequency-domain system identification method.

04 · A free diagnostic

The frequency response theorem as a linearity test

The whole method assumes the plant is linear and time-invariant. Unusually, the experiment that collects the data also checks that assumption — at no extra cost.

key

Frequency response theorem

For an LTI system, if the input is a sinusoid of frequency $\omega_i$ , the steady-state output must be a sinusoid of the same frequency $\omega_i$ — differing only in amplitude and phase.

The theorem is normally used as a computational tool. Turn it around and it becomes a diagnostic: an LTI system cannot create energy at a frequency that was not in its input. So if you inject a clean tone at $\omega_i$ and the output contains a visible component at $2\omega_i$ or $3\omega_i$ , the plant is not LTI, and every result downstream is void.

✓

Two ways an output fails to be a perfect sinusoid — only one matters

Real outputs are never exactly sinusoidal, for two quite different reasons:

Noise — measurement noise, internal disturbances. Broadband, uncorrelated with the input, and precisely what §02’s least-squares fit is designed to average away.
Nonlinearity — saturation, friction, hysteresis. This shows up as energy at harmonics of $\omega_i$ , and it is not noise: no amount of averaging removes it, because it is deterministically locked to the input.

Small harmonics are tolerated — the method is explicitly looking for “the best LTI approximation” of a mildly nonlinear plant. Large ones mean the method does not apply and you should be linearising around an operating point instead.

You inject a 5 Hz sinusoid and the measured output has a strong component at 5 Hz and a clearly visible one at 10 Hz. What should you conclude?

05 · Step ①, output

One experiment, one point of the response

Each experiment has now been reduced to a clean amplitude and phase. Turning that pair into a piece of the frequency response takes one division and one exponential.

If the output is not a pure sinusoid at the input frequency, the plant is not linear — the same experiment is the test.

Experiment $i$ injected $A_i\sin(\omega_i t)$ and, after cleaning, returned $\hat B_i\sin(\omega_i t+\hat\varphi_i)$ . Comparing input and output:

quantity	value
amplitude gain	$\hat B_i / A_i$
phase shift	$\hat\varphi_i - 0 = \hat\varphi_i$

and a gain together with a phase is a complex number. The estimated frequency response at $\omega_i$ is therefore

One frequency-response point

$\hat W(e^{\,j\omega_i})=\frac{\hat B_i}{A_i}\,e^{\,j\hat\varphi_i}.$

Note that the amplitude is a ratio — dividing by $A_i$ is what undoes the amplitude schedule of §03, so the decreasing amplitudes leave no trace in the result.

Repeat $H$ times and step ① is finished. Its output is

$\left\{\frac{\hat B_1}{A_1}e^{\,j\hat\varphi_1},\ \frac{\hat B_2}{A_2}e^{\,j\hat\varphi_2},\ \dots,\ \frac{\hat B_H}{A_H}e^{\,j\hat\varphi_H}\right\}$

— $H$ complex numbers indexed by frequency. Plotted as magnitude and phase against $\omega$ , they are measured points of a Bode plot.

key

The dataset was manufactured, not collected

This is the conceptual payoff of the chapter. In time-domain identification the dataset is the raw $\{u(t),y(t)\}$ you recorded. Here step ① transforms the raw recordings into an entirely different object: a set of complex numbers indexed by frequency, with the time axis gone.

Everything downstream — choose a model family, define an index, optimise — is then ordinary parametric identification against that new dataset. This is why only steps ① and ③ are described as “special”: they are the only two that know the data are complex.

06 · Steps ②–④

A complex-plane index that is not quadratic

With a dataset in hand, the remaining three steps are the familiar parametric ones — with one consequential difference in the performance index.

Step ② — choose a parametric transfer-function family:

The model family

$W(z;\theta)=\frac{b_0+b_1z^{-1}+\cdots+b_pz^{-p}}{1+a_1z^{-1}+\cdots+a_nz^{-n}}\cdot z^{-1},\qquad \theta=\begin{bmatrix}a_1&\cdots&a_n&b_0&\cdots&b_p\end{bmatrix}^{\top}.$

Step ③ — define the performance index. Model and measurement are both complex numbers at each grid frequency, so the error is a distance in the complex plane:

The frequency-domain index

$J_H(\theta)=\frac1H\sum_{i=1}^{H}\left|\,W(e^{\,j\omega_i};\theta)-\frac{\hat B_i}{A_i}e^{\,j\hat\varphi_i}\,\right|^{2}$

— the variance of the frequency-response modelling error. Each term is the squared length of the segment joining the modelled point to the measured one, so the index penalises gain and phase errors together, in one number, without ever choosing a weighting between them.

Step ④ — optimise, $\hat\theta_H=\arg\min_\theta J_H(\theta)$ , giving the final model $\hat W(z;\hat\theta_H)$ .

Two optimisations, and only one of them is easy

This chapter contains two least-squares problems and they behave completely differently. Conflating them is the classic error.

The inner fit (§02), estimating $\{a_i,b_i\}$ for one experiment, is quadratic in its parameters — a $2\times2$ solve, one shot, no iteration, global optimum guaranteed.
The outer fit (here), estimating $\theta$ , is not quadratic. The parameters $a_1,\dots,a_n$ sit in the denominator of $W(z;\theta)$ , so the model is a nonlinear function of them and the cost surface has no closed-form minimum.

Step ④ therefore needs iterative numerical optimisation — gradient, Newton, quasi-Newton, or randomised methods — with all that this implies: an initial guess, convergence that is not guaranteed, and the possibility of local minima.

Deep dive Why the denominator is what breaks quadraticity

It is worth being precise about which parameters cause the trouble, because the numerator ones do not.

Write $W(z;\theta)=B(z;b)/A(z;a)$ . Holding $a$ fixed, $W$ is a linear function of the numerator coefficients $b$ — each $b_k$ multiplies a fixed known quantity $z^{-k}/A(z;a)$ . If $a$ were known, the whole problem would collapse back to linear least squares in $b$ .

The denominator coefficients enter through a reciprocal, and $1/A(z;a)$ is not linear in $a$ — nor is its magnitude-squared, so no re-parameterisation of the §02 kind rescues it. The residual is a ratio of polynomials in the unknowns, and the index is a sum of squared magnitudes of such ratios: smooth, but neither convex nor quadratic.

This asymmetry is the reason a family of classical methods exists that alternate between a linear solve for $b$ and an update of $a$ — and it is the same structural difficulty as the ARMAX identification of MIDA1, where the moving-average part is what forces iteration while the ARX case stays closed-form.

07 · In practice

Weighting, non-uniform grids, and unstable plants

Three refinements that come up whenever the method meets a real plant.

Focusing the fit where it matters

Sometimes one band matters far more than the rest — a lightly damped resonance the controller must not excite, say. A fit that spreads error uniformly will smooth right over such a peak. Two levers correct this, and they are worth contrasting:

Weight the index

Attach a weight $\gamma_i$ to each grid point, giving

$\tilde J_H(\theta)=\frac1H\sum_{i=1}^{H}\gamma_i\left|W(e^{\,j\omega_i};\theta)-\frac{\hat B_i}{A_i}e^{\,j\hat\varphi_i}\right|^{2}$

with $\gamma_i=1$ outside the band of interest and larger inside it. Applied after the fact, costs nothing, and reuses the data you already have — but it cannot add information the experiments did not gather.

Redesign the grid

Drop the uniform $\Delta\omega$ and concentrate experiments in the band of interest. Decided a priori and costs experiment time — but it genuinely collects more information there, which weighting cannot fake.

The unweighted index is the special case $\gamma_i\equiv1$ .

Identifying an unstable plant

An open-loop experiment on an unstable system is impossible: drive it with a sinusoid and the output diverges, so there is nothing to fit. The fix is to stabilise it first.

✓

Closed-loop identification

Design a simple stabilising controller, close the loop, and drive the reference $y^{\circ}(t)$ with the sinusoid. The loop keeps the plant bounded, and you record $u(t)$ and $y(t)$ inside the loop — the plant’s own input and output. Those two datasets go into the identification exactly as before.

The controller is scaffolding, present only so the experiment can run. It is not part of what you identify, and the resulting model is a model of the open-loop plant — poles outside the unit circle and all.

Expecting a stable model out of a closed-loop experiment

The measured signals are well behaved, so it is tempting to expect the identified model to be stable too. It should not be, and if it is, something has gone wrong — most likely you identified the closed-loop transfer function by feeding the algorithm $y^{\circ}$ instead of $u$ . The plant is unstable; a correct identification says so.

Load-bearing ideas

The four-step recipe with only ① and ③ special: step ① manufactures a dataset of complex numbers, step ③ measures error in the complex plane.
The parameterisation trick: fit $a\sin(\Omega t)+b\cos(\Omega t)$ , never $B\sin(\Omega t+\varphi)$ . Linear in $\{a,b\}$ ⟹ quadratic cost ⟹ a $2\times2$ closed-form solve. Recover $\hat\varphi=\arctan(\hat b/\hat a)$ and $\hat B=\sqrt{\hat a^{2}+\hat b^{2}}$ . $\Omega$ known is what makes this linear.
Design chain: $\omega_H\approx3\omega_c$ , $\omega_N\approx10\omega_c$ , $\omega_s=2\omega_N$ ; and the amplitude schedule $A_i=\sqrt{c/(K\omega_i)}$ from a constant actuator-power budget.
The frequency response theorem doubles as a linearity test — harmonics in the output mean the plant is not LTI, and no amount of averaging will fix that.
One experiment, one point: $(\hat B_i/A_i)e^{\,j\hat\varphi_i}$ . $H$ of them make a measured Bode plot.
The outer index is not quadratic — the denominator coefficients enter nonlinearly — so step ④ needs iterative optimisation, unlike the one-shot inner fit.
Refinements: weight $\gamma_i$ post-hoc or redesign the grid a priori; and identify an unstable plant in closed loop, remembering that the model that comes out is still unstable.

Exam radar. This chapter is tested as theory, not calculation — both its questions are “describe the procedure” prompts worth 3 marks each, and between them they span 5 of the 12 sittings. Be able to write §02 end to end from memory: the two parameterisations, why one is linear, the normal equations, and the conversion back. Then be able to list step ①’s four decisions — grid, amplitudes, independent experiments, sinusoid fit — and say what step ① produces.