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 
% |   |     |          | |____________________________ data sequence 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!\n");
 else
   fprintf("H1 integrity check FAILED!!!\n");
 end

if norm(H2b-H2)<1e-12
   fprintf("H2 integrity check OK!\n");
 else
   fprintf("H2 integrity check FAILED!!!\n");
 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

References

Allemang, R.J., Brown, D.L., & Rost, R.W. (1987, December). Experimental modal analysis and dynamic component synthesis: Measurement techniques for experimental modal analysis (Vol. 2). (Report No. ADA195145). Department of Mechanical and Industrial Engineering at the University of Cincinnati, Ohio. https://apps.dtic.mil/dtic/tr/fulltext/u2/a195145.pdf.

Goldman, S. (1999). Vibration spectrum analysis: A practical approach (2nd ed.). (pp. 179–182). Industrial Press.

Greenhoe, D.J. (2019). A book concerning statistical signal processing.

Hammond, J.K. (1999). Introduction to signal processing. In: Silva J.M.M., Maia N.M.M. (Eds.), Modal Analysis and Testing: NATO Science Series (Series E: Applied Sciences) (Vol. 363, pp. 35–64). Springer, Dordrecht. https://doi.org/10.1007/978-94-011-4503-9_2.

Leclere, Q., Roozen, B., & Sandier, C. (2012). Experimental estimation of transmissibility matrices. Proceedings of ISMA2012. http://past.isma-isaac.be/downloads/isma2012/papers/isma2012_0728.pdf.

Leuridan, J., De Vis, D., Van der Auweraer, H. & Lembregts, F. (1986). A comparison of some frequency response function measurement techniques. Proceedings of the Fourth International Modal Analysis Conference. (pp. 908-918).

Otnes, R.K. & Enochson, L. (1978). Applied time series analysis: Basic techniques. (Vol. 1). John Wiley and Sons.

Rocklin, G.T., Crowley, J. & Vold, H. (1985). A comparison of h1, h2, and hv frequency response functions. Proceedings of the 3rd International Modal Analysis Conference (pp. 272-278). Union College.

Wagstaff, Peter & Henrio, J.C. (2004). Hα – a consistent estimator for frequency response functions with input and output noise. IEEE Transactions on Instrumentation and Measurement, 53(2), 457-465. 10.1109/TIM.2004.823314.

Wicks, A.L. & Vold, H. (1986). The hs frequency response function estimator. Proceedings of the 4th International Modal Analysis Conference. Union College.