For a stable LTI system with transfer function \( G(s) \), the frequency response is defined as the complex function

\[ G(\mathrm{j}\omega) = G(s)\big|_{s=\mathrm{j}\omega}, \quad \omega \in \mathbb{R}_{\ge 0}. \]

In previous lessons, we obtained \( G(\mathrm{j}\omega) \) analytically from a known transfer function. In practice, especially for robotic joints, flexible structures, or poorly known actuators, only experimental input–output data \( u(t), y(t) \) are available. The goal of this lesson is:

\[ \text{Given sampled data } \{u(kT_s),y(kT_s)\}_{k=0}^{N-1},\ \text{estimate } G(\mathrm{j}\omega_\ell) \ \text{for selected frequencies } \omega_\ell. \]

A canonical experiment uses a sinusoidal input at a single angular frequency \( \omega \):

\[ u(t) = U_0 \sin(\omega t), \quad t \ge 0, \]

and after transients, the steady-state output of an LTI system is also sinusoidal at the same frequency, with a different amplitude and phase:

\[ y_{\text{ss}}(t) = Y_0 \sin(\omega t + \varphi(\omega)). \]

The experimentally measured frequency response at \( \omega \) is then

\[ \hat{G}(\mathrm{j}\omega) = \frac{Y_0}{U_0}\,\mathrm{e}^{\mathrm{j}\varphi(\omega)}. \]

Repeating this measurement for a grid of frequencies yields discrete samples of \( G(\mathrm{j}\omega) \), which can later be plotted or used for controller design.

Assume we apply a sinusoidal input with known amplitude \( U_0 \) and angular frequency \( \omega \):

\[ u(t) = U_0 \sin(\omega t). \]

For a stable LTI system, the steady-state output can be written as

\[ y_{\text{ss}}(t) = A\cos(\omega t) + B\sin(\omega t), \]

for some real coefficients \( A \) and \( B \). The complex steady-state output phasor is

\[ Y(\mathrm{j}\omega) = B - \mathrm{j}A \quad\text{(up to a conventional sign choice)}. \]

If we denote the input phasor by \( U(\mathrm{j}\omega) \), then by definition

\[ G(\mathrm{j}\omega) = \frac{Y(\mathrm{j}\omega)}{U(\mathrm{j}\omega)}. \]

In an experiment, we only observe noisy samples \( y(kT_s) \), \( k = 0,\dots,N-1 \), where \( T_s \) is the sampling period. We model

\[ y(kT_s) = A\cos(\omega kT_s) + B\sin(\omega kT_s) + v(kT_s), \]

where \( v(kT_s) \) is additive measurement noise with zero mean. Define \( \phi_1(k) = \cos(\omega kT_s) \) and \( \phi_2(k) = \sin(\omega kT_s) \). We can write in vector form

\[ \mathbf{y} = \mathbf{\Phi}\boldsymbol{\theta} + \mathbf{v}, \quad \boldsymbol{\theta} = \begin{bmatrix} A \\ B \end{bmatrix},\quad \mathbf{\Phi} = \begin{bmatrix} \phi_1(0) & \phi_2(0)\\ \phi_1(1) & \phi_2(1)\\ \vdots & \vdots\\ \phi_1(N-1) & \phi_2(N-1) \end{bmatrix}. \]

The least-squares estimator (LS) of \( \boldsymbol{\theta} \) is

\[ \hat{\boldsymbol{\theta}} = \begin{bmatrix} \hat{A} \\ \hat{B} \end{bmatrix} = \left(\mathbf{\Phi}^\top \mathbf{\Phi}\right)^{-1} \mathbf{\Phi}^\top \mathbf{y}. \]

If the sampling interval covers an integer number of input periods and the noise is zero-mean and uncorrelated with the regressors, one can show that the LS estimator is unbiased, i.e. \( \mathbb{E}[\hat{A}] = A \), \( \mathbb{E}[\hat{B}] = B \).

