%SU modulation
% v2: no pwelch, add fft normalized and dif
% xcorr: with cross correlation version
% combine: combines multiple units & time frames of the same state
% combine v2: refined save .mat
% SW: draw each individual figures and save to folders(label1 &label2 instead of no tremor/tremor)
%   : 2022.04.18 v2: small bug fixed(displaying tremor time frames for no_tremor time frames on figs)
% SU-modulation_vs_FPfreq: 
%      2022.05.11, modified from xcorr combine, but added firing freq & instFR intp and all compare against FP tremor freq

clear all; close all;
%% Input Parameters
kilo=1; %from kilosort or not
NEVchannels=[1,1,1,1,1,1,1]; % Can be multiple channels, matrix length match NEVunits
NEVunits=[6,7,8,9,10,11,12]; % Corresponding detected units of each channel, can be more than one unit in each channel

FPchannel=129; 
SumWindow=[5,55]; % Define the frequency range for normalize the data 
ExtractFreqRange=[5,40]; % Define the frequency range to extract the data

%time_frame_cell={[0,120;3031,3112];[0,120;1639,1776];[0,120;1945,2059];[0,120;2340,2400]}
fft_ylim=0.05;

dsfs=250; %down-sampled SU (hz)
maxlag=1; % xcorr maxlag (sec)

L1=[177,252;503,538]; % actually means Label 1
L2=[68,119;409,455;565,621]; % actually means Label 2

max_range=[5,50]; % in hz, range to find local maximum

label1='Rest';
label2='Move';

