Compute and plot the transfer function estimate between two sequences, x and y. The sequence x consists of white Gaussian noise. y results from filtering x with a 30th-order lowpass filter with normalized cutoff frequency $$0.2\pi$$ rad/sample. Use a rectangular window to design the filter. Specify a sample rate of 500 Hz and a Hamming window of length 1024 for the transfer function estimate.

h = fir1(30,0.2,rectwin(31));
x = randn(16384,1);
y = filter(h,1,x);

fs = 500;
tfestimate(x,y,1024,[],[],fs)

Use fvtool to verify that the transfer function approximates the frequency response of the filter.

fvtool(h,1,'Fs',fs)

Obtain the same result by returning the transfer function estimate in a variable and plotting its absolute value in decibels.

[Txy,f] = tfestimate(x,y,1024,[],[],fs);

plot(f,mag2db(abs(Txy)))

Estimate the transfer 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 of unit elastic constant. A sensor samples the acceleration, $\mathit{a}$, of the mass at ${\mathit{F}}_{\mathit{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 2000 time samples. Define the sampling interval $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$.

Fs = 1;
dt = 1/Fs;
N = 2000;
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 position and velocity of the mass, $\mathit{u}$ is the driving force, and $\mathit{y}=\mathit{a}$ 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(size(A)))*Bc;

C = [-1 -b];
D = 1;

The mass is driven by random input for half of the measurement interval. Use the state-space model to compute the time evolution of the system starting from an all-zero initial state. Plot the acceleration of the mass as a function of time.

rng default

u = zeros(1,N);
u(1:N/2) = randn(1,N/2);

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 transfer function of the system as a function of frequency. Use 2048 DFT points and specify a Kaiser window with a shape factor of 15. Use the default value of overlap between adjoining segments.

nfs = 2048;
wind = kaiser(N,15);

[txy,ft] = tfestimate(u,y,wind,[],nfs,Fs);

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. Verify that the estimate computed by tfestimate coincides with this definition.

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

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

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

Plot the estimate using the built-in functionality of tfestimate.

tfestimate(u,y,wind,[],nfs,Fs)

Estimate the transfer function for 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{a}}_1$ and ${\mathit{a}}_2$, the accelerations of the masses, at ${\mathit{F}}_{\mathit{s}}=50$ Hz.

Generate 30000 time samples, equivalent to 600 seconds. Define the sampling interval $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{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}{a}_1& a_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$, and $\mu =1/10$.

k = 400;
b = 0;
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 = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];
D = [1 0;0 1/m];

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

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. Specify the 'mimo' option to produce all four transfer functions. Use a 5000-sample Hann window to divide the signals into segments. Specify 2500 samples of overlap between adjoining segments and $2^{14}$ DFT points. Plot the estimates.

wind = hann(5000);
nov = 2500;

[q,fq] = tfestimate(u',y',wind,nov,2^14,Fs,'mimo');

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

nfs = 2^14;

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 theoretical transfer functions and their corresponding estimates.

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(fq,20*log10(abs(q(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        grid
        title(['Input ' int2str(kj) ', Output ' int2str(jk)])
        axis([0 Fs/2 -50 100])
    end
end

The transfer functions have maxima at the expected values, ${\omega }_{1,2}/2\pi$, where the $\omega$ are the eigenvalues of the modal matrix.

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

ans = 2×1

    3.8470
   14.4259

Add damping to the system by setting $\mathit{b}=0.1$. Compute the time evolution of the damped system with the same driving forces. Compute the ${\mathit{H}}_2$ estimate of the MIMO transfer function using the same window and overlap. Plot the estimates using the tfestimate functionality.

b = 0.1;

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);
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];

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

clf
tfestimate(u',y',wind,nov,[],Fs,'mimo','Estimator','H2')
legend('I1, O1','I1, O2','I2, O1','I2, O2')

yl = ylim;

Compare the estimates to the theoretical predictions.

[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(fz*Fs,20*log10(abs(reshape(permute(frf,[2 1 3]),[nfs/2 4]))))
legend('I1, O1','I1, O2','I2, O1','I2, O2')
ylim(yl)
grid