Once we have \( \hat{A} \) and \( \hat{B} \), we construct the complex output phasor estimate as

\[ \hat{Y}(\mathrm{j}\omega) = \hat{B} - \mathrm{j}\hat{A}, \]

and similarly obtain the input phasor (for example, from \( u(t) = U_0 \sin(\omega t) \) we have \( U(\mathrm{j}\omega) = U_0/\mathrm{j} \) up to phase), giving the experimental estimate

\[ \hat{G}(\mathrm{j}\omega) = \frac{\hat{Y}(\mathrm{j}\omega)}{U(\mathrm{j}\omega)}. \]

From \( \hat{G}(\mathrm{j}\omega) \), the magnitude and phase are

\[ \left|\hat{G}(\mathrm{j}\omega)\right| = \sqrt{\left(\operatorname{Re}\hat{G}(\mathrm{j}\omega)\right)^2 +\left(\operatorname{Im}\hat{G}(\mathrm{j}\omega)\right)^2},\quad \angle \hat{G}(\mathrm{j}\omega) = \operatorname{atan2} \left(\operatorname{Im}\hat{G}(\mathrm{j}\omega), \operatorname{Re}\hat{G}(\mathrm{j}\omega)\right). \]

In practice, it is often convenient to treat the measured signals as finite-length discrete-time sequences and use the Discrete Fourier Transform (DFT). Let \( u[k] = u(kT_s) \), \( y[k] = y(kT_s) \), with \( k = 0,\dots,N-1 \). The DFTs are

\[ U[m] = \sum_{k=0}^{N-1} u[k]\mathrm{e}^{-\mathrm{j}\frac{2\pi mk}{N}},\quad Y[m] = \sum_{k=0}^{N-1} y[k]\mathrm{e}^{-\mathrm{j}\frac{2\pi mk}{N}}, \]

where the discrete frequency index \( m \) corresponds to the analog angular frequency

\[ \omega_m = \frac{2\pi m}{N T_s}, \quad m = 0,\dots,N-1. \]

If the data window contains an integer number of periods of the sinusoid and leakage effects are negligible (or windows are used appropriately), we can estimate the frequency response at \( \omega_m \) by

\[ \hat{G}(\mathrm{j}\omega_m) = \frac{Y[m]}{U[m]}. \]

In the presence of noise and possible nonlinear distortions, it is often better to average multiple periods or multiple realizations, which leads to the notion of cross-spectrum based estimators:

\[ \hat{G}(\mathrm{j}\omega_m) = \frac{\hat{S}_{yu}(\omega_m)}{\hat{S}_{uu}(\omega_m)}, \]

where \( \hat{S}_{yu} \) is the estimated cross-spectrum between output and input, and \( \hat{S}_{uu} \) the auto-spectrum of the input. This kind of estimator naturally appears when we excite the system with broadband periodic signals (multi-sines, pseudo-random sequences), but the derivation still relies on the LTI assumption.

A typical experimental frequency response measurement on a robotic joint (for instance, a DC motor driving a link) follows several steps: choice of excitation frequencies, safe amplitude selection, data acquisition, and numerical estimation.

• Choice of frequencies: select a grid \( \{\omega_\ell\} \) covering the expected closed-loop bandwidth and important resonances.
• Amplitude selection: large enough to obtain a good signal-to-noise ratio (SNR), but small enough to avoid saturation and nonlinear behavior.
• Dwell time: at each frequency, allow transients to decay, then record several periods for averaging.
• Sampling: use sampling period \( T_s \) such that the Nyquist frequency \( \pi/T_s \) is well above the highest excited frequency to avoid aliasing.

The following flow diagram summarizes the procedure from experiment planning to a set of discrete frequency response samples \( \{\hat{G}(\mathrm{j}\omega_\ell)\} \).