lim=40; % plot freq lim, in hz
%%
GetFilePath=[pwd '\SU modulations & FR vs Tremor freq\']
% SubPath=[sprintf('Ch%d U%d modulation %d-%d_%d-%d',NEVchannels,NEVunits,time_frame(1,1),time_frame(1,2),time_frame(2,1),time_frame(2,2))]; 
% SubPath=[SubPath,'\']
% FullFilePath=[GetFilePath,SubPath];
mkdir(GetFilePath)

%% NEV units spike time frame
if kilo==0
    NEV = F_openNEX;
elseif kilo==1
    NEV= F_openKilosort2;
end

fs=NEV.MetaTags.SampleRes;

for h=1:length(NEVunits);  
    SU_Name(h)={sprintf('%d Chan%d-U%d',h,NEVchannels(h),NEVunits(h))};
    C=find(NEV.Data.Spikes.Electrode==NEVchannels(h)); % find spikes index that match the selected ChannelID
    D=find(NEV.Data.Spikes.Unit(C)==NEVunits(h)); % find spikes index that match the sorted UnitID. 
    SpikeFrames{h,1}=NEV.Data.Spikes.TimeStamp(C(D)); % find spikes that match the sorted channel-unit.
    dsSpikeFrames{h,1}=floor(SpikeFrames{h,1}/(fs/dsfs));
    
    nt_corr_dsraw{1,h+1}=SU_Name(h);
    nt_lags_dsraw{1,h+1}=SU_Name(h);
    nt_ds_pseudo_data{1,h+1}=SU_Name(h);
    nt_fft{1,h+1}=SU_Name(h);
    nt_fr{1,h+1}=SU_Name(h);
    
    t_corr_dsraw{1,h+1}=SU_Name(h);
    t_lags_dsraw{1,h+1}=SU_Name(h);
    t_ds_pseudo_data{1,h+1}=SU_Name(h);
    t_fft{1,h+1}=SU_Name(h);
    t_fr{1,h+1}=SU_Name(h);
end
nt_fft{1,length(NEVunits)+2}='Combined';
t_fft{1,length(NEVunits)+2}='Combined';
nt_fft{length(L1)+2,1}='Combined';
t_fft{length(L2)+2,1}='Combined';

Fs=dsfs;
L=maxlag*dsfs*2+1;
f = Fs*(0:(L/2))/L;
max_range_index=[find(f>max_range(1),1),find(f<max_range(2),1,'last')];

tempsum_intime=zeros(1,length(f));
for x=1:1:length(L1)
    time_frame=L1(x,:);
    nt_corr_dsraw{x+1,1}=time_frame;
    nt_lags_dsraw{x+1,1}=time_frame;
    nt_ds_pseudo_data{x+1,1}=time_frame;
    nt_fft{x+1,1}=time_frame;
    
    tempsum=zeros(1,length(f));
    for i=1:length(NEVunits)
    nt_ds_pseudo_data{x+1,i+1}=zeros(1,(time_frame(2)-time_frame(1))*dsfs);
    ds_range=[find(dsSpikeFrames{i,1}>time_frame(1)*dsfs,1), find(dsSpikeFrames{i,1}<time_frame(2)*dsfs,1,"last")];
    nt_ds_pseudo_data{x+1,i+1}(dsSpikeFrames{i,1}(ds_range(1):ds_range(2))-time_frame(1)*dsfs)=1;
    [nt_corr_dsraw{x+1,i+1} nt_lags_dsraw{x+1,i+1}]=xcorr(nt_ds_pseudo_data{x+1,i+1},maxlag*dsfs,'coeff');
    
    Y=fft(nt_corr_dsraw{x+1,i+1});
    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    nt_fft{x+1,i+1}=P1;
    tempsum=tempsum+nt_fft{x+1,i+1};
    
    end
    nt_fft{x+1,length(NEVunits)+2}=tempsum/length(NEVunits);
    tempsum_intime=tempsum_intime+nt_fft{x+1,length(NEVunits)+2};
end
nt_fft{length(L1)+2,length(NEVunits)+2}=tempsum_intime/length(L1);

tempsum_intime=zeros(1,length(f));
for x=1:1:length(L2)
    time_frame=L2(x,:);
    t_corr_dsraw{x+1,1}=time_frame;
    t_lags_dsraw{x+1,1}=time_frame;
    t_ds_pseudo_data{x+1,1}=time_frame;
    t_fft{x+1,1}=time_frame;
    
    tempsum=zeros(1,length(f));
    for i=1:length(NEVunits)
    t_ds_pseudo_data{x+1,i+1}=zeros(1,(time_frame(2)-time_frame(1))*dsfs);
    ds_range=[find(dsSpikeFrames{i,1}>time_frame(1)*dsfs,1), find(dsSpikeFrames{i,1}<time_frame(2)*dsfs,1,"last")];
    t_ds_pseudo_data{x+1,i+1}(dsSpikeFrames{i,1}(ds_range(1):ds_range(2))-time_frame(1)*dsfs)=1;
    [t_corr_dsraw{x+1,i+1} t_lags_dsraw{x+1,i+1}]=xcorr(t_ds_pseudo_data{x+1,i+1},maxlag*dsfs,'coeff');
    
    Y=fft(t_corr_dsraw{x+1,i+1});
    P2 = abs(Y/L);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    t_fft{x+1,i+1}=P1;
    tempsum=tempsum+t_fft{x+1,i+1};
    
    end
    t_fft{x+1,length(NEVunits)+2}=tempsum/length(NEVunits);
    tempsum_intime=tempsum_intime+t_fft{x+1,length(NEVunits)+2};
end
t_fft{length(L2)+2,length(NEVunits)+2}=tempsum_intime/length(L2);

for i=1:length(L1)
    for j=1:length(NEVunits)+1
        [max_nt(i,j) index_nt(i,j)]=max(nt_fft{i+1,j+1}(max_range_index(1):max_range_index(2)));
        index_nt(i,j)=f(index_nt(i,j)+max_range_index(1)-1);
    end
end

for i=1:length(L2)
    for j=1:length(NEVunits)+1
        [max_t(i,j) index_t(i,j)]=max(t_fft{i+1,j+1}(max_range_index(1):max_range_index(2)));
        index_t(i,j)=f(index_t(i,j)+max_range_index(1)-1);
    end
end

%% firing freq
for x=1:length(L1)
    time_frame=L1(x,:);
    nt_fr{x+1,1}=time_frame;
    for i=1:length(NEVunits)
        ds_range=[find(dsSpikeFrames{i,1}>time_frame(1)*dsfs,1), find(dsSpikeFrames{i,1}<time_frame(2)*dsfs,1,"last")];
        nt_fr{x+1,i+1}=(ds_range(2)-ds_range(1)+1)/(time_frame(2)-time_frame(1));   
    end
end
for x=1:length(L2)
    time_frame=L2(x,:);
    t_fr{x+1,1}=time_frame;
    for i=1:length(NEVunits)
        ds_range=[find(dsSpikeFrames{i,1}>time_frame(1)*dsfs,1), find(dsSpikeFrames{i,1}<time_frame(2)*dsfs,1,"last")];
        t_fr{x+1,i+1}=(ds_range(2)-ds_range(1)+1)/(time_frame(2)-time_frame(1));   
    end
end

%% instFR intp SU modulation

%% tremor freq
[NsFileName NsFilePath]=uigetfile('*.ns*'); % get filename and path
NSx=openNSx([NsFilePath NsFileName]);
data=double(NSx.Data');
NSx=rmfield(NSx,'Data');
channel=find(NSx.MetaTags.ChannelID==FPchannel);
for x=1:length(L1)
    time_frame=L1(x,:);
    NT_FPfreqMax{x,1}=strcat(string(time_frame(1)),"-",string(time_frame(2)));
    NT_FPvalMax{x,1}=time_frame;
    [NT_FPfreqMax{x,2},NT_FPvalMax{x,2}]=find_freqmax(1,data,channel,time_frame,1,SumWindow,ExtractFreqRange);
end

for x=1:length(L2)
    time_frame=L2(x,:);
    T_FPfreqMax{x,1}=strcat(string(time_frame(1)),"-",string(time_frame(2)));
    T_FPvalMax{x,1}=time_frame;
    [T_FPfreqMax{x,2},T_FPvalMax{x,2}]=find_freqmax(1,data,channel,time_frame,1,SumWindow,ExtractFreqRange);
end
%% FSI (frequency similarity index)
FSI{1,2}="FR";
FSI{1,3}="index";
FSI{1,4}="Xcorr Mod";
FSI{1,5}="index";
FSI{1,6}="Xcorr Mod(combined)";
FSI{1,7}="index";
FSI{2,1}=label1;
FSI{3,1}=label2;
%% plot(FR)
cd(GetFilePath)

figure()
hold on
tempx=reshape(repmat(cell2mat(NT_FPfreqMax(:,2)),1,length(NEVunits))',[1,length(NEVunits)*length(L1)]);
tempy=reshape(cell2mat(nt_fr(2:end,2:end))',[1,length(NEVunits)*length(L1)]);
scatter(tempx,tempy)

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])
hold off

FSI{2,2}=(tempy-tempx).^2;
FSI{2,3}=mean(FSI{2,2});

xlabel('Tremor Frequency (hz)')
ylabel('Firing Frequency (hz)')
suptitle([label1,'-',NsFileName])
title(sprintf('FSI:%.2f Firing Frequency vs. Tremor Frequency-%d units %d time frames',FSI{2,3},length(NEVunits),length(L1)))
saveas(gcf,sprintf('%s_FRvsTremorFreq.fig',label1))
saveas(gcf,sprintf('%s_FRvsTremorFreq.jpg',label1))

figure()
tempx=reshape(repmat(cell2mat(T_FPfreqMax(:,2)),1,length(NEVunits))',[1,length(NEVunits)*length(L2)]);
tempy=reshape(cell2mat(t_fr(2:end,2:end))',[1,length(NEVunits)*length(L2)]);
scatter(tempx,tempy)

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])

FSI{3,2}=(tempy-tempx).^2;
FSI{3,3}=mean(FSI{3,2});

xlabel('Tremor Frequency (hz)')
ylabel('Firing Frequency (hz)')
suptitle([label2,'-',NsFileName])
title(sprintf('FSI:%.2f Firing Frequency vs. Tremor Frequency-%d units %d time frames',FSI{3,3},length(NEVunits),length(L2)))
saveas(gcf,sprintf('%s_FRvsTremorFreq.fig',label2))
saveas(gcf,sprintf('%s_FRvsTremorFreq.jpg',label2))

%% plot(xcorr)
figure()
tempx=reshape(repmat(cell2mat(NT_FPfreqMax(:,2)),1,length(NEVunits))',[1,length(NEVunits)*length(L1)]);
tempy=reshape(index_nt(:,1:end-1)',[1,length(NEVunits)*length(L1)]);
scatter(tempx,tempy);

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])

FSI{2,4}=(tempy-tempx).^2;
FSI{2,5}=mean(FSI{2,4});

xlabel('Tremor Frequency (hz)')
ylabel('Modulation Peak Frequency (hz)')
suptitle([label1,'-',NsFileName])
title(sprintf('FSI:%.2f Modulation Peak Frequency vs. Tremor Frequency-%d units %d time frames',FSI{2,5},length(NEVunits),length(L1)))
saveas(gcf,sprintf('%s_XcorrModvsTremorFreq.fig',label1))
saveas(gcf,sprintf('%s_XcorrModvsTremorFreq.jpg',label1))

figure()
tempx=reshape(repmat(cell2mat(T_FPfreqMax(:,2)),1,length(NEVunits))',[1,length(NEVunits)*length(L2)]);
tempy=reshape(index_t(:,1:end-1)',[1,length(NEVunits)*length(L2)]);
scatter(tempx,tempy);

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])

FSI{3,4}=(tempy-tempx).^2;
FSI{3,5}=mean(FSI{3,4});

xlabel('Tremor Frequency (hz)')
ylabel('Modulation Peak Frequency (hz)')
suptitle([label2,'-',NsFileName])
title(sprintf('FSI:%.2f Modulation Peak Frequency vs. Tremor Frequency-%d units %d time frames',FSI{3,5},length(NEVunits),length(L2)))
saveas(gcf,sprintf('%s_XcorrModvsTremorFreq.fig',label2))
saveas(gcf,sprintf('%s_XcorrModvsTremorFreq.jpg',label2))
%% plot(xcorr, combined)
figure()
tempx=reshape(cell2mat(NT_FPfreqMax(:,2)),[1,length(L1)]);
tempy=reshape(index_nt(:,end)',[1,length(L1)]);
scatter(tempx,tempy);

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])

FSI{2,6}=(tempy-tempx).^2;
FSI{2,7}=mean(FSI{2,6});

xlabel('Tremor Frequency (hz)')
ylabel('Modulation Peak Frequency(Combined) (hz)')
suptitle([label1,'-',NsFileName])
title(sprintf('FSI:%.2f Modulation Frequency(Combined) vs. Tremor Frequency-%d units %d time frames',FSI{2,7},length(NEVunits),length(L1)))
saveas(gcf,sprintf('%s_XcorrMod(combined)vsTremorFreq.fig',label1))
saveas(gcf,sprintf('%s_XcorrMod(combined)vsTremorFreq.jpg',label1))

figure()
tempx=reshape(cell2mat(T_FPfreqMax(:,2)),[1,length(L2)]);
tempy=reshape(index_t(:,end)',[1,length(L2)]);
scatter(tempx,tempy);

X = [ones(length(tempx),1) tempx'];
b = X\tempy';
line([0 lim],[b(1) b(1)+lim*b(2)],'Color','blue')
text(lim*0.55,lim*0.4,sprintf('y=%.4fx+%.4f',b(2),b(1)))
line([0 lim],[0 lim],'Color','red','LineStyle','--')
xlim([0 lim])
ylim([0 100])

FSI{3,6}=(tempy-tempx).^2;
FSI{3,7}=mean(FSI{3,6});

xlabel('Tremor Frequency (hz)')
ylabel('Modulation Peak Frequency(Combined) (hz)')
suptitle([label2,'-',NsFileName])
title(sprintf('FSI:%.2f Modulation Frequency(Combined) vs. Tremor Frequency-%d units %d time frames',FSI{3,7},length(NEVunits),length(L2)))
saveas(gcf,sprintf('%s_XcorrMod(combined)vsTremorFreq.fig',label2))
saveas(gcf,sprintf('%s_XcorrMod(combined)vsTremorFreq.jpg',label2))
%% Results(important values used to draw figures)
Report.FR.label1=nt_fr;
Report.FR.label2=t_fr;
Report.Xcorr.label1_fftraw=nt_fft;
Report.Xcorr.label1_maxfreq=index_nt;
Report.Xcorr.label2_fftraw=t_fft;
Report.Xcorr.label2_maxfreq=index_t;
Report.TremorFreq.label1_freqMax=NT_FPfreqMax;
Report.TremorFreq.label1_valMax=NT_FPvalMax;
Report.TremorFreq.label2_freqMax=T_FPfreqMax;
Report.TremorFreq.label2_valMax=T_FPvalMax;
Report.FSI=FSI;
save('Freq_vs_TremorFreq.mat','Report')
%% Max value and freq for all combined(units & timeframes), not useful for now
% [M_nt,I_nt]=max(nt_fft{length(L1)+2,length(NEVunits)+2}(max_range_index(1):max_range_index(2)));
% [M_t,I_t]=max(t_fft{length(L2)+2,length(NEVunits)+2}(max_range_index(1):max_range_index(2)));
% I_nt=f(I_nt+max_range_index(1)-1);
% I_t=f(I_t+max_range_index(1)-1);
%% Xcorr combine ploting, disabled for now
% %% plot all combine
% figure()
% subplot(2,1,1)
% plot(f,nt_fft{length(L1)+2,length(NEVunits)+2})
% title([NEV.MetaTags.Filename,sprintf(' %s combined Max=%.2f hz, %.3f',label1,I_nt,M_nt)]);
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.6,string(index_nt(:,1:end-1)),'FontSize',8) 
% hold on
% s=size(index_nt);
% for i=1:s(1)
%     for j=1:(s(2)-1)
%       line([index_nt(i,j) index_nt(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
%     end
% end
% line([I_nt I_nt],[0 fft_ylim],'Color','red')
% hold off
% 
% subplot(2,1,2)
% plot(f,t_fft{length(L2)+2,length(NEVunits)+2})
% title([sprintf('%s Combined Max=%.2f hz, %.3f, %d units %d/%d time frames',label2,I_t,M_t,length(NEVunits),length(L1),length(L2))])
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.5,string(index_t(:,1:end-1)),'FontSize',7)
% hold on
% s=size(index_t);
% for i=1:s(1)
%     for j=1:(s(2)-1)
%       line([index_t(i,j) index_t(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
%     end
% end
% line([I_t I_t],[0 fft_ylim],'Color','red')
% hold off
% 
% 
% saveas(gcf,[GetFilePath sprintf('Combine_nt%d_t%d_u%d_xcorr_fft',length(L1),length(L2),length(NEVunits)),'.fig'])
% saveas(gcf,[GetFilePath sprintf('Combine_nt%d_t%d_u%d_xcorr_fft',length(L1),length(L2),length(NEVunits)),'.jpg'])
% 
% %% plot all combine (-NT)
% figure()
% subplot(3,1,1)
% plot(f,nt_fft{length(L1)+2,length(NEVunits)+2})
% title([NEV.MetaTags.Filename,sprintf(' %s combined Max=%.2f hz, %.3f',label1,I_nt,M_nt)]);
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.6,string(index_nt(:,1:end-1)),'FontSize',8) 
% hold on
% s=size(index_nt);
% for i=1:s(1)
%     for j=1:(s(2)-1)
%       line([index_nt(i,j) index_nt(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
%     end
% end
% line([I_nt I_nt],[0 fft_ylim],'Color','red')
% hold off
% 
% subplot(3,1,2)
% plot(f,t_fft{length(L2)+2,length(NEVunits)+2})
% title([sprintf('%s Combined Max=%.2f hz, %.3f, %d units %d/%d time frames',label2,I_t,M_t,length(NEVunits),length(L1),length(L2))])
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.5,string(index_t(:,1:end-1)),'FontSize',7)
% hold on
% s=size(index_t);
% for i=1:s(1)
%     for j=1:(s(2)-1)
%       line([index_t(i,j) index_t(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
%     end
% end
% line([I_t I_t],[0 fft_ylim],'Color','red')
% hold off
% 
% subplot(3,1,3)
% a=t_fft{length(L2)+2,length(NEVunits)+2}-nt_fft{length(L1)+2,length(NEVunits)+2};
% [max_minus index_minus]=max(a(max_range_index(1):max_range_index(2)));
% index_minus=f(index_minus+max_range_index(1)-1);
% plot(f,t_fft{length(L2)+2,length(NEVunits)+2}-nt_fft{length(L1)+2,length(NEVunits)+2})
% title(sprintf('%s-%s, max %.2f hz',label2,label1,index_minus))
% xlim([0 60])
% ylim([0 fft_ylim/10])
% line([index_minus index_minus],[0 fft_ylim/10],'Color','red')
% 
% saveas(gcf,[GetFilePath sprintf('Combine_nt%d_t%d_u%d_xcorr_fft(-nt)',length(L1),length(L2),length(NEVunits)),'.fig'])
% saveas(gcf,[GetFilePath sprintf('Combine_nt%d_t%d_u%d_xcorr_fft(-nt)',length(L1),length(L2),length(NEVunits)),'.jpg'])
% 
% %% combine units only
% for i=1:length(L1)
% figure()
% plot(f,nt_fft{i+1,length(NEVunits)+2})
% title([NEV.MetaTags.Filename,sprintf(' %s %d-%ds All units combined(%dUnits)Max=%.2f hz, %.3f',label1,L1(i,1),L1(i,2),length(NEVunits),max_nt(i,end),index_nt(i,end))]);
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.6,string(index_nt(i,1:end-1)),'FontSize',8) 
% 
% hold on
% s=size(index_nt);
% for j=1:(s(2)-1)
%   line([index_nt(i,j) index_nt(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
% end
% line([index_nt(i,end) index_nt(i,end)],[0 fft_ylim],'Color','red')
% 
% hold off
% 
% saveas(gcf,[GetFilePath sprintf('%s_AllUnitsCombined_%dunits_%d-%d_xcorr_fft',label1,length(NEVunits),L1(i,1),L1(i,2)),'.fig'])
% saveas(gcf,[GetFilePath sprintf('%s_AllUnitsCombined_%dunits_%d-%d_xcorr_fft',label1,length(NEVunits),L1(i,1),L1(i,2)),'.jpg'])
% end
% 
% for i=1:length(L2)
% figure()
% plot(f,t_fft{i+1,length(NEVunits)+2})
% title([NEV.MetaTags.Filename,sprintf(' %s %d-%ds All units combined(%dUnits)Max=%.2f hz, %.3f',label2,L2(i,1),L2(i,2),length(NEVunits),max_t(i,end),index_t(i,end))]);
% xlim([0 60])
% ylim([0 fft_ylim])
% text(50,fft_ylim*0.6,string(index_t(i,1:end-1)),'FontSize',8) 
% 
% hold on
% s=size(index_t);
% for j=1:(s(2)-1)
%   line([index_t(i,j) index_t(i,j)],[0 fft_ylim],'Color','red','LineStyle','--') 
% end
% line([index_t(i,end) index_t(i,end)],[0 fft_ylim],'Color','red')
% 
% hold off
% 
% saveas(gcf,[GetFilePath sprintf('%s_AllUnitsCombined_%dunits_%d-%d_xcorr_fft',label2,length(NEVunits),L2(i,1),L2(i,2)),'.fig'])
% saveas(gcf,[GetFilePath sprintf('%s_AllUnitsCombined_%dunits_%d-%d_xcorr_fft',label2,length(NEVunits),L2(i,1),L2(i,2)),'.jpg'])
% end
% 
% 
% %% individual figures
% subpath=[GetFilePath sprintf('%s',label1)];
% subpath=[subpath '\'];
% mkdir(subpath)
% 
% for i =1:length(NEVunits)
%     for j=1:length(L1)
%         figure()
%         plot(f,nt_fft{j+1,i+1})
%         xlim([0,60])
%         ylim([0,fft_ylim])
%         xlabel('hz')
%         title(sprintf('%s Ch%d U%d %d-%d %s',label1,NEVchannels(i),NEVunits(i),L1(j,1),L1(j,2),NEV.MetaTags.Filename))
% 
%         saveas(gcf,[subpath sprintf('%s_Ch%d_U%d_%d-%ds',label1,NEVchannels(i),NEVunits(i),L1(j,1),L1(j,2)),'.fig'])
%         saveas(gcf,[subpath sprintf('%s_Ch%d_U%d_%d-%ds',label1,NEVchannels(i),NEVunits(i),L1(j,1),L1(j,2)),'.jpg'])
% 
%     end
% end
% 
% subpath=[GetFilePath sprintf('%s',label2)];
% subpath=[subpath '\'];
% mkdir(subpath)
% 
% for i =1:length(NEVunits)
%     for j=1:length(L2)
%         figure()
%         plot(f,t_fft{j+1,i+1})
%         xlim([0,60])
%         ylim([0,fft_ylim])
%         xlabel('hz')
%         title(sprintf('%s Ch%d U%d %d-%ds %s',label2,NEVchannels(i),NEVunits(i),L2(j,1),L2(j,2),NEV.MetaTags.Filename))
%         
%         saveas(gcf,[subpath sprintf('%s_Ch%d_U%d_%d-%ds',label2,NEVchannels(i),NEVunits(i),L2(j,1),L2(j,2)),'.fig'])
%         saveas(gcf,[subpath sprintf('%s_Ch%d_U%d_%d-%ds',label2,NEVchannels(i),NEVunits(i),L2(j,1),L2(j,2)),'.jpg'])
% 
%     end
% end
% 
% cd(GetFilePath)
% save([sprintf('Combine_nt%d_t%d_u%d_xcorr',length(L1),length(L2),length(NEVunits)) '.mat'],'index_minus','I_nt','I_t','index_nt','index_t','nt_fft','nt_corr_dsraw','nt_lags_dsraw','t_fft','t_corr_dsraw','t_lags_dsraw','tremor','no_tremor')