Estimate the frequency-response function for a simple single-input/single-output system and compare it to the definition.

A one-dimensional discrete-time oscillating system consists of a unit mass, $\mathit{m}$, attached to a wall by a spring with elastic constant $\mathit{k}=1$. A sensor samples the displacement of the mass at ${\mathit{F}}_{\mathrm{s}}=1$ Hz. A damper impedes the motion of the mass by exerting on it a force proportional to speed, with damping constant $\mathit{b}=0.01$.

Generate 3000 time samples. Define the sampling interval $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$.

Fs = 1;
dt = 1/Fs;
N = 3000;
t = dt*(0:N-1);
b = 0.01;

The system can be described by the state-space model

where $\mathit{x}={\left[\begin{array}{cc}\mathit{r}& \mathit{v}\end{array}\right]}^{\mathit{T}}$ is the state vector, $\mathit{r}$ and $\mathit{v}$ are respectively the displacement and velocity of the mass, $\mathit{u}$ is the driving force, and $\mathit{y}=\mathit{r}$ is the measured output. The state-space matrices are

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

Ac = [0 1;-1 -b];
A = expm(Ac*dt);

Bc = [0;1];
B = Ac\(A-eye(2))*Bc;

C = [1 0];
D = 0;

The mass is driven by random input for the first 2000 seconds and then left to return to rest. Use the state-space model to compute the time evolution of the system starting from an all-zero initial state. Plot the displacement of the mass as a function of time.

rng default
u = randn(1,N)/2;
u(2001:end) = 0;

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

plot(t,y)

Estimate the modal frequency-response function of the system. Use a Hann window half as long as the measured signals. Specify that the output is the displacement of the mass.

wind = hann(N/2);

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

The frequency-response function of a discrete-time system can be expressed as the Z-transform of the time-domain transfer function of the system, evaluated at the unit circle. Compare the modalfrf estimate with the definition.

[b,a] = ss2tf(A,B,C,D);

nfs = 2048;
fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);
ztf = polyval(b,z)./polyval(a,z);

plot(f,20*log10(abs(frf)))
hold on
plot(fz*Fs,20*log10(abs(ztf)))
hold off
grid
ylim([-60 40])

Estimate the natural frequency and the damping ratio for the vibration mode.

[fn,dr] = modalfit(frf,f,Fs,1,'FitMethod','PP')

fn = 0.1593

dr = 0.0043

Compare the natural frequency to $1/2\pi$, which is the theoretical value for the undamped system.

theo = 1/(2*pi)

theo = 0.1592

Compute the modal parameters of a Space Station module starting from its frequency-response function (FRF) array.

Load a structure containing the three-input/three-output FRF array. The system is sampled at 320 Hz.

load modaldata SpaceStationFRF

frf = SpaceStationFRF.FRF; 
f = SpaceStationFRF.f; 
fs = SpaceStationFRF.Fs;

Extract the modal parameters of the lowest 24 modes using the least-squares rational function method.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,24,'FitMethod','lsrf');

Compare the reconstructed FRF array to the measured one.

for ij = 1:3
 for ji = 1:3
    subplot(3,3,3*(ij-1)+ji)
    loglog(f,abs(frf(:,ji,ij)))
    hold on
    loglog(f,abs(ofrf(:,ji,ij)))
    hold off
    axis tight
    title(sprintf('In%d -> Out%d',ij,ji))
    if ij==3
        xlabel('Frequency (Hz)')
    end
 end
end

Estimate the frequency-response function and modal parameters of a simple multi-input/multi-output system.

An ideal one-dimensional 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$ and ${\mathit{m}}_2=\mu$. Each mass is attached to the nearest wall by a spring with an elastic constant $\mathit{k}$. 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}$. Sensors sample ${\mathit{r}}_1$ and ${\mathit{r}}_2$, the displacements of the masses, at ${\mathit{F}}_{\mathrm{s}}=50$ Hz.

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

Fs = 50;
dt = 1/Fs;
N = 30000;
t = dt*(0:N-1);

The system can be described by the state-space model

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

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

Set $\mathit{k}=400$, $\mathit{b}=0.1$, and $\mu =1/10$.

k = 400;
b = 0.1;
m = 1/10;

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

The masses are driven by random input throughout the measurement. Use the state-space model to compute the time evolution of the system starting from an all-zero initial state.

rng default
u = randn(2,N);

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

Use the input and output data to estimate the transfer function of the system as a function of frequency. Use a 15000-sample Hann window with 9000 samples of overlap between adjoining segments. Specify that the measured outputs are displacements.

wind = hann(15000);
nove = 9000;
[FRF,f] = modalfrf(u',y',Fs,wind,nove,'Sensor','dis');

Compute the theoretical transfer function as the Z-transform of the time-domain transfer function, evaluated at the unit circle.

nfs = 2048;
fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);

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

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);
frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

Plot the estimates and overlay the theoretical predictions.

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(f,20*log10(abs(FRF(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        axis([0 Fs/2 -100 0])
        title(sprintf('Input %d, Output %d',jk,kj))
    end
end

Plot the estimates by using the syntax of modalfrf with no output arguments.

figure
modalfrf(u',y',Fs,wind,nove,'Sensor','dis')

Estimate the natural frequencies, damping ratios, and mode shapes of the system. Use the peak-picking method for the calculation.

[fn,dr,ms] = modalfit(FRF,f,Fs,2,'FitMethod','pp');
fn

fn = 
fn(:,:,1) =

    3.8466    3.8466
    3.8495    3.8495


fn(:,:,2) =

    3.8492    3.8490
    3.8552   14.4684

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

undamped = sqrt(eig([2*k -k;-k/m 2*k/m]))/2/pi

undamped = 2×1

    3.8470
   14.4259

Compute the natural frequencies, the damping ratios, and the mode shapes for a two-input/three-output system excited by several bursts of random noise. Each burst lasts for 1 second, and there are 2 seconds between the end of each burst and the start of the next. The data are sampled at 4 kHz.

Load the data file. Plot the input signals and the output signals.

load modaldata

subplot(2,1,1)
plot(Xburst)
title('Input Signals')
subplot(2,1,2)
plot(Yburst)
title('Output Signals')

Compute the frequency-response functions. Specify a rectangular window with length equal to the burst period and no overlap between adjoining segments.

burstLen = 12000;
[frf,f] = modalfrf(Xburst,Yburst,fs,burstLen);

Visualize a stabilization diagram and return the stable natural frequencies. Specify a maximum model order of 30 modes.

figure
modalsd(frf,f,fs,'MaxModes',30);

Zoom in on the plot. The averaged response function has maxima at 373 Hz, 852 Hz, and 1371 Hz, which correspond to the physical frequencies of the system. Save the maxima to a variable.

phfr = [373 852 1371];

Compute the modal parameters using the least-squares complex exponential (LSCE) algorithm. Specify a model order of 6 modes and specify physical frequencies for the 3 modes determined from the stabilization diagram. The function generates one set of natural frequencies and damping ratios for each input reference.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,6,'PhysFreq',phfr);

Plot the reconstructed frequency-response functions and compare them to the original ones.

for k = 1:2
    for m = 1:3
        subplot(2,3,m+3*(k-1))
        plot(f/1000,10*log10(abs(frf(:,m,k))))
        hold on
        plot(f/1000,10*log10(abs(ofrf(:,m,k))))
        hold off
        text(1,-50,[['Output ';' Input '] num2str([m k]')])
        ylim([-100 -40])
    end
end
subplot(2,3,2)
title('Frequency-Response Functions')