function [Hs,Tz,Te,Tpc]=Wavestatsfromspectra(S,f,df);

% define basic spectral moments
m0 = sum(S.*df);  %variance
m1 = sum(f.*S.*df); %mean
m2 = sum((f.^2).*S.*df);
m3 = sum((f.^3).*S.*df);
m4 = sum((f.^4).*S.*df);
mn1 = sum((f.^-1).*S.*df);
mn2 = sum((f.^-2).*S.*df);

% calculate mean spectral stats
Hs = 4*sqrt(m0);        % significant wave height
Tz = sqrt(m0/m2);       % mean zero crossing period
Te = mn1/m0;            % Energy period
Tpc = (mn2.*m1)./(m0^2);  % calculated peak period (an approx)(noisy)