%% Behav states AFR1

dir = 'D:\Data\AFR1\AFR1_140118\AFRtriple1_180114_'; 
file1 = readmulti([dir 'cage.eeg'],64);
%file2 = readmulti([dir 'field.eeg'],64);
%eegcat = cat(1,file1,file2);
Fs = 1250;
%t1 = size(file1,1); % length sleep session
%LFPtheta = eegcat(:,15);

%tfrag1 = 2000*1250; %2000
tfrag2 = 1150*1250; %end
LFPtheta = file1(1:tfrag2,15);

%% Find theta epochs
TDratio_thresh = 1.5; % lo mas optimo me parece 1.0 para MEC

t = (0:(size(LFPtheta,1)-1))'/Fs;
Fs_spec = Fs/10;
t_spec = downsample(t,10);
[thetaSpec,theta_f] = long_wavespec(downsample(LFPtheta, 10), t_spec, 10, [1 20], 1);
nthetaAvg = 2*round(0.25*Fs_spec/2) - 1;
smooth_sec = 2;
theta_Finds = find_inds(5, theta_f):find_inds(8, theta_f);
delta_Finds = find_inds(1, theta_f):find_inds(4, theta_f);
thetaDeltaRatio = runavg([zeros(round(smooth_sec*Fs_spec/2),1); mean(thetaSpec(theta_Finds,:),1)'./mean(thetaSpec(delta_Finds,:),1)'; ...
    zeros(round(smooth_sec*Fs_spec/2),1)], round(smooth_sec*Fs_spec/2)*2+1);
thetaDeltaInterp = [lin_interp(thetaDeltaRatio, 10); ones(9,1)*thetaDeltaRatio(end)];
thetaBools = thetaDeltaInterp > TDratio_thresh;
thetaInds = find(thetaBools);
if length(thetaBools) > length(LFPtheta)
thetaBools=thetaBools(1:length(LFPtheta));
end
NOthetaInds = find(not(thetaBools));

% plot para comprobar resultado
%     params = struct ('tapers',{[3 5]},'pad',{0},'Fs',{1250},'fpass',{[1 30]},'err',{[1 0.95]},'trialave',{1});
%     [spect1,fbins1,Serr1] = mtspectrumc(LFPtheta(thetaInds,:),params);
%     [spect2,fbins2,Serr2] = mtspectrumc(LFPtheta(NOthetaInds,:),params);
%     figure;plot_vector(spect1,fbins1,'n',Serr1,'b',2);hold on;
%            plot_vector(spect2,fbins2,'n',Serr2,'r',2);hold on;
%     set(gca,'XLim',[params.fpass(1) params.fpass(2)]);

%% Create STATES variable 

% tsleep = t1; % sep sleep session / run session
% tsT = find(thetaInds >= tsleep);
% tsT = tsT(1);
% tsNT = find(NOthetaInds >= tsleep);
% tsNT = tsNT(1);
% 
% STATES = struct;
% STATES.REM = 'REM';
% STATES.REM.inds_all = thetaInds(1:tsT);
% STATES.SWS = 'SWS';
% STATES.SWS.inds_all = NOthetaInds(1:tsNT);
% STATES.RUN = 'RUN';
% STATES.RUN.inds_all = thetaInds(tsT:end);
% STATES.QUIET = 'QUIET';
% STATES.QUIET.inds_all = NOthetaInds(tsNT:end);

%save([dir filename '_STATES.mat'], 'STATES');

%% Detec theta peaks to make averaged CSD
ca1_theta = filtsig(LFPtheta, 1000/Fs, [2 4 12 14]/1000);
ca1_hilb = hilbert(ca1_theta);
thetaP = angle(ca1_hilb(thetaInds));
thetaH = abs(ca1_hilb(thetaInds));

p=[];
for i=1:length(thetaP)
    if thetaP(i,1) > 3.1 
        p(i) = 1 ;
    end
end
tet = find(p);
ventana = floor(0.05*Fs); % 50ms
for j=1:length(p)-ventana
    if p(j)==1
%        for k=j+1:j+ventana
       p(j+find(p((j+1):(j+ventana)))) = 0;
    end
end
tetok=find(p);
tetokT= tetok/Fs; 

% STATES.REM.theta_peaks = intersect(STATES.REM.inds_all,tetok');
% STATES.RUN.theta_peaks = intersect(STATES.RUN.inds_all,tetok');
% save([dir filename '_STATES.mat'], 'STATES');

% %salvar como tren spikes
Trains=struct();
Trains.Spikes = {tetokT'};
%Trains.Spikes = {STATES.REM.theta_peaks/1250,STATES.RUN.theta_peaks/1250,};
Trains.Titles = char('01');
