%..........................Wave Spectral Analysis..........................
%Author: James Bramante
%Date: September 3, 2015
%
% This function performs spectral analysis on water level/depth output at
% constant frequency from models. It was written to be nearly identical
% with a similar function used to perform spectral analysis on wave
% pressure sensor data, and thus has some extraneous features and variable
% names relating more to actual pressure measurements than model output

function [hsall,hrmsall,hs,hsig,hslf,tp,tpig,tplf,tdom,hrms,hrmsig,hrmslf,tm01all,tm01,tm01ig,tm01lf,freq,EDS1] = spectral_wave_artificial(pressBurst,pressAir,fs,sensorHeight,welchSize,coord,cutoffs)


%% Defaults / Final variables
FREQ = 1; %Frequency in hertz
SENSOR_HEIGHT = 0; %Height of sensor above bed in m (included to model this function on one using pressure sensor output)
WELCH = 'all'; %Length of hanning window / frequency
g = 9.80665; %gravitational acceleration in m/s^2
RHO = 1027; %Water density in kg m^-3
TEMP = 28; %Water temperature in degrees C
CUTOFFS = {[1/3,1/25],[1/25,1/250],[1/250,1/1000]};

% Fill in non-supplied input
if nargin < 3
    fs = FREQ;
    sensorHeight = SENSOR_HEIGHT;
    welchSize = WELCH;
    cutoffs = CUTOFFS;
else if nargin < 4
        sensorHeight = SENSOR_HEIGHT;
        welchSize = WELCH;
        cutoffs = CUTOFFS;
    else if nargin < 5
            welchSize = WELCH;
            cutoffs = CUTOFFS;
        else if (nargin < 6) & (length(coord) < 3) & (~isempty(coord))
                coord(3) = TEMP;
            else if nargin < 7
                    cutoffs = CUTOFFS;
                end
            end
        end
    end
end

if isempty(welchSize)
    welchSize = WELCH;
end
if isempty(sensorHeight)
    sensorHeight = SENSOR_HEIGHT;
end
if isempty(fs)
    fs = FREQ;
end
if isempty(cutoffs)
    cutoffs = CUTOFFS;
end

%% Defining cutoffs for frequency groups (wind, IG, long wave)

wvf1=cutoffs{1}(1);       % 0.33 Hz
wvf2=cutoffs{1}(2);      % 0.04 Hz

igf1=cutoffs{2}(1);      % 0.04 Hz
igf2=cutoffs{2}(2);     % 0.004 Hz

lff1=cutoffs{3}(1);     % 0.004 Hz
lff2=cutoffs{3}(2);    % 0.001 Hz

[a, b] = size(pressBurst);

if strcmpi(welchSize,'all')
    welchSize = b / fs;
end

bnum = 1:a;

%% Average over every hour
pressHour = nanmean(pressBurst,2);


%% Pre-allocate output variables
if 2^nextpow2(b/2+1) < b
    EDS1 = NaN(a,2^nextpow2(b/2+1));
else
    EDS1 = NaN(a,max((welchSize * fs) / 2,1));
end
hs=NaN(length(bnum),1);   % short-period waves
hrms=NaN(length(bnum),1);
tp=NaN(length(bnum),1);
tm01=NaN(length(bnum),1);
Ewv=NaN(length(bnum),1);
tm01all=NaN(length(bnum),1);

hsig=NaN(length(bnum),1); % infragravity waves
hrmsig=NaN(length(bnum),1);
tpig=NaN(length(bnum),1);
Eig=NaN(length(bnum),1);
tm01ig=NaN(length(bnum),1);

hslf = NaN(length(bnum),1);  % very-low frequency waves
hrmslf = NaN(length(bnum),1);
Elf = NaN(length(bnum),1);
tplf = NaN(length(bnum),1);
tm01lf=NaN(length(bnum),1);

tdom = NaN(length(bnum),1);  % dominant period overall
hrmsall = NaN(length(bnum),1);  % Root-mean-square wave height for all freq
hsall = NaN(length(bnum),1); % Significant wave height for all frequencies
eta = NaN(size(pressBurst)); % Free surface elevation

