clear all
close all

addpath(genpath('/Users/Mau/Documents/MATLABhelp_new/'));

% Magdalena Carranza
% Testing pwelch function against my code to compute spectra

% path_mat='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-11-16/matfiles/';
path_mat='/Users/Mau/Documents/MBARIprojects/CPF-MR/ParagonTests/2021-11-16/matfiles/'; % updated file!
path_plots='/Users/Mau/Documents/MBARIprojects/MR/TurbulenceData/ParagonTests/2021-11-16/plots/';

load([path_mat,'w_CTD_MR_allprofiles.mat']) % has profile data for CTD (1Hz and 2m-binned) and MR data
% load([path_mat,'w_CTD_MR_allprofiles_v2.mat']) % saved version 2 in dpdt_ctdmr_comparion.m (just in case I modified fs)

% Use pwelch to compare W from the 3 signals using default values

%% Define some color for plotting
cbar=cbrewer('seq','Blues',101);
% cbar=cbar(5:5:end,:);
lblue=cbar(20,:);
blue=cbar(60,:);

Reds=cbrewer('seq','Reds',101);
lred=Reds(20,:);
red=Reds(50,:);

Greens=cbrewer('seq','Greens',101);
lgreen=Greens(20,:);
green1=Greens(40,:);
green2=Greens(60,:);

lgrey=[0.7 0.7 0.7];

% Get default color order
for p=1:8
h = plot(1:10,1:10); hold on
Colors = get(h,'Color');
end

%% Truncate data, to use 10-80m
zi=10; zf=80;

di=searchclosest(P_1hz(:,1),zi); % 1Hz data
df=searchclosest(P_1hz(:,1),zf);
pres1=P_1hz(di:df,:);
data1=W_1hz(di:df,:);

di=searchclosest(D_2mb(:,1),zi); % 2-m binned (1Hz data)
df=searchclosest(D_2mb(:,1),zf);
pres2=D_2mb(di:df,:);
data2=W_2mb(di:df,:);
% Ndata2=N_interp(di:df,:);

di=searchclosest(mr_paxis(:,1),zi); % MR 64 Hz)
df=searchclosest(mr_paxis(:,1),zf);
pres3=mr_paxis(di:df,:);
data3=W_mr(di:df,:);

fs1=1; % 1 Hz;
fs2=1/10; % for mean dt ~10s (when binning data int0 2-m bins, and profiling at 0.15 m/s)
fs3=1./(0.3/0.15); % from MR, sampled at 64Hz, but 1 sample every 0.3dbar profiling at 0.15 dbar/s

dz=nanmean(diff(mr_paxis(di:df,1))); % better estimate for MR data, but still close to 0.49 (note butterworth filter cuttoff is 0.5)
PV=nanmean(data3(:));
fs3=1./(dz/PV);

input='raw';
method='raw';
% NOTE: in this version, plotting spectra for RAW data each with its of
% sampling frequency

%% PLot raw PROFILE data (for a given mission)

close all
plot(W_1hz,P_1hz); 
axis ij
ylabel('Depth');
xlabel('w_{1hz} [m/s]');

xlim([0 0.24])
ylim([0 100])
print_setup(20,40)
print([path_plots,'Profiles_w1hz_',input,'_',method,'.jpg'],'-djpeg');

close all
plot(W_2mb,D_2mb); 
axis ij
ylabel('Depth');
xlabel('w_{2mb} [m/s]');

ylim([0 100])
xlim([0 0.24])
print_setup(20,40)
print([path_plots,'Profiles_w2mb_',input,'_',method,'.jpg'],'-djpeg');

close all
plot(W_mr,mr_paxis); 
axis ij
ylabel('Depth');
xlabel('w_{mr} [m/s]');

xlim([0 0.24])
ylim([0 100])
print_setup(20,40)
print([path_plots,'Profiles_wmr_',input,'_',method,'.jpg'],'-djpeg');


%%  Truncate data, to use 10-80m bin, but over interpolated data instead
zi=10; zf=80;

di=searchclosest(p1hz_interp(:,1),zi); % 1Hz data
df=searchclosest(p1hz_interp(:,1),zf);
pres1=p1hz_interp(di:df,:);
data1=W_1hz_interp(di:df,:);

di=searchclosest(p2mb_interp(:,1),zi); % 2-m binned (1Hz data)
df=searchclosest(p2mb_interp(:,1),zf);
pres2=p2mb_interp(di:df,:);
data2=W_2mb_interp(di:df,:);

di=searchclosest(pmr_interp(:,1),zi); % MR 64 Hz)
df=searchclosest(pmr_interp(:,1),zf);
pres3=pmr_interp(di:df,:);
data3=W_mr_interp(di:df,:);

dz=nanmean(diff(pres3)); PV=nanmean(data3(:));
dtt=dz./PV;
fs3=1./dtt;
fs1=fs3; fs2=fs3; % ~69 Hz %% NOTE!!! all timeseries interpolated to the same grid. thus, fs is the same (~for the highes freq s)

fs3=64;
fs1=fs3; fs2=fs3;

input='interp';
method='raw';

%% Processing w signal (as for glider data)

fs=fs3; % sampling rate; 1/2.6;?
% fc_lowpass=0.1; % 0.01 (for glider data)

fc_lowpass=0.004; % 4x10(-2) where spectral energy starts to increase in raw w spectra signal
fc_highpass=1/700;  % for glider data: 1/5400; % 1800s ~ 0.5h (profiling

% % Fd=designfilt('highpassiir', 'PassbandFrequency', 0.011, 'StopbandFrequency', 0.009, 'PassbandRipple', 1, 'StopbandAttenuation', 60, 'SampleRate', 0.125);
% Fd_bp=designfilt('bandpassiir', 'FilterOrder', 4, 'PassbandFrequency1', 1/5400, 'PassbandFrequency2', 1/100, 'PassbandRipple', 0.1, 'SampleRate', fs);

Fd_lp=designfilt('lowpassiir', 'FilterOrder', 4, 'PassbandFrequency', fc_lowpass,'PassbandRipple', 0.1, 'SampleRate', fs); %low-pass filter to remove high-freq noise <100s
Fd_hp=designfilt('highpassiir', 'FilterOrder', 4, 'PassbandFrequency', fc_highpass,'PassbandRipple', 0.1, 'SampleRate', fs); %low-pass filter to remove low freq from inacuracies
% of the model 5400s cutoff

W=cell(3,1);
D=cell(3,1);
P=cell(3,1);

D{1}=data1; D{2}=data2; D{3}=data3;
P{1}=pres1; P{2}=pres2; P{3}=pres3;

% Fs{1}=fs1; % if plotting spectra from non-interpolated data. specify freq
% for each
% Fs{2}=fs2;
% Fs{3}=fs3;

W_detrend=cell(3,1);
W_detrend_lp=cell(3,1);
W_detrend_lp_hp=cell(3,1);

for d=1:3
    
    data=D{d};
    pres=P{d};
    
    [nZ,nP]=size(data);
    pfit=nan(nZ,nP);
    w_meas_dtrend=nan(nZ,nP);
%     w_meas_lp=nan(nZ,nP);
%     w_meas_filt=nan(nZ,nP);
%     w_water_hpfilt=nan(nZ,nP);
    
    for p=1:nP
        
        indnani=find(~isnan(data(:,p)),1,'first');
        indnanf=find(~isnan(data(:,p)),1,'last');
        
        % Quadratic detrending
        x=pres(indnani:indnanf,1);
        y=data(indnani:indnanf,p);
        
        if length(y)>15
            
            % % Quadratic polynomial fit
            [p_coef, S, mu] = polyfit(x,y,4);
            pfit(indnani:indnanf,p)=polyval(p_coef,x,[],mu);
            
            w_meas_dtrend(indnani:indnanf,p)=y-pfit(indnani:indnanf,p); % detrend low-frequencies introduced by buoyancy control?
            
            W_detrend{d}=w_meas_dtrend;
            
            % Low Pass (to remove high freq noise >1/100s)
            w_meas_lp(indnani:indnanf,p)=filtfilt(Fd_lp,w_meas_dtrend(indnani:indnanf,p)); %low-pass to remove high-freq noise (T>100s)
            w_meas_filt(indnani:indnanf,p)=filtfilt(Fd_hp,w_meas_lp(indnani:indnanf,p)); % high-pass with cutoff 1/5400s (~1.5h), after low-pass was applied

            W_detrend_lp{d}=w_meas_lp;
            W_detrend_lp_hp{d}=w_meas_filt;

            %     % Band-Pass filter (1/100s < f < 1/5400);
            %      w_water_hpfilt(indnani:indnanf,p)=filtfilt(Fd_bp,w_meas_dtrend(indnani:indnanf,p));
            %      % agrees with above , except on edges but above method agrees better
            %      with the model based w
            
            
            %     % 4) Finally, data above 125 m and within 75m of the maximum dive depth are excluded to remove accelerations/decelarations not due to water motion
            %     % and to avoid filter edge effects
            %     w_meas_filt(indnani:indnani+125,p)=NaN;
            %     w_meas_filt(indnanf-25:indnanf,p)=NaN;
            
        else
            continue
        end
    end
end

%%  SELECT data to use in SPECTRA plot
% Estimated error from R2013, w_water=0.005 m/s
% should compute spectra to estimate the noise level
% w_error=zeros(nZ-1,nP);
% w_error(w_meas_dtrend<=0.002)=1;

%% Detrended + high-passed (to remove surface wave effectts)
data1=W_detrend{1}; 
data2=W_detrend{2}; 
data3=W_detrend{3}; method='detrend';

%% Detrended + high-passed (to remove surface wave effectts)
data1=W_detrend_lp{1}; 
data2=W_detrend_lp{2}; 
data3=W_detrend_lp{3}; method='detrend_lp';

%% Compute spectral estimates using pwelch (and default values
% NOTE: I'm subtracting a mean.. my spectra code does it inside the
% function.

figure(1)
x1=data1(:); 
nfft1=length(data1);
ns1=length(x1)/length(data1);
window1=hann(nfft1); noverlap1=nfft1/2;
% fs1=1; % 1 Hz
% [pxx1,f1] = pwelch(x1,window,noverlap,nfft,fs);
% pwelch(x1,window1,noverlap1,nfft1,fs1);
[pxx1,f1] = pwelch(x1,window1,noverlap1,nfft1,fs1);
% hold on

x2=data2(:);
nfft2=length(data2);
ns2=length(x2)/length(data2);
window2=hann(nfft2); noverlap2=nfft2/2;
% fs2=1/10; % for mean dt ~10s (when binning data int0 2-m bins, and profiling at 0.15 m/s)
% pwelch(x2,window2,noverlap2,nfft2,fs2);
[pxx2,f2] = pwelch(x2,window2,noverlap2,nfft2,fs2);

x3=data3(:);
nfft3=length(data3);
ns3=length(x3)/length(data3);
window3=hann(nfft3); noverlap3=nfft3/2;
% fs3=1./(0.3/0.15); % from MR, sampled at 64Hz, but 1 sample every 0.3dbar profiling at 0.15 dbar/s
% pwelch(x3,window3,noverlap3,nfft3,fs3);
[pxx3,f3] = pwelch(x3,window3,noverlap3,nfft3,fs3);

% legend('CTD 1hz','CTD 2m-binned','MR 64 Hz')

% Plot with freq axis in NON-log scale
close all
semilogy(f1,pxx1); hold on
semilogy(f2,pxx2)
semilogy(f3,pxx3)
legend('CTD 1hz','CTD 2m-binned','MR P 64 Hz');

xlabel('Frequency (cps)')
ylabel('Power Spectra (m/s)^2/cps')

ylim([10^(-4) 1])

xlim([10^(-3) 10^(-1)])

print_setup(40,30);
% print([path_plots,'Spectra_3ws_semilogy_pwelch_',input,'_',method,'.jpg'],'-djpeg');

%% Plot with freq axis in a log-log scale
close all

% % if plotting with original matlab colors
% loglog(f1,pxx1); hold on
% loglog(f2,pxx2)
% loglog(f3,pxx3)
% legend('CTD 1hz','CTD 2m-binned','MR P 64 Hz')

% overlay colors in dpdt_ctdmr_comparison.m spectra

% h3=loglog(kn_W1hz,meanSE_W1hz,'color',lgreen);
% h2=loglog(kn_2mbin,meanSE_2mbin,'color',lblue);
% h1=loglog(kn_MRw,meanSE_MRw,'color',lred);

feff=1/10;

h1=loglog(f1,pxx1,'color',lgreen); hold on
h2=loglog(f2,pxx2,'color',lblue);
h3=loglog(f3,pxx3,'color',lred);

h1=loglog(f1(f1<feff),pxx1(f1<feff),'color',green2);
h2=loglog(f2(f2<feff),pxx2(f2<feff),'color',blue);
h3=loglog(f3(f3<feff),pxx3(f3<feff),'color',red);

legend('','','','CTD 1Hz','2-m binned CTD 1Hz','P from MR 64Hz');

grid on

xlabel('Frequency (cps)')
ylabel('Power Spectra (m/s)^2/cps')

ylim([10^(-4) 1]) % limits for w spectra in Merckelbach et al
ylim([10^(-5) 1]) % limits for w spectra in Rudnick
ylim([10^(-5) 10^(-2)])
% xlim([2*10^(-3) 10^(-1)]) % % upper limit for w spectra in Rudnick et al 2013.


ylim([10^(-8) 10^(-1)])
xlim([10^(-3) 10^(0)])

print_setup(40,30);
% print([path_plots,'Spectra_3ws_loglog_pwelch_',input,'_',method,'_cps_ncb.jpg'],'-djpeg');

%% plot Spectra in cpm
close all
PV=0.18;
loglog(f1./PV,pxx1); hold on
loglog(f2./PV,pxx2)
loglog(f3./PV,pxx3)
legend('CTD 1hz','CTD 2m-binned','MR P 64 Hz')

ylim([10^(-5) 1]) % 
ylim([10^(-5) 10^(-2)])

xlim([10^(-2) 10^(0)])

xlabel('Frequency (cpm)')
ylabel('Power Spectra (m/s)^2/cpm')

print_setup(40,30);
% print([path_plots,'Spectra_3ws_loglog_pwelch_',input,'_',method,'_cpm.jpg'],'-djpeg');

%% plot Spectra in rad/s
close all
c2rad=2*pi;
loglog(f1*c2rad,pxx1); hold on
loglog(f2*c2rad,pxx2)
loglog(f3*c2rad,pxx3)
legend('CTD 1hz','CTD 2m-binned','MR P 64 Hz')

xlim([10^(-2) 10^(0)])

xlabel('Frequency (rad/s)')
ylabel('Power Spectra (m/s)^2/rad/s')

print_setup(40,30);
% print([path_plots,'Spectra_3ws_loglog_pwelch_',input,'_',method,'_rad.jpg'],'-djpeg');

%% Check PROFILE data
%% Plot FILTERING technique for a given profile?

p=15;

close all
subplot(131)
plot(data1(:,p),pres1); 
axis ij
ylabel('Depth');
xlabel('w_{1hz} [m/s]');

xlim([-0.015 0.015])
ylim([0 100])
print_setup(20,40)
print([path_plots,'Profiles_w1hz_',input,'_',method,'.jpg'],'-djpeg');

subplot(132)
plot(data2(:,p),pres2); 
axis ij
ylabel('Depth');
xlabel('w_{2mb} [m/s]');

ylim([0 100])
xlim([-0.015 0.015])
print_setup(20,40)
print([path_plots,'Profiles_w2mb_',input,'_',method,'.jpg'],'-djpeg');

subplot(133)
plot(data3(:,p),pres3); 
axis ij
ylabel('Depth');
xlabel('w_{mr} [m/s]');

xlim([-0.015 0.015])
ylim([0 100])
print_setup(60,40)
% print([path_plots,'Profiles_filtering_technique_',input,'_',method,'_p',num2str(p),'.jpg'],'-djpeg');

%% 
% %%  compare my spectra.m code vs pwelch
% W_MR=W_mr(10:200,:);
% x=detrend(W_MR,1);
%
%
% Nx = length(x);
% nsc = floor(Nx/18);
%
% nu=18; % degrees of freedom
% ns=nu/2; % number of segments
% nsc = floor(Nx/ns);
%
% nov = floor(nsc/2); % for 50% overlap
% nfft = 382; %nsc; %max(256,2^nextpow2(nsc));
%
% fs=64;
% fs=0.3750
% % nfft=191; % (fft length)
%
% probability=0.95;
%
% figure(1)
% % [pxx,f] = pwelch(x,window,noverlap,f,fs)
% [pxx,f] = pwelch(detrend(x));
% [pxx,f] = pwelch(detrend(x),hanning(nsc),'',nfft,fs,'ConfidenceLevel',probability);
%
%
% figure(2)
% loglog(f,pxx)
%
% hold on
%
% % Compute Spectra for 64Hz dp/dt (smoothed by RSI)
% dst=1/64; % % in sec for 64Hz sampling rate.
% % % dst=0.0982; % % Effective resolution? for a generic low-pass filter: pi* DeltaR_raw/w_cuttoff
% % NOTE: RSI applies a butterworth LP filter! and the data seems to be
% % 'binned' in 0.3-0.4m bins (that's the effective resolution, 1 sample per
% dst=(0.3/0.15);
%
% % Spectra code uses a mean dt
% % time=dst; %10; %/3600;
% % data=MR.speed(:,~isnan(MR.speed(1,:))); % remove profiles w/NaN
% % data=data(:,:); % truncate to omit bottom layer (BC ON)
%
% % if using interpolated data
% time=dst;
% data=W_mr(10:200,:); [nz,np]=size(data);
%
% nu=np*2; % to compute spectra for each profile, and average across profiles
% varunits='m/s';
% timedt='s';
% color='r';
% [kn_MRw,meanSE_MRw,dk,varcheckSat,h,SE_MRw]=spectra(time,data,nu,varunits,timedt,color);
%
% % % Change color of lines overlaying plot
% % h3=loglog(kn_W1hz,meanSE_W1hz,'color',green1);
% % h2=loglog(kn_2mbin,meanSE_2mbin,'color',blue);
%
% meanSE_pwelch=pxx;
% kn_pwelch=f;
%
% h1=loglog(kn_pwelch,meanSE_pwelch,'color','b');
% h2=loglog(kn_MRw,meanSE_MRw,'color','r');
%
% legend('','','pwelch','fft code');
%
% % maximize;
% print_setup(40,30)
% print([path_plots,'Spectra_wctd_MR_pwelch_vs_fft.jpg'],'-djpeg');

%% PDF comparison of w from the three signals
close all

% define colors;
cmap=colormap(m_colmap('jet','step',10));
grey=[0.3 0.3 0.3];

blue=cmap(1+20,:);
green=cmap(75,:)
red=cmap(end-30,:);

w3=detrend(data3,1);
pdffit6(w3(:),0.0025,green,0.4)
m_w3=nanmean(w3);
rms_w3=nanmean(nanstd(w3,1)); % rms w3 = 0.0071 m/s

w1=detrend(data1,1);
pdffit6(w1(:),0.0025,blue,0.6)
rms_w1=nanmean(nanstd(w1,1)); % rms w1 = 0.0118 m/s
hold on

w2=detrend(data2,1);
pdffit6(w2(:),0.0025,red,0.8)
rms_w2=nanmean(nanstd(w2,1)); % rms w2 = 0.0084 m/s

xlim([-0.06 0.06])
xlabel([' w ',method,' (m/s)']);
legend('MR P 64Hz','CTD 1Hz','CTD 1Hz, 2-m binned','Location','NW');

print_setup(40,30)
print([path_plots,'PDF_wctd_MR_',method,'.jpg'],'-djpeg');


% % %%% From Glider processing code:
% % % Fd=designfilt('highpassiir', 'PassbandFrequency', 0.011, 'StopbandFrequency', 0.009, 'PassbandRipple', 1, 'StopbandAttenuation', 60, 'SampleRate', 0.125);
% % Fd_bp=designfilt('bandpassiir', 'FilterOrder', 4, 'PassbandFrequency1', 1/5400, 'PassbandFrequency2', 1/100, 'PassbandRipple', 0.1, 'SampleRate', 0.125);
% % 
% % Fd_lp=designfilt('lowpassiir', 'FilterOrder', 4, 'PassbandFrequency', 0.01,'PassbandRipple', 0.1, 'SampleRate', 0.125); %low-pass filter to remove high-freq noise <100s
% % Fd_hp=designfilt('highpassiir', 'FilterOrder', 4, 'PassbandFrequency', 1/5400,'PassbandRipple', 0.1, 'SampleRate', 0.125); %low-pass filter to remove low freq from inacuracies of the model 5400s cutoff
% % 
% % pfit=nan(nZ-1,nP);
% % w_meas_dtrend=nan(nZ-1,nP);
% % w_meas_lp=nan(nZ-1,nP);
% % w_meas_filt=nan(nZ-1,nP);
% % w_water_hpfilt=nan(nZ-1,nP);
% % 
% % for p=1:nP
% %     
% %     indnani=find(~isnan(w_meas_interp(:,p)),1,'first');
% %     indnanf=find(~isnan(w_meas_interp(:,p)),1,'last');
% %     
% %     % Quadratic detrending
% %     x=du.depth(indnani:indnanf,p);
% %     y=w_meas_interp(indnani:indnanf,p);
% %     
% %     if length(y)>15
% %      
% % %     % Quadratic polynomial fit
% %     [p_coef, S, mu] = polyfit(x,y,4);
% %     pfit(indnani:indnanf,p)=polyval(p_coef,x,[],mu);
% %     
% %     w_meas_dtrend(indnani:indnanf,p)=y-pfit(indnani:indnanf,p); % detrend low-frequencies introdced by buoyancy control
% %     
% %     % High Pass (to remove high freq noise >1/100s)
% %      w_meas_lp(indnani:indnanf,p)=filtfilt(Fd_lp,w_meas_dtrend(indnani:indnanf,p)); %low-pass to remove high-freq noise (T>100s)
% %      w_meas_filt(indnani:indnanf,p)=filtfilt(Fd_hp,w_meas_lp(indnani:indnanf,p)); % high-pass with cutoff 1/5400s (~1.5h), after low-pass was applied
% %      
% % %     % Band-Pass filter (1/100s < f < 1/5400);
% % %      w_water_hpfilt(indnani:indnanf,p)=filtfilt(Fd_bp,w_meas_dtrend(indnani:indnanf,p));
% % %      % agrees with above , except on edges but above method agrees better
% % %      with the model based w
% % 
% % 
% % %     % 4) Finally, data above 125 m and within 75m of the maximum dive depth are excluded to remove accelerations/decelarations not due to water motion 
% % %     % and to avoid filter edge effects
% % %     w_meas_filt(indnani:indnani+125,p)=NaN;
% % %     w_meas_filt(indnanf-25:indnanf,p)=NaN;
% %     
% %     else
% %         continue
% %     end
% % end
% % 
% % % % Estimated error from R2013, w_water=0.005 m/s
% % % % should compute spectra to estimate the noise level
% % % w_error=zeros(nZ-1,nP);
% % % w_error(w_meas_dtrend<=0.002)=1;
