# Transfer Function MATLAB Resources

March 22, 2019

MATLAB® offers operators supporting several statistical estimators using Welch’s method. These include the following:

 tfestimate(…) Estimate H1 and H2. https://www.mathworks.com/help/signal/ref/tfestimate.html pwelch(…) Estimate power spectral density. https://www.mathworks.com/help/signal/ref/pwelch.html cpsd(…) Estimate cross power spectral density. https://www.mathworks.com/help/signal/ref/cpsd.html

### Example Code from MATLAB

```  clear                    % clear all variables
format long g;           % set format of text output
N       = 16;            % length of input data sequence x
Fs      = 1;             % sample rate (used for scaling)
percent = 50;            % percent overlap
Nfft = N * (1-percent/100);
x= [0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0]';
y= [0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0]';
%   0     3       7       11  13  15
%----------------------------------------------------------------------------
% smoothing window
%----------------------------------------------------------------------------
window = hanning(Nfft);  % Hanning     window
%----------------------------------------------------------------------------
% calculate estimates
%----------------------------------------------------------------------------
[H1, freq]=tfestimate(x,y,window,Noverlap,Nfft,Fs,'Estimator','H1');
[H2, freq]=tfestimate(x,y,window,Noverlap,Nfft,Fs,'Estimator','H2');
[Sxy,freq]=cpsd(      x,y,window,Noverlap,Nfft,Fs);  % estimate Sxy
[Sxx,freq]=pwelch(    x,  window,Noverlap,Nfft,Fs);  % estimate Sxx
[Syy,freq]=pwelch(      y,window,Noverlap,Nfft,Fs);  % estimate Syy
% |   |     |         | |  |      |        |    |_____ sample rate
% |   |     |         | |  |      |        |__________ estimate length
% |   |     |         | |  |      |___________________ num. overlapping positions
% |   |     |         | |  |__________________________ window type
% |   |     |         | |_____________________________ datasequence y
% |   |     |         |_______________________________ data sequence x
% |   |     |_________________________________________ estimation operation
% |   |_______________________________________________ x-axis sequence
% |___________________________________________________ estimate sequence Syx = conj(Sxy);
% estimate Syx H1b = Syx ./ Sxx;
% estimate alternate H1
H2b = Syy ./ Sxy;
% estimate alternate H2
%----------------------------------------------------------------------------
% estimate scaling parameter s
%----------------------------------------------------------------------------
% s = (norm(y)-(mean(y)*mean(y)))/(norm(x)-(mean(x)*mean(x)))
s = var(y)/var(x)
%----------------------------------------------------------------------------
% estimate Hs and Hv
%----------------------------------------------------------------------------
Hs = estimate_Hs(s,Sxx,Syy,Sxy);
% estimate Hs
Hv = estimate_Hv(Sxx,Syy,Sxy);
% estimate Hv
%----------------------------------------------------------------------------
% integrity checks
%----------------------------------------------------------------------------
if norm(H1b-H1)<1e-12
fprintf("H1 integrity check OK!
"); else
fprintf("H1 integrity check FAILED!!!");
end
if norm(H2b-H2)<1e-12
fprintf("H2 integrity check OK!
"); else
fprintf("H2 integrity check FAILED!!!");
end
```
```function Hv = estimate_Hv(Sxx,Syy,Sxy)
s=1;
Hv = estimate_Hs(s,Sxx,Syy,Sxy);
end

function Hs = estimate_Hs(s,Sxx,Syy,Sxy)
u1 = Syy;u2 = s * Sxx;
u3 = (s * Sxx - Syy) .* (s * Sxx - Syy);
u4 = abs(Sxy) .* abs(Sxy);
v1 = 2.0 * Sxy;
Hs = (u1 - u2 + sqrt(u3 + 4 * s * u4)) ./ v1;
end
```