%% Transfer Function Estimate
% One application of Welch's method is nonparametric system identification. Assume that H is a linear, time invariant system, and x(n) and y(n) are the input to and output of H, respectively. Then the power spectrum of x(n) is related to the CPSD of x(n) and y(n) by
% P_{yx}(\omega)=H(\omega)P_{xx}(\omega).
% An estimate of the transfer function between x(n) and y(n) is
% \hat H(\omega)={{\hat P_{yx}(\omega)}\over{\hat P_{xx}(\omega)}}.
% This method estimates both magnitude and phase information. The tfestimate function uses Welch's method to compute the CPSD and power spectrum, and then forms their quotient for the transfer function estimate. Use tfestimate the same way that you use the cpsd function.
% Generate a signal consisting of two sinusoids embedded in white Gaussian noise.

rng('default')

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

%
% Filter the signal xn with an FIR moving-average filter. Compute the actual magnitude response and the estimated response.

h = ones(1,10)/10;	            % Moving-average filter
yn = filter(h,1,xn);

[HEST,f] = tfestimate(xn,yn,256,128,256,fs);
H = freqz(h,1,f,fs);

%
% Plot the results.

subplot(2,1,1)
plot(f,abs(H))
title('Actual Transfer Function Magnitude')
yl = ylim;
grid
subplot(2,1,2)
plot(f,abs(HEST))
title('Transfer Function Magnitude Estimate')
xlabel('Frequency (Hz)')
ylim(yl)
grid

% Copyright 2012 The MathWorks, Inc.