flowchart TD
  A["Define frequency grid {omega_l}"] --> B["Check actuator limits and safety"]
  B --> C["For each omega_l: generate sinusoidal input"]
  C --> D["Wait for transients to decay"]
  D --> E["Record u(kTs), y(kTs) for several periods"]
  E --> F["Estimate complex response Ghat(j omega_l)"]
        

In robotic applications, the plant is often operated in closed loop for safety. In that case, it is possible to inject small sinusoidal perturbations on the reference or torque command and use specialized formulas (derived later in the course when feedback loops and sensitivity functions are fully introduced) to recover the open-loop or closed-loop frequency response.

We simulate a simple DC-motor-like plant in Python and then estimate its frequency response from simulated data using the DFT approach. Libraries like python-control and numpy are common in robotic control research.

import numpy as np
from control.matlab import tf, lsim  # python-control library
import matplotlib.pyplot as plt

# Continuous-time plant: approximate DC motor (position dynamics)
# G(s) = K / (tau s + 1) / s
K = 10.0
tau = 0.05
num = [K]
den = [tau, 1.0, 0.0]  # tau s^2 + s
G = tf(num, den)

# Experiment setup
Ts = 0.001   # sampling period [s]
T_total = 5.0
t = np.arange(0.0, T_total, Ts)

# Frequencies to test (rad/s)
omegas = np.array([1.0, 5.0, 10.0, 20.0])
U0 = 0.1  # input amplitude

def sine_input(omega, t, U0):
    return U0 * np.sin(omega * t)

def estimate_G_at_frequency(u, y, Ts, omega):
    """
    Estimate G(j*omega) from time data using DFT at the corresponding frequency bin.
    For simplicity we assume the record has an integer number of periods.
    """
    N = len(u)
    # Compute DFT via FFT
    U = np.fft.fft(u)
    Y = np.fft.fft(y)
    # Frequency bins
    freqs = 2.0 * np.pi * np.fft.fftfreq(N, d=Ts)  # angular frequencies
    # Find the index closest to omega
    idx = np.argmin(np.abs(freqs - omega))
    Ghat = Y[idx] / U[idx]
    return Ghat, freqs[idx]

Ghat_list = []

for omega in omegas:
    u = sine_input(omega, t, U0)
    # continuous-time simulation, then sampled
    tout, y, _ = lsim(G, U=u, T=t)
    # discard transient: keep last half of the record
    half = len(u) // 2
    u_ss = u[half:]
    y_ss = y[half:]
    Ghat, omega_eff = estimate_G_at_frequency(u_ss, y_ss, Ts, omega)
    Ghat_list.append(Ghat)
    print(f"omega target = {omega:.2f}, omega_eff = {omega_eff:.2f}, "
          f"|Ghat| = {np.abs(Ghat):.3f}, angle = {np.angle(Ghat):.3f} rad")

# Optional: compare with theoretical G(j*omega)
from control.matlab import bode

# bode returns magnitude and phase on a grid; we evaluate at omegas
mag, phase, _ = bode(G, omegas, Plot=False)
for i, omega in enumerate(omegas):
    print(f"Theoretical |G(j{omega:.2f})| = {mag[i][0]:.3f}, phase = {phase[i][0]:.3f} rad")

# Plot estimated vs theoretical magnitude
plt.figure()
plt.loglog(omegas, [np.abs(g) for g in Ghat_list], "o-", label="Estimated")
plt.loglog(omegas, [m[0] for m in mag], "x--", label="Theoretical")
plt.xlabel("omega [rad/s]")
plt.ylabel("|G(j omega)|")
plt.legend()
plt.grid(True)
plt.show()
      

In real experiments, the measured signals would come from sensors on the robot (e.g., encoder positions, motor currents) instead of a simulation. The same estimation routine can be applied to that data, provided the LTI approximation is reasonable over the chosen amplitude range.

In embedded robotic systems, it is common to estimate local frequency responses on a microcontroller or industrial PC. Libraries like Eigen can be used to perform the required linear algebra and discrete Fourier computations. Below is a minimal example computing the complex frequency response for a single test frequency using a direct least-squares fit to sinusoids.

