%% Frequency-Response Function of SISO System
% 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, m, attached to a wall by a spring with elastic constant k=1. A sensor samples the displacement of the mass at 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 b=0\ldotp 01.
%
%
% Generate 3000 time samples. Define the sampling interval \Delta t=1/F_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
% \matrix{x(k+1)=Ax(k)+Bu(k),\cr y(k)=Cx(k)+Du(k),\cr}
% where x={\left\lbrack \begin{array}{cc}
% r & v
% \end{array}\right\rbrack }^T is the state vector, r and v are respectively the displacement and velocity of the mass, u is the driving force, and y=r is the measured output. The state-space matrices are
% A=\exp(A_{\rm c}\Delta t),\quad B=A_{\rm c}^{-1}(A-I)B_{\rm c},\quad
% C=\left[\matrix{1&0}\right],\quad D=0,
% I is the 2\times 2 identity, and the continuous-time state-space matrices are
% A_{\rm c}=\left[\matrix{0&1\cr-1&-b\cr}\right],\quad B_{\rm c}=\left[\matrix{0\cr1\cr}\right].

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')

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


theo = 1/(2*pi)

%
%
% Copyright 2012 The MathWorks, Inc.