Accelerometer Wave Products MATLAB Implementation

function [hmo, hspc, dominant_period] = calculate_wave_parameters ( volts )
% John Evan's processing software from Sept 2002
%
% Modified by Jim Irish 11 October 2002
%    to reduce problems due to accelerometer and spectral noise
%    to improve the agreement between GoMOOS Hmo and DWP and NDBC's
%
% CALC_SIGNIFICANT_WAVE_HEIGHT: calculate significant wave height from first moment of spectral density
%
%  USAGE:  [hmo,hspc,dominant_period,average_period,y] = calculate_wave_parameters (volts)
%
%  INPUT:
%     volts:  Accelerometer millivolts.
%
%  OUTPUT:
%     hmo:  significant wave height
%     hspc: spectra of the vertical dispalcement
%     dominant_wave_period:
%
%  U.Maine GoMOOS Buoy C - fall 2001 deployment
%        Instrument sampled at 2 Hz for 1024 samples
%        A/D output was in millivolts.
%        Accelerometer output divided by 2 before sampled, 
%                 so converstion to g's is 1.9960 g's/v
%                 - that is 0.501 volts per g acceleration.  
%        FFT by Campbell system, 1024 points with 16 bin averages
%        With sample rate of 2 Hz, xform length is 512 seconds
%        Resolution is 0.0020 Hz before bin averaging.  
%        Bin averaging, resolution is 16 time this or 0.03125 Hz
%
   df = 0.03125; 
%
%  Bin Averaged xform has 31 points running from 0.0313 to 1 Hz 
	k=1:31;
	freq(k)=df.*k;
	fsq=freq.^2;
%
% First processing step is to get normalization factor
% to convert Campbell FFT output to acceleration spectra
%   divide by 1000² to convert sample from mv² to volts²
%   multiply by 1.996² to convert from volts² to g's² (0.501 v/g)
%   multiply by 9.8² to convert from g's² to (m/s²)²
%   divide by 2 as xform values sum to twice the record variance
%   divide by 512 points as the average sums to the variance
%   multiploy by 16 as 16 bins were averaged, this gets the sum rather than the average.  
%   the result is the acceleration spectrum in (m/sec²)² that sums to the variance.  
% 
% note that some GoMOOS buoys need an additional 2500 here!
%
	anorm=16.*1.9960784.^2.*9.8.^2./2.0./512./1000.^2 * 2500;
%   
	[r,c] = size(volts);
%   
	for i=1:r


		% first step is to select acceleration spectra and normalize 
		accspec(i,:)=anorm.*volts(i,:); 
%      
		% second step is to apply a correction for accelerometer and spectral spillage noises
		% we will not do this step as it doesn't imrove things with the current algorithm
		% used below to separate accelerometer noise and wave energy.  
      
		% third step is to integrate acceleration to displacement spectra
		%
		hspc(i,:)=1./4.^2./pi.^4.*accspec(i,:)./fsq.^2;

		%
		% Next we determine the point where accelerometer noise = wave signal
		% The index of the first relative minumum in the spectrun is the notch between 
		% the low frequncy noise and the wave energy.  Use this as the start for summing 
		% the variance and the start for seeking for the index of the maximum energy for 
		% the dominant wave period. 
		s=1;
		if hspc(i,1) == 0 
			s = 2; 
		end

		%
		% find first inflection point, or relative min in file.       
		%for k=s:31
		%	if (hspc(i,k+1)-hspc(i,k)) > 0.0, break, end       
		%end   
		ddiff = diff(hspc(i,:));
		ind = find(ddiff>0);
		if isempty(ind)
			msg = sprintf ( ' row = %d, no inflection point found', i );
			disp ( msg );
			hmo(i) = NaN;
			dominant_period(i) = NaN;
			average_period(i) = NaN;
			indexy(i) = NaN;
			continue
		else
			k = ind(1);
		end


		% k is the index we want for the summation
		%     
		% The next to final step is to sum spectra to variance and estimate Hmo
		% square root variance times 4 to get NDBC first estimate of H1/3
		% This value is really first moment of spectral density function (Hmo) which is
		% about 5-10% higher than H1/3, but is calculatable.  
		%
		hmo(i)=4.*sqrt(sum(hspc(i,k:31)));
		%
		%The final step is to determine the period of the maximum spectral ampliude as the
		% dominant wave period
		[y,I] = max(hspc(i,k:31));
		wfreq = (I+k-1) * df;
		dominant_period(i) = 1 / wfreq;
		ave_freq = sum(hspc(i,k:31) .* freq(k:31)) / sum(hspc(i,k:31)); 
		indexy(i) = I+k-1;
		% end of main processing loop
	end

	%
	y  = indexy;

hmo = hmo(:);
dominant_period = dominant_period(:);
y = y(:);
return;