#include <iostream>
#include <vector>
#include <complex>
#include <Eigen/Dense>

struct FRFResult {
    std::complex<double> Ghat;
    double omega;
};

FRFResult estimateFRF_singleSine(const std::vector<double> &u,
                                 const std::vector<double> &y,
                                 double Ts,
                                 double omega,
                                 double U0)
{
    std::size_t N = y.size();
    Eigen::MatrixXd Phi(N, 2);
    Eigen::VectorXd yvec(N);

    for (std::size_t k = 0; k < N; ++k) {
        double tk = k * Ts;
        Phi(k, 0) = std::cos(omega * tk);
        Phi(k, 1) = std::sin(omega * tk);
        yvec(k)   = y[k];
    }

    // Least squares: theta = (Phi^T Phi)^(-1) Phi^T y
    Eigen::Vector2d theta =
        (Phi.transpose() * Phi).ldlt().solve(Phi.transpose() * yvec);

    double Ahat = theta(0);
    double Bhat = theta(1);

    // Output phasor estimate Y(j omega) = Bhat - j Ahat
    std::complex<double> Yhat(Bhat, -Ahat);

    // Input phasor for u(t) = U0 sin(omega t) is U0 / j = -j U0
    std::complex<double> Uhat(0.0, -U0);

    std::complex<double> Ghat = Yhat / Uhat;

    FRFResult res;
    res.Ghat = Ghat;
    res.omega = omega;
    return res;
}

int main()
{
    double Ts = 0.001;
    double omega = 5.0;
    double U0 = 0.1;

    // Example data (in practice, fill from sensors)
    std::vector<double> u;
    std::vector<double> y;
    int N = 5000;
    u.reserve(N);
    y.reserve(N);

    // Here we just create synthetic data: u = U0 sin(omega t), y = 2 u (gain 2)
    for (int k = 0; k < N; ++k) {
        double tk = k * Ts;
        double uk = U0 * std::sin(omega * tk);
        double yk = 2.0 * uk; // ideal gain 2
        u.push_back(uk);
        y.push_back(yk);
    }

    // Discard first half as "transient"
    std::vector<double> u_ss(u.begin() + N / 2, u.end());
    std::vector<double> y_ss(y.begin() + N / 2, y.end());

    FRFResult r = estimateFRF_singleSine(u_ss, y_ss, Ts, omega, U0);
    std::cout << "Estimated |G(j omega)| = " << std::abs(r.Ghat)
              << ", angle = " << std::arg(r.Ghat) << " rad\n";
    return 0;
}
      

This style of implementation can be integrated into a ROS-based robotic control node, where input signals are commanded through torque or velocity interfaces and outputs are read from joint encoders.

Java is commonly used in some robotic frameworks and simulation environments. Below is a basic Java implementation that computes the FRF magnitude and phase at a single frequency using a DFT-like formula. The code assumes that u and y are populated from the data-logging mechanism of a robot.

import java.util.List;

public class FRFEstimator {

    public static class Complex {
        public double re;
        public double im;
        public Complex(double re, double im) {
            this.re = re;
            this.im = im;
        }
        public Complex add(Complex other) {
            return new Complex(this.re + other.re, this.im + other.im);
        }
        public Complex multiply(Complex other) {
            return new Complex(
                this.re * other.re - this.im * other.im,
                this.re * other.im + this.im * other.re
            );
        }
        public Complex divide(Complex other) {
            double denom = other.re * other.re + other.im * other.im;
            return new Complex(
                (this.re * other.re + this.im * other.im) / denom,
                (this.im * other.re - this.re * other.im) / denom
            );
        }
        public double abs() { return Math.hypot(re, im); }
        public double arg() { return Math.atan2(im, re); }
    }

