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; |