for x = 1 : a
    pressBurst(x,:) = (pressBurst(x,:) - pressHour(x));
    
    %% First try Welch method and chop off first value
    if b >= 256
        nfft = [];
    else
        nfft = min(2^nextpow2(b/2+1),b);
    end
    if (~welchSize) 
        [EDS, freq] = pwelch(pressBurst(x,:),[],[],nfft,fs);
    else
        [EDS,freq] = pwelch(pressBurst(x,:),hanning(max(welchSize * fs,1)), [], nfft, fs);
    end
    
    EDS = EDS(2:end);
    freq = freq(2:end);

    EDS1(x,:) = EDS;
    
    %% Separate frequencies into pre-defined bins
    %All frequencies and dominant frequency
    iall=find(freq>=lff2 & freq<=wvf1);
    peaks = localpeaks_cheriton(EDS(iall));  % instead of searching for the max, I used a peak-finder code (trying to find "real" peaks, and not just the tail-end maximas)
    ipeaks=find(peaks==1);
    if ~isempty(ipeaks)
        [ps,is]=sort(EDS(iall(ipeaks)),'descend');
        ip=ipeaks(is(1));
        tdom(x) = 1./freq(iall(ip));   % dominant period over entire freq spectrum
    end
    
    %All waves
    mo = trapz(freq,EDS);
    m1 = trapz(freq,EDS.*freq);
    tm01all(x) = mo ./ m1;
    hsall(x) = 4.*sqrt(mo);
    hrmsall(x) = sqrt(8*mo);
    
    %Short-period band
    iwv = find(freq>=wvf2 & freq<wvf1);
    mo = trapz(freq(iwv),EDS(iwv));
    m1 = trapz(freq(iwv),EDS(iwv).*freq(iwv));
    tm01(x) = mo ./ m1;        % 1st moment mean wave period
    hs(x) = 4.*sqrt(mo);    % significant wave height
    hrms(x) = sqrt(8*mo);   % rms wave height
    Ewv(x) = (1/2)*(hrms(x)/2)^2;   % wave energy
    
    peaks = localpeaks_cheriton(EDS(iwv));
    ipeaks=find(peaks==1);
    if ~isempty(ipeaks)
        [ps,is]=sort(EDS(iwv(ipeaks)),'descend');
        ip=ipeaks(is(1));
        tp(x) = 1./freq(iwv(ip));  % peak period
    else 
        peaks = localpeaks_cheriton(EDS(iwv),'peaks',1);
        ipeaks = find(peaks==1);
        if ~isempty(ipeaks)
            [ps,is]=sort(EDS(iwv(ipeaks)),'descend');
            ip=ipeaks(is(1));
            tp(x) = 1./freq(iwv(ip));  % peak period
            warning(strcat('No local peaks found for sw. Frequency ',...
                'spectrum endpoints included in peak search. ',...
                'One or both endpoints identified as peaks/troughs.'));
        end
    end
    
    % infragravity band
    iig = find(freq>=igf2 & freq<igf1);
    mo = trapz(freq(iig),EDS(iig));
    m1 = trapz(freq(iig),EDS(iig).*freq(iig));
    tm01ig(x) = mo ./ m1;        % 1st moment mean wave period
    hsig(x) = 4.*sqrt(mo);    % significant wave height
    hrmsig(x) = sqrt(8*mo);   % rms wave height
    Eig(x) = (1/2)*(hrmsig(x)/2)^2;   % wave energy
    
    peaks = localpeaks_cheriton(EDS(iig));
    ipeaks=find(peaks==1);
    if ~isempty(ipeaks)
        [ps,is]=sort(EDS(iig(ipeaks)),'descend');
        ip=ipeaks(is(1));
        tpig(x) = 1./freq(iig(ip));  % peak period
    else 
        peaks = localpeaks_cheriton(EDS(iig),'peaks',1);
        ipeaks = find(peaks==1);
        if ~isempty(ipeaks)
            [ps,is]=sort(EDS(iig(ipeaks)),'descend');
            ip=ipeaks(is(1));
            tpig(x) = 1./freq(iig(ip));  % peak period
            warning(strcat('No local peaks found for ig. Frequency ',...
                'spectrum endpoints included in peak search. ',...
                'One or both endpoints identified as peaks/troughs.'));
        end
    end
    
    % very low freq band
    ilf = find(freq>=lff2 & freq<lff1);
    mo = trapz(freq(ilf),EDS(ilf));
    m1 = trapz(freq(ilf),EDS(ilf).*freq(ilf));
    tm01lf(x) = mo ./ m1;        % 1st moment mean wave period
    hslf(x) = 4.*sqrt(mo);    % significant wave height
    hrmslf(x) = sqrt(8*mo);   % rms wave height
    Elf(x) = (1/2)*(hrmslf(x)/2)^2;   % wave energy
    
    peaks = localpeaks_cheriton(EDS(ilf));
    ipeaks=find(peaks==1);
    if ~isempty(ipeaks)
        [ps,is]=sort(EDS(ilf(ipeaks)),'descend');
        ip=ipeaks(is(1));
        tplf(x) = 1./freq(ilf(ip));  % peak period
    else 
        peaks = localpeaks_cheriton(EDS(ilf),'peaks',1);
        ipeaks = find(peaks==1);
        if ~isempty(ipeaks)
            [ps,is]=sort(EDS(ilf(ipeaks)),'descend');
            ip=ipeaks(is(1));
            tplf(x) = 1./freq(ilf(ip));  % peak period
            warning(strcat('No local peaks found for lf. Frequency ',...
                'spectrum endpoints included in peak search. ',...
                'One or both endpoints identified as peaks/troughs.'));
        end
        
    end
    
    
    
end

end