    public static Complex estimateDFT(List<Double> x, double Ts, double omega) {
        int N = x.size();
        double re = 0.0;
        double im = 0.0;
        for (int k = 0; k < N; ++k) {
            double tk = k * Ts;
            double angle = -omega * tk;
            double c = Math.cos(angle);
            double s = Math.sin(angle);
            re += x.get(k) * c;
            im += x.get(k) * s;
        }
        return new Complex(re, im);
    }

    public static Complex estimateFRF(List<Double> u,
                                      List<Double> y,
                                      double Ts,
                                      double omega)
    {
        Complex U = estimateDFT(u, Ts, omega);
        Complex Y = estimateDFT(y, Ts, omega);
        return Y.divide(U);
    }

    public static void main(String[] args) {
        // Example usage; in practice, fill from robot logs.
        // ...
    }
}
      

This structure separates the complex arithmetic from the estimation logic and can be embedded into a larger Java-based framework for automated plant identification and controller tuning.

MATLAB and Simulink are standard tools in control and robotics. The following script simulates a plant, generates sinusoidal input, and estimates the frequency response using the FFT. In a Simulink model, the same estimation routine can be applied to logged signals from a plant model or hardware-in-the-loop setup.

% Continuous-time plant: G(s) = K / (tau s + 1)
K = 5;
tau = 0.1;
s = tf('s');
G = K / (tau*s + 1);

Ts = 0.001;
Ttotal = 5;
t = 0:Ts:Ttotal-Ts;
omegas = [1 5 10 20];
U0 = 0.2;

Ghat = zeros(size(omegas));

for idx = 1:length(omegas)
    omega = omegas(idx);
    u = U0 * sin(omega * t);
    y = lsim(G, u, t);

    % Discard transient
    N = numel(t);
    u_ss = u(N/2:end);
    y_ss = y(N/2:end);

    % FFT-based estimation
    Nss = numel(u_ss);
    U = fft(u_ss);
    Y = fft(y_ss);
    w = 2*pi*(0:Nss-1)/(Nss*Ts);

    [~, kmin] = min(abs(w - omega));
    Ghat(idx) = Y(kmin) / U(kmin);

    fprintf('omega = %.2f rad/s, |Ghat| = %.3f, angle = %.3f rad\n', ...
            omega, abs(Ghat(idx)), angle(Ghat(idx)));
end

% Compare with theoretical frequency response
[mag, phase] = bode(G, omegas);
mag = squeeze(mag);
phase = squeeze(phase);

disp('Theoretical vs estimated magnitudes:');
disp([mag(:), abs(Ghat(:))]);
      

In Simulink, one can use a Sine Wave block as excitation, Simulink-PS/Ps-Simulink blocks for interfacing with physical models, and log signals to the workspace for post-processing with scripts like the one above. Simulink Control Design and System Identification Toolbox provide higher-level blocks for automated frequency response estimation.

Wolfram Mathematica can symbolically and numerically handle transfer functions and frequency responses. The following code constructs a first-order plant and evaluates the theoretical and estimated frequency response.

(* Define transfer function model *)
K = 3;
tau = 0.2;
G = TransferFunctionModel[K/(tau*s + 1), s];

(* Frequencies of interest *)
omegas = {1, 5, 10, 20};

(* Theoretical frequency response *)
Gtheoretical = Table[
  FrequencyResponse[G, I*omega],
  {omega, omegas}
];

(* Simulate time-domain response to sinusoidal input and estimate FRF *)
Ts = 0.001;
Ttotal = 5.0;
t = Range[0, Ttotal - Ts, Ts];

U0 = 0.1;
estimateFRF[omega_] := Module[{u, y, sol, uSS, ySS, N, U, Y, freqs, idx},
  u = U0*Sin[omega*t];
  sol = OutputResponse[G, u, t];
  y = sol;

  N = Length[t];
  uSS = Take[u, {Round[N/2], N}];
  ySS = Take[y, {Round[N/2], N}];

  U = Fourier[uSS];
  Y = Fourier[ySS];

  freqs = 2*Pi*Range[0, Length[uSS] - 1]/(Length[uSS]*Ts);
  idx = First@Ordering[Abs[freqs - omega], 1];

  Y[[idx]]/U[[idx]]
];

Gest = estimateFRF /@ omegas;

Transpose@{omegas, Gtheoretical, Gest}
      

Mathematica is particularly useful for symbolic derivations, such as calculating exact expressions for \( G(\mathrm{j}\omega) \) from \( G(s) \), and then comparing them with experimentally estimated responses for validation of the LTI model.

Problem 1 (LS Estimator for Single Frequency): Let an LTI system be excited by a sinusoidal input \( u(t) = U_0 \sin(\omega t) \), and assume that the steady-state output is sampled at times \( t_k = kT_s \), \( k=0,\dots,N-1 \). The model is

\[ y[k] = A\cos(\omega kT_s) + B\sin(\omega kT_s) + v[k], \]

where \( v[k] \) is zero-mean noise. Derive the normal equations for the least-squares estimates \( \hat{A} \) and \( \hat{B} \).

Solution: Define \( \phi_1(k) = \cos(\omega kT_s) \), \( \phi_2(k) = \sin(\omega kT_s) \). The LS criterion is

\[ J(A,B) = \sum_{k=0}^{N-1}\left(y[k] - A\phi_1(k) - B\phi_2(k)\right)^2. \]

Setting partial derivatives to zero:

\[ \frac{\partial J}{\partial A} = -2\sum_{k=0}^{N-1} \phi_1(k)\left(y[k]-A\phi_1(k)-B\phi_2(k)\right) = 0, \]

\[ \frac{\partial J}{\partial B} = -2\sum_{k=0}^{N-1} \phi_2(k)\left(y[k]-A\phi_1(k)-B\phi_2(k)\right) = 0. \]

Rearranging gives the normal equations:

\[ \begin{aligned} \left(\sum_{k}\phi_1^2(k)\right)\hat{A} +\left(\sum_{k}\phi_1(k)\phi_2(k)\right)\hat{B} &= \sum_{k}\phi_1(k)y[k],\\ \left(\sum_{k}\phi_1(k)\phi_2(k)\right)\hat{A} +\left(\sum_{k}\phi_2^2(k)\right)\hat{B} &= \sum_{k}\phi_2(k)y[k]. \end{aligned} \]

In matrix form, this is exactly \( (\mathbf{\Phi}^\top\mathbf{\Phi})\hat{\boldsymbol{\theta}} = \mathbf{\Phi}^\top\mathbf{y} \).

Problem 2 (Unbiasedness of the LS Estimator): Using the same model as Problem 1, suppose \( v[k] \) is zero-mean and uncorrelated with the regressors. Show that \( \mathbb{E}[\hat{A}] = A \) and \( \mathbb{E}[\hat{B}] = B \), i.e. the LS estimator is unbiased.

Solution: In vector form, \( \mathbf{y} = \mathbf{\Phi}\boldsymbol{\theta} + \mathbf{v} \) with \( \boldsymbol{\theta} = [A\;B]^\top \). The LS estimator is

\[ \hat{\boldsymbol{\theta}} = (\mathbf{\Phi}^\top\mathbf{\Phi})^{-1}\mathbf{\Phi}^\top\mathbf{y} = (\mathbf{\Phi}^\top\mathbf{\Phi})^{-1}\mathbf{\Phi}^\top (\mathbf{\Phi}\boldsymbol{\theta}+\mathbf{v}) = \boldsymbol{\theta} + (\mathbf{\Phi}^\top\mathbf{\Phi})^{-1}\mathbf{\Phi}^\top\mathbf{v}. \]

Taking expectations and using \( \mathbb{E}[\mathbf{v}] = \mathbf{0} \) and independence from \( \mathbf{\Phi} \), we obtain

\[ \mathbb{E}[\hat{\boldsymbol{\theta}}] = \boldsymbol{\theta} + (\mathbf{\Phi}^\top\mathbf{\Phi})^{-1}\mathbf{\Phi}^\top \mathbb{E}[\mathbf{v}] = \boldsymbol{\theta}, \]

so each component is unbiased, in particular \( \mathbb{E}[\hat{A}] = A \) and \( \mathbb{E}[\hat{B}] = B \).

Problem 3 (Analytical Frequency Response of a First-Order Plant): Consider the transfer function

\[ G(s) = \frac{K}{1 + \tau s}, \]

where \( K > 0 \) and \( \tau > 0 \). Derive closed-form expressions for the magnitude and phase of \( G(\mathrm{j}\omega) \).

Solution: By substitution,

\[ G(\mathrm{j}\omega) = \frac{K}{1 + \tau\mathrm{j}\omega} = \frac{K}{1 + \mathrm{j}\tau\omega}. \]

We can write the denominator in polar form:

\[ 1 + \mathrm{j}\tau\omega = \sqrt{1 + (\tau\omega)^2}\, \mathrm{e}^{\mathrm{j}\arctan(\tau\omega)}. \]

Hence,

\[ G(\mathrm{j}\omega) = \frac{K}{\sqrt{1 + (\tau\omega)^2}}\, \mathrm{e}^{-\mathrm{j}\arctan(\tau\omega)}. \]

Therefore,

\[ |G(\mathrm{j}\omega)| = \frac{K}{\sqrt{1 + (\tau\omega)^2}},\quad \angle G(\mathrm{j}\omega) = -\arctan(\tau\omega). \]

These expressions can be directly compared with experimentally obtained \( \hat{G}(\mathrm{j}\omega) \) to validate the parametric model.

Problem 4 (Sampling and Aliasing Constraint): You want to measure experimental frequency response up to a maximum angular frequency \( \omega_{\max} \) using a sampling period \( T_s \). Derive a condition on \( T_s \) to avoid aliasing, assuming no frequency content above \( \omega_{\max} \) is excited.

Solution: The Nyquist angular frequency is \( \omega_N = \pi/T_s \). To avoid aliasing, the maximum excited frequency must satisfy \( \omega_{\max} < \omega_N \), i.e.

\[ \omega_{\max} < \frac{\pi}{T_s} \quad\Longleftrightarrow\quad T_s < \frac{\pi}{\omega_{\max}}. \]

Thus, for a desired \( \omega_{\max} \), the sampling period must satisfy \( T_s < \pi/\omega_{\max} \).

Problem 5 (Experimental Workflow for a Robotic Joint): Sketch the logical steps required to obtain experimental frequency response data of a robot joint driven by a torque-controlled motor, in which the position of the joint is measured by an encoder.

Solution (flow):

flowchart TD
  S["Start experiment"] --> F["Select safe frequency range and amplitudes"]
  F --> I["Inject sinusoidal torque commands"]
  I --> M["Measure encoder position signals"]
  M --> E["Estimate Ghat(j omega) from u,y data"]
      

The resulting \( \hat{G}(\mathrm{j}\omega) \) data can be used to refine the joint model or design a controller that meets bandwidth and robustness requirements.

In this lesson, we developed the mathematical and practical foundations of experimental frequency response measurement for SISO LTI systems. Starting from the sinusoidal steady-state behavior of LTI systems, we derived least-squares estimators for amplitude and phase at a single frequency and connected them with DFT-based formulations. We discussed experimental design issues such as frequency grids, amplitude selection, dwell times, and sampling constraints. Finally, we implemented practical estimators in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, emphasizing relevance to robotic actuators and servo drives. The resulting \( \hat{G}(\mathrm{j}\omega) \) data constitute a bridge between time-domain experiments and frequency-domain analysis that will be further exploited in subsequent chapters on Bode plots, Nyquist analysis, and frequency-domain controller design.