%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WAVES 1 - LAB 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PURPOSE : compute the frequency surface elevation spectrum from surface
% elevation time series. Uses Welch method.
% HISTORY : - origination in matlab
% - 01.2019 comments and rewriting, C. Peureux, cpeureux@ifremer.fr
% REQUIREMENTS : katsiveli_data_sample.mat

clear all ;
%close all ;
clc ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. Load surface elevation time series (wave gauge)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load katsiveli_data_sample   % reads previous extracted time series
elevation=squeeze(elevation1);
nmes=length(elevation); % Number of measurements. Here, full time series are taken

%%% Plot
figure(1) % first plot : time series 
clf
% set(gca,'FontSize',16,'LineWidth',1);
plot(times(1:nmes),elevation(1:nmes),'k-','LineWidth',2);
xlabel('time (s)');
ylabel('surface elevation (m)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2. Process time series into frequency spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Smoothing
%%% Below is 3-point filter with weights [ 0.25 0.5 0.25 ] -> data smoothing
%%% (empirical correction sometimes required)
elevation=elevation1;
elevation(2:nmes-1)=0.25.*((elevation(3:nmes)+elevation(1:nmes-2))+2.*elevation(2:nmes-1));

%%%% Fourier vectors calculation
%%% parameters
nfft=1000; % Fourier transform window length (number of samples) -> preferably a power of 2 (fft algorithm more efficient)
overlap = 50 ; % overlap between windows [%]
%%% constants
Fs=1./(times(2)-times(1)); % sampling frequency [Hz]
df=Fs/nfft; % Frequency resolution [Hz]
numoverlap = overlap*nfft/100 ; % overlap between windows [number of samples]
Nf=nfft/2+1; % number of positive frequencies
Nfo=(Nf-1)/2;
freq1=linspace(0,df*(nfft-1),nfft);
freq=linspace(df,df*(Nf-1),Nf-1)';
ov = overlap/100. ; % More convenient

%%%% Spectrum calculation
%%% Windowing
hanning=transpose(0.5 * (1-cos(2*pi*linspace(0,nfft-1,nfft)/(nfft-1)))); % Hanning window in time
% hanning=hanning.*0+1; % decomment this line to remove window
wc2=1/mean(hanning.^2); % Window correction factor (will be useful later)  
NS1=floor(nmes/nfft); % Number of windows without overlap
NS=NS1*(1+2*ov)-2*ov; % Number of windows with overlap -> the more overlap, the more windows
Eh=ones(1,NS);
H=hanning*Eh;% array with Hanning windows on each column

%%%% Data preconditioning
%%% All spectra will be computed at the same time, for each column of this
%%% elevation matrix :
elevmat=zeros(nfft,NS);
%%% Now fills array with the contents of each window (Welch method)
%%% Loop over windows
for iw=1:NS
    nstart = floor(1+(iw-1)*(1-ov)*nfft) ; % floor : integer part
    nend = nstart+nfft-1 ;
    elevmat(:,iw) = elevation(nstart:nend) ;
end
elevmat2=elevmat ; %detrend(elevmat,'linear'); % We remove the linear trend from each column
Hs1=4.*sqrt(mean(var(elevmat2))) % We compute and display the significant wave height from each time serie
 
elevmat2=H.*detrend(elevmat2,'constant'); % remove mean and apply window

Hs2=4.*sqrt(mean(var(elevmat2))) % Check the impact of windowing and linear detrend on the wave height ... windwowing must be corrected
   
%%%% Computes FFT
zspec=fft(elevmat2,nfft,1)/nfft; % Complex Fourier spectral amplitudes [m].%%% Repeated for all columns. nfft = normalization factor.
%%% computes Power Spectral Density (PSD) 
Szz=(abs(zspec).^2)*wc2/df; % [m^2/Hz]
%%% Folding spectrum onto positive frequencies only
Ef=2.*mean(Szz(2:Nf,:),2); % = single-sided spectrum [m^2/Hz]
%%% Double check that variance is conserved 
Hs3=4.*sqrt(sum(Ef.*df)) % must be close to Hs1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 3. Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7) % first plot : amplitudes
clf
mpdc10 = hsv(NS);
ha = axes; hold(ha,'on')
% set(ha,'ColorOrder',mpdc10) 
plot(abs(zspec),'-','LineWidth',1);
set(gca,'FontSize',16,'LineWidth',1);
xlabel('frequency index');
ylabel('|Z| (m)');

figure(8) % second : phases
clf
ha = axes; hold(ha,'on')
% set(ha,'ColorOrder',mpdc10) 
plot(angle(zspec)./pi,'.','LineWidth',1);
set(gca,'FontSize',16,'LineWidth',1);
xlabel('frequency index');
ylabel('phase / \pi');
xlim([1,nfft])

%%%% Now a few plots %%%% 
figure(10)
clf
loglog(freq,Ef,'k-','LineWidth',2);
set(gca,'FontSize',16,'LineWidth',1);
xlabel('frequency (Hz)');
ylabel('E(f)  ({m^2/Hz})');
axis([0.01 8 1E-8 1.])
grid on
hold on ;
lwb = chi2inv(0.025,2*NS1)*Ef/(2*NS1) ; % 95% interval lower bound
upb = chi2inv(0.975,2*NS1)*Ef/(2*NS1) ; % 95% interval upper bound
loglog(freq,lwb,'r-','LineWidth',0.5)
loglog(freq,upb,'r-','LineWidth',0.5)

%%%% plots spectra before averaging %%%% 
figure(30)
clf
loglog(freq,2.*Szz(2:Nf,:),'LineWidth',2);
% set(gca,'FontSize',16,'LineWidth',1);
xlabel('frequency (Hz)');
ylabel('E(f)  ({m^2/Hz})');
hold on
plot(freq,Ef,'k-','LineWidth',3)
axis([0.01 8 1E-7 0.1])
grid on

%%%% plots non-dimensional spectrum %%%% 
figure(20)
clf
fac=(2*pi)^4/(9.81^2);
plot(freq,Ef.*freq.^5.*fac,'LineWidth',2);

% set(gca,'FontSize',14);
xlabel('frequency (Hz)');
ylabel('E(f) {\times \sigma^5/(2 \pi g^2)}  (no dimensions)')

%%%%%%%%%%%%%% ANSWERS TO QUESTIONS
I = freq>1.5 ; % boolean
beta = trapz(freq(I),Ef(I))/trapz(freq,Ef)

Ff = (2*pi*freq).^4.*Ef/9.8^2 ;
gam = trapz(freq(I),Ff(I))/trapz(freq,Ff)