# 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, ${\mathit{m}}_{1}$ and ${\mathit{m}}_{2}$, confined between two walls. The units are such that ${\mathit{m}}_{1}=1$ kg and ${\mathit{m}}_{2}=\mu $ kg. Each mass is attached to the nearest wall by a spring with elastic constant $\mathit{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 $\mathit{b}$ (in kg/s). Two sensors sample ${\mathit{r}}_{1}$ and ${\mathit{r}}_{2}$ (in $\mathrm{m}$), the displacements of the masses, at a rate ${\mathit{F}}_{\mathrm{s}}=40$ Hz.

Generate 24,000 time samples, equivalent to 10 minutes. Define the sampling interval $\Delta \mathit{t}=1/{\mathit{F}}_{\mathrm{s}}$.

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

The system can be described by the state-space model

$$\begin{array}{c}x(k+1)=Ax(k)+Bu(k),\\ y(k)=Cx(k)+Du(k),\end{array}$$

where $$x={\left[\begin{array}{cccc}{r}_{1}& {v}_{1}& {r}_{2}& {v}_{2}\end{array}\right]}^{T}$$ is the state vector, $${r}_{i}$$ and $${v}_{i}$$ are respectively the location and the velocity of the $$i$$th mass, $$u={\left[\begin{array}{cc}{u}_{1}& {u}_{2}\end{array}\right]}^{T}$$ is the vector of input driving forces, and $$y={\left[\begin{array}{cc}{r}_{1}& {r}_{2}\end{array}\right]}^{T}$$ is the output vector. The state-space matrices are

$$A=\mathrm{exp}({A}_{c}\Delta t),\phantom{\rule{1em}{0ex}}B={A}_{c}^{-1}(A-I){B}_{c},\phantom{\rule{1em}{0ex}}C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],\phantom{\rule{1em}{0ex}}D=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],$$

$\mathit{I}$ is the $4\times 4$ identity matrix, and the continuous-time state-space matrices are

$${A}_{c}=\left[\begin{array}{cccc}0& 1& 0& 0\\ -2k& -2b& k& b\\ 0& 0& 0& 1\\ k/\mu & b/\mu & -2k/\mu & -2b/\mu \end{array}\right],\phantom{\rule{1em}{0ex}}{B}_{c}=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 0\\ 0& 1/\mu \end{array}\right].$$

Set $\mathit{k}=400$ N/m, $\mathit{b}=0.1$ kg/s, and $\mu =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); end

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) set(gca,Ydir="reverse") xticks([1 2]) xticklabels("Mass m_"+[1 2]) ylabel("Time (s)") zlabel("Displacement (m)") xlim([0.5 2.5]) grid

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.

`plotComparison(fz,ztf,ft,tXY,"tfestimate")`

### 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.

`plotComparison(fz,ztf,f,frf,"modalfrf")`

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.

```
figure
modalsd(frf,f,Fs,MaxModes=10,FitMethod="lsrf")
```

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.

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

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

fn_theo = sqrt(eig([2*k -k;-k/mu 2*k/mu]))/(2*pi); disp(table(fn,fn_theo,dr,RowNames= "Mode "+(1:nModes)', ... VariableNames=["fn" "fn_theo" "dr"]))

fn fn_theo dr ______ _______ _________ Mode 1 3.8476 3.847 0.0035235 Mode 2 14.426 14.426 0.011308

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`

.

`plotComparison(fz,ztf,f,ofrf,"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`

.

### Comparison

Compare all transfer function estimation methods discussed in this example.

plotComparison(fz,ztf,ft,tXY,"tfestimate",... f,frf,"modalfrf",f,ofrf,"modalfit")

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) figure tiledlayout flow for jk = 1:2 for kj = 1:2 nexttile plot(varargin{1}, ... mag2db(abs(varargin{2}(jk,:,kj))),"--",LineWidth=2) hold on for it = 1:fix(nargin/3) plot(varargin{3*it},mag2db(abs(varargin{3*it+1}(:,jk,kj)))) end hold off grid on axis tight title("Input "+jk+", Output "+kj) xlabel("Frequency (Hz)") ylabel("Magnitude (dB)") end legend("tf definition",varargin{5:3:end},Location="EastOutside") end end