tfestimate

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)

Verify that the transfer function approximates the frequency response of the filter.

freqz(h,1,[],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}$ (in kg), 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$ kg/s.

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,mag2db(abs(txy)))
hold on
plot(fz,mag2db(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 multi-input/multi-output (MIMO) system.

Two masses connected to a spring and a damper on each side form an ideal one-dimensional discrete-time oscillating system. The system input array u consists of random driving forces applied to the masses. The system output array y contains the observed displacements of the masses from their initial reference positions. The system is sampled at a rate Fs of 40 Hz.

Load the data file containing the MIMO system inputs, the system outputs, and the sample rate. The example Frequency-Response Analysis of MIMO System analyzes the system that generated the data used in this example.

load MIMOdata

Estimate and plot the frequency-domain transfer functions of the system using the system data and the function tfestimate. Select the "mimo" option to produce all four transfer functions. Use $2^{14}$ sampling points to calculate the discrete Fourier transform, divide the signal into 5000-sample segments, and window each segment with a Hann window. Specify 2500 samples of overlap between adjoining segments.

wind = hann(5000);
nfs = 2^14;
nov = 2500;

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

tiledlayout flow
for jk = 1:2
    for kj = 1:2
        nexttile
        plot(ft,mag2db(abs(tXY(:,jk,kj))))
        grid on
        ylim([-120 0])
        title("Input "+jk+", Output "+kj)
        xlabel("Frequency (Hz)")
        ylabel("Magnitude (dB)")
    end
end