Main Content

Frequency-Response Analysis of MIMO System

This example shows how to estimate the multi-input/multi-output (MIMO) transfer function for a two-body oscillator. This example starts by computing the system inputs and outputs in the time domain, emulating measurements. After applying the definition of the Z-transform to calculate the frequency-domain transfer function of the system, the example uses the functions tfestimate, modalfrf, and modalfit to estimate the transfer function from the simulated measurement data. The comparison between the outputs of these functions shows how you can retrieve the transfer function from system input/output data.

Two-Body Oscillating System

An ideal one-dimensional discrete-time oscillating system consists of two masses, m1 and m2, confined between two walls. The units are such that m1=1 kg and m2=μ kg. Each mass is attached to the nearest wall by a spring with elastic constant k (in N/m). An identical spring connects the two masses. Three dampers impede the motion of the masses by exerting on them forces proportional to speed, with damping constant b (in kg/s). Two sensors sample r1 and r2 (in m), the displacements of the masses, at a rate Fs=40 Hz.

Generate 24,000 time samples, equivalent to 10 minutes. Define the sampling interval Δt=1/Fs.

Fs = 40;
dt = 1/Fs;
N = 24000;
t = dt*(0:N-1);

The system can be described by the state-space model


where x=[r1v1r2v2]T is the state vector, ri and vi are respectively the location and the velocity of the ith mass, u=[u1u2]T is the vector of input driving forces, and y=[r1r2]T is the output vector. The state-space matrices are


I is the 4×4 identity matrix, and the continuous-time state-space matrices are


Set k=400 N/m, b=0.1 kg/s, and μ=0.1.

k = 400;
b = 0.1;
mu = 0.1;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/mu b/mu -2*k/mu -2*b/mu];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/mu];
B = Ac\(A-eye(4))*Bc;
C = [1 0 0 0;0 0 1 0];
D = zeros(2);

A random input drives the masses for the first 390 seconds and then the masses are left to rest. Set the initial conditions of location and velocity to zero. Use the state-space model to compute the time evolution of the system. Plot the displacements of the masses as a function of time.

rng default
u = randn(2,N);
u(:,t>390) = 0;

y = [0;0];
x = [0;0;0;0];
for iter = 1:N
    y(:,iter) = C*x + D*u(:,iter);
    x = A*x + B*u(:,iter);

Set the system time, input, and response data to vertical vector format.

t = t';
u = u';
y = y';

Plot the system response over time for the two masses.

plot3([1 2].*ones(size(t)),t,y)
xticks([1 2])
xticklabels("Mass m_"+[1 2])
ylabel("Time (s)")
zlabel("Displacement (m)")
xlim([0.5 2.5])

Save the system data. The examples Frequency-Response Function of Two-Body Oscillator, Modal Analysis of Two-Body Oscillator, and Transfer Function Estimation of MIMO System use this data file.

save MIMOdata

Theoretical Frequency Response

The frequency-response function of a discrete-time system is the Fourier transform of its time-domain transfer function or, equivalently, the Z-transform of the transfer function evaluated at the unit circle.

Compute the theoretical frequency-response functions. Use 2048 sampling points to calculate the discrete Fourier transform.

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

nfs = 2048;
fz = 0:Fs/nfs:Fs/2;
ztf(1,:,1) = freqz(b1(1,:),a1,fz,Fs);
ztf(1,:,2) = freqz(b1(2,:),a1,fz,Fs);
ztf(2,:,1) = freqz(b2(1,:),a2,fz,Fs);
ztf(2,:,2) = freqz(b2(2,:),a2,fz,Fs);

Transfer Function Estimate using tfestimate

Use the system data and the function tfestimate without output arguments to plot the estimate of the MIMO transfer functions. Select the "mimo" option to produce all four transfer functions. Use a periodic 5000-sample Hann window to divide the signals into segments. Specify 2500 samples of overlap between adjoining segments. Store the transfer function of the system as a function of frequency.

wind = hann(5000,"periodic");
nov = 2500;

[tXY,ft] = tfestimate(u,y,wind,nov,nfs,Fs,"mimo");

Verify that the estimate computed by tfestimate coincides with the definition.


Transfer Function Estimate Using Modal Analysis

Estimate the modal frequency response function of the system using the function modalfrf. Specify the sensor data type as measured displacements.

[frf,f] = modalfrf(u,y,Fs,wind,nov,Sensor="dis");

Plot the estimates and overlay the theoretical predictions.


Use modalsd with no output arguments to generate a stabilization diagram for modal analysis and estimate the minimum number of modes. Use a maximum of 10 modes and pick the least-squares rational function estimation method for the calculation.


The stabilization diagram shows two "+" marks for 2, 4, 5 and higher model order, indicating a stable modal fit for both frequency and damping using at least two modes.

Estimate the natural frequencies, damping ratios, and mode shapes of the system. Use the function modalfit with two modes and pick the least-squares rational function estimation method for the calculation. Obtain the reconstructed transfer functions of the system in the frequency domain.

[fn,dr,ms,ofrf] = modalfit(frf,f,Fs,2,FitMethod="lsrf");

Compare the natural frequencies to the theoretical predictions for an undamped system.

fn = 2×1


fn_theo = sqrt(eig([2*k -k;-k/mu 2*k/mu]))/(2*pi)
fn_theo = 2×1


dr = 2×1


The proximity in value of the damping ratio (dr) to zero indicates an underdamped system. The natural frequencies of the system (fn) approach the equivalent undamped natural frequency (fn_theo).

Compare the transfer-function definition with the recovered transfer function from the output of modalfit.


The natural frequencies agree between the recovered and theoretical frequency-response functions despite the differences outside the natural frequencies, correlating with the two modes calculated with modalfit.


Compare all transfer function estimation methods discussed in this example.


This example uses the functions tfestimate, modalfrf, modalsd, and modalfit to perform frequency-response analysis based on the space-state model of a mass-spring-damper oscillator MIMO system. The functions tfestimate and modalfrf estimate and plot the frequency response of the system, given the information about the system input and output in the time domain. The function modalsd helps to identify the number of modes to use for modal fit. The modal parameters calculated using the function modalfit also help retrieve the frequency response of the system.

Helper Function

The function plotComparison plots and formats a comparison between theoretical and estimated transfer functions.

function plotComparison(varargin)
tiledlayout flow
for jk = 1:2
    for kj = 1:2
        hold on
        for it = 1:fix(nargin/3)
        hold off
        grid on
        axis tight
        title("Input "+jk+", Output "+kj)
        xlabel("Frequency (Hz)")
        ylabel("Magnitude (dB)")
    legend("tf definition",varargin{5:3:end},Location="EastOutside")