function [mag1sig,phs1sig,coh1sig,cps1sig,real1sig,imag1sig] = genstat(G,Gxy,Coh,navg)
%
%  Name:   genstat
%
%  Usage:	[mag1sig,phs1sig,coh1sig,cps1sig,real1sig,imag1sig] = genstat(G,Gxy,Coh,navg)
%
%  Purpose:  Generate the statistical (1 sigma) error bounds on the
%	spectral estimates from a modal test:
%
%	Inputs:
%		G = measured FRF estimate
%		Gxy = measured CPS estimate
%		Coh = measured coherence estimate
%		navg = number of experimental averages
%
%	Outputs:
%		mag1sig = 1 sigma std for FRF magnitude
%		phs1sig = 1 sigma std for FRF phase
%		coh1sig = 1 sigma std for coherence function
%		cps1sig = 1 sigma std for CPS magnitude
%

%  Version SWD970423
%
%  Author(s): Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%               ESA-EA
%
%  Copyright (c) 1996-1997, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)
%

%
%  Interpolate points in Coherence where the value is zero
%

[zpveci,zpvecj] = find(Coh == 0);

for i = 1:length(zpveci),

   Coh(zpveci(i),zpvecj(i)) = (Coh(zpveci(i)+1,zpvecj(i)) + Coh(zpveci(i)-1,zpvecj(i)))/2;

end


%
%   FRF magnitude
%

mag1sig = sqrt(1 - Coh) ./ sqrt(Coh * 2 * navg) .* abs(G);

%
%   FRF phase
%

%phs1sig = 0.5 * sqrt(1 - Coh) ./ sqrt(Coh * 2 * navg) .* angle(G);

ephs = sqrt(1 - Coh) ./ sqrt(Coh * 2 * navg);

%
%   Make sure argument for asin is always <=1

gt1 = find(ephs>1);

ephs(gt1)=ones(length(gt1),1);

phs1sig = asin(ephs);

%
%   Coherence
%

coh1sig = sqrt(2) * (1 - Coh) ./ sqrt(Coh * navg) .* Coh;

%
%   Cross-Power Spectrum
%

cps1sig = 1 ./ sqrt(Coh * navg) .* Gxy;

%
%   Compute bounds on FRF based on mag and phase bounds
%

gmag = abs(G);
gph = angle(G);
sgph = sin(gph);
cgph = cos(gph);

real1sig = sqrt( (cgph .* mag1sig).^2 + (gmag .* sgph .* phs1sig).^2 ); 

imag1sig = sqrt( (sgph .* mag1sig).^2 + (gmag .* cgph .* phs1sig).^2 ); 

return