Transfer Function MATLAB Resources

  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