Hi: the three routines I am sending use zp, pressure already
converted to equivalent depth.  There is a routine around
somewhere which looks at non-hydrostatic depth attenuation in
the pressure, but I can't find it right now.  I guess hydrostatic
is ok for peak period and longer, so it may not be an issue.

Hope these are legible.

Someone else you might ask about this sort of thing is Michael Li.
He works in matlab, and does very similar sort of work.
-- 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  Douglas J. Wilson                                    +
+  Oregon Graduate Institute of Science and Technology  +
+  20000 NW Walker Rd. Beaverton, OR 97006              +
+  phone: 503 748 1099		       Fax: 503 748 1273         +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function [Usig,twoA,thetas,period] = zeroups(vx,vy,zp,F)
% This routine takes two axis velocity and a pressure height
% m/s m/s and m(H2O) and does counting-type wave statistics

dfifty = 0.00017;
vabs = sqrt(vx.*vx + vy.*vy);
vabsmean = mean(vabs);
vabs2mean = mean(vabs.*vabs);

vabsstd = std(vabs);

% CUMSUM of the velocity is a way to convert velocity to 
% excursion distance s = vdt.  

xpos  = cumsum(vx)./F;
ypos  = cumsum(vy)./F;
excur = sqrt(xpos.^2 + ypos.^2);

% ---------------------------------------------------------------------- %
%  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
%  To get a distribution of wave heights, I find the ZERO UP CROSSINGS!
zpstd  = std(zp);		 % RMS wave height.
zpmean = mean(zp);		 % mean water depth above sensor
iof = find(zp<zpmean);  % iof is the locations of all z < mean.

len_iof = length(iof);

zptest = zp;		 		 % create a fool around vector.
zpl    = length(zptest);
zptest(iof) = ones(len_iof,1).*1000;		 % make all zp<mean = 1000.
dzp_dt = zptest(2:zpl) - zptest(1:zpl-1); 		 % differentiate by 
		 		 		 		 		 		 % subtracting i+1th from ith
cr = find(dzp_dt<-900);		 		 % find all really negative slopes
% (since realzp-1000 is -999.?)  this differentiates negative slopes
% due to noise in the data from negative slopes due to zero-up crossings
		 
lcros = length(cr);		 		 % number of up crossings
%		 stf     = zeros(200,1);
%        stf(cr) = ones(lcros,1);
%		 plot([1:200],zp(1:200),'*',[1:200],stf(1:200)-1.2,'-')

% 		 Now calc. all wave heights since we have indices for all the up
% 		 crossings.  Find the max and min of zp within crossing(q) and (q+1)
%		 the p is just for multiple sections of the file (which was done
% 		 because the mean zp tends to wander, and crossings is somewhat 
% 		 sensitive to zpmean.
p = 1;
for q=1:lcros-1		 
		 pera(p) = cr(q+1) - cr(q);		 		 
		 hite(p) = max(zp(cr(q):cr(q+1)))- min(zp(cr(q):cr(q+1)));   
		 vtmax(p) = max(vabs(cr(q):cr(q+1)));
   vtmin(p) = min(vabs(cr(q):cr(q+1)));
   tvec     = detrend(excur(cr(q):cr(q+1)));
   exmax(p) = max(tvec);
   exmin(p) = min(tvec);
   twoa(p)  = exmax(p) -exmin(p);
   p = p+1;
end

% and find H1/3 by sorting hites and averaging the highest one third
ll       = length(hite); 
hsort    = sort(hite);
vmaxsort = sort(vtmax);
vminsort = sort(-vtmin);
twoasort = sort(twoa);
% persort  = sort(pera);
twothird = ll - fix((ll)/3)+1;
Honethird = mean(hsort(twothird:ll));
vmaxthird = mean(vmaxsort(twothird:ll));
vminthird = mean(vminsort(twothird:ll));
twoathird = mean(twoasort(twothird:ll));
% perathird = mean(persort(twothird:ll));
%DeltaV = vmaxthird - vminthird;
Usig   = vmaxthird;
twoA   = twoathird;
period = mean(pera)./F;
%B(2) = Usig(k);

%twoA = Usig*peraray/pi;
fw = exp(5.213*(2.5*dfifty*2./twoA).^0.194 - 5.977);
thetas = fw*Usig.^2/(32.34*dfifty);











function [peraray, twoA, Usig, alpham,alphap,alphas,thetas] = wave_stats(vx,vy,zp,F)
% Spectral processing of UVP data for wave statistics
% Input parameters needed are: VX, VY near-bottom velocity
%                              ZP  pressure converted to wave height 
%                              F   sampling frequency.
% The spectrum resolution for the wave period found by spectral peak method 
% is 512 samples.  

spctres = 512;
dfifty  = 0.00017;      %  Sand size used for Shields parameter (nondimensionalized bottom stress)

ppres = spectrum(zp,spctres);
    % [Y,II] = max(ppres(:,1));
     %peraray = spctres./(II.*F);
    % if(peraray>15)              % incident wave periods>15 sec invalid, try again.
    [Y,II] = max(ppres(9:length(ppres),1));
    peraray = spctres./(II.*F)   
     %end
         M = cov([zp(:) vx(:) vy(:)]);

         M00 = M(1,1);

         M20 = M(2,2);
         vxstd = sqrt(M20);
         vxmean = mean(vx);
         M02 = M(3,3);
         vystd = sqrt(M02);
         vymean = mean(vy);
         M01 = M(1,2);
         M10 = M(1,3);
         M11 = M(2,3);

         %vxskew = skew(vx);
         %vyskew = skew(vy);
         %pskew = skew(zp);
         
         % Wave directionality factors,  GAMLH is Longuet-Higgin's wave directionality
         
         Gam2 = (M20 +M02); 
         Gam3 = sqrt((M20-M02)^2 + 4*M11^2);
         GamLH = sqrt((Gam2 - Gam3)/(Gam2 + Gam3));

         alpham = atan(M10/M01).*(360/pi);            %  Mean wave direction
         alphap = atan(2*M11/(M20-M02)).*(180/pi);    %  Principal wave direction
         Usig = 2*sqrt(vxstd.*vxstd + vystd.*vystd);  %  Significant near-bottom wave velocity
    

         afit = polyfit(vx,vy,1);
         alphas = atan(afit(1))*(180/pi);             %  Orientation of wave orbital ellipse.
%        Now calculate the Jonsson friction factor from the orbital amplitude.
         twoA = Usig*peraray/(2.*pi); 
         fw = exp(5.213*(2.5*dfifty*2./twoA).^0.194 - 5.977);
         thetas = fw*Usig.^2/(32.34*dfifty);
 function zeta = vxtozeta(vx,F,h,z)
%  FUNCTION VXTOZETA(VX,F,H,Z)
%  This function converts current meter velocities to surface elevations
%  using linear theory, shallow water wave velocity, with all units in meters,
%  seconds and m/s.
%
%  DECLARATIONS
%		 vx = vector of velocity data to be converted.
%		 F  = sampling frequency in samples/sec.
%       h  = is water depth
%       z  = distance current meter is from the bottom.
%
%  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %  %


		 recl = length(vx);
		 spctres = 2048;
		 if(recl/spctres < 1)
		     spctres = recl;
		 end
		 vxsp   =  spectrum(vx,spctres);
		 [Y,III] = max(vxsp(:,1));
		 Peak_per = spctres./(III*F);
		  
%		  convert velocities into surface heights
		 g = 9.8;		 		 		 %  Gravity in m/s
		 sigma = 2*pi/Peak_per;		 		 %  Wave angular velocity (mean sort of)
		 kappa = sigma/sqrt(h*g);		 %  Wave number = 2pi/L (linear shallow
		 		 		 		 		 %  water wave theory for convenience.


		  zeta = sinh(kappa*h)*vx./(sigma*cosh(kappa*(h-z)));