% Phase spike phase analysis (spike-time vs waveform phase by hilbert)
% Execute 'A_NsxFilter_ALLchannels_v2.m' first if requires filtering of data
% v2: add polarity index (Check PolarityProfiles 2)
% for SW: continuous .ns6, able to choose time frame(one by one)
% for SW v2: create folders & polarity index(PI) ** need validation
%        v2.1: polarity index=abs(polarity sum)/length ** "abs" of complex sum
%              => NOT the problem (changed back to without abs in this step, cuz will be caculated at polarity profile 2 eventually)
%              normalized to the wrong length(normalized to all spikes, not spikes in that time frame)
%        v2.2: check NSx&NEV sampling rate
%        v2_3: 2022.01.17 add 'kilo' option, uses F_openKilosort2 func./ cd before loading filtered data
%        v2_3b: 2022.01.18 modified a bit for filtered data from batch processing (line 71)
%        v2_3b2: 2022.02.07 compatible with decimated data
%        v2_3b3: 2022.04.07 multiple time frame & label on fig
%        v2_3b4: 2022.04.24 output shuffled spike phase & save a profile compilation
%        v2_3b5: 2022.05.16 compatible with .nev file 
%                2022.05.18 discard frames with too few units (below a threshold)
%        v2_3b6: 2022.05.22 allow filtered data with reference subtracted

clear all, close all;
%% Input parameters %%
% LFP channel profiles
LFPChannels=[10,11,129]; % can be multiple channels
RefSubtracted=1; % include filtered data with Ref subtracted or not
RefSubCh=[201,202,203]; % starts with 201, pseudo channel name, for indexing the Ref pair 
LoadFilterFile=1; % if 1, loadfilterfile for the LFP waveform
kilo=1; % load kilosort data or not
boss=0; % load .nev from boss 
%%% timeframe input at line 75

% NEV channel profiles    
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
% Bin spike-to-phase profiles
BinNumber=36;
  
Alltimeframes=[177,252;503,538;68,119;409,455;565,621];% in secs  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
frametype=[2,2,1,1,1]; % 1: move/ 2: rest/ 3: baseline move/ 4: baseline rest

HistThr=150; % if there are fewer than 150 spikes for a unit in a certain time frame, it will be saved seperately and won't be added to the total compilation

%%%%%%%%% Put input parameters into structure %%%%%%%%%%%%%
save temp.mat;
a_InputParameters=load('temp.mat');
delete temp.mat;
%%  Main function
%% open NEV file %%
%openNEV;
if kilo==1
    NEV=F_openKilosort2;
    NEVchannels=ones(size(NEVunits));
elseif kilo==0
    if boss==0
      NEV=F_openNEX;
    elseif boss==1
      NEV=openNEV;
    end
end

%% Define workspace 
SpikeFreqCell=cell(6,length(NEVunits)+1);
s_SpikeFreqCell=cell(6,length(NEVunits)+1);
A=num2cell(NEVchannels);
B=num2cell(NEVunits);
SpikeFreqCell(1,:)=['NEVchannels',A];
SpikeFreqCell(2,:)=['NEVchannels',B];
SpikeFreqCell(3,1)={'Phase of SU'};
SpikeFreqCell(4,1)={'SpikeHistCounts'};
SpikeFreqCell(5,1)={'PolarityCounts'};
SpikeFreqCell(6,1)={'PolarityIndex'};

s_SpikeFreqCell(1,:)=['NEVchannels',A];
s_SpikeFreqCell(2,:)=['NEVchannels',B];
s_SpikeFreqCell(3,1)={'Phase of SU'};
s_SpikeFreqCell(4,1)={'SpikeHistCounts'};
s_SpikeFreqCell(5,1)={'PolarityCounts'};
s_SpikeFreqCell(6,1)={'PolarityIndex'};
clear A B;

names=["Channels";"Units";"LFPChannels";"PolarityCounts";"PolarityIndex";"Degree(-pi:pi)";"Degree(0:2pi)"];

Failed_timeframes={};
%% Load LFP file%%
if LoadFilterFile==1;
    [A,B]=uigetfile('*.mat');
    cd (B);
    load(A);
    a_InputParameters_LFPdataSource=[A,B];
    clear A B;
else
    NSx=openNSx;
    SourceFile=[NSx.MetaTags.Filename NSx.MetaTags.FileExt];
    data=double(NSx.Data');
    NSx=rmfield(NSx,'Data'); 
    fs=NSx.MetaTags.SamplingFreq;
    Title1=SourceFile;
end
a_InputParameters.SourceFile=SourceFile;
%A=find(ismember(NSx.MetaTags.ChannelID,LFPChannels))';

% modified for filtered data from batch
clear i
if RefSubtracted==0
    A=find(ismember(filteronly,LFPChannels));
    data=data(:,A);
else
    A=find(ismember(filteronly,LFPChannels));
    refA=find(ismember(str2double(withref_pseudo_ch(2:end,1)),RefSubCh))';
    A=[A,refA+length(filteronly)]
    data=data(:,A);
    LFPChannels=[LFPChannels RefSubCh];
end

% %% NSx SourceFile identification %%
%  % check if the .nev & .nsx file matched
% if NSx.MetaTags.DateTime==NEV.MetaTags.DateTime;
%     NEV.MetaTags.Filename=NSx.MetaTags.Filename;

if NEV.MetaTags.SampleRes ~= NSx.MetaTags.SamplingFreq
    error(sprintf('Warning: Different sampling rate, NSx:%d, NEV:%d',NSx.MetaTags.SamplingFreq,NEV.MetaTags.SampleRes));
end

% if NSx.MetaTags.FileExt=='.ns6';
% else % correct the error of NEV if not processed by .ns6 file.
%     NEV.MetaTags.SampleRes=NSx.MetaTags.SamplingFreq;
%     NEV.MetaTags.DataDurationSec=NEV.MetaTags.DataDuration/NEV.MetaTags.SampleRes;
% end
% else
%     error('NEV & NSx file mismatch');
% end

compilation{1,1}=SourceFile;
compilation{6,1}="Shuffled";
%% Hilbert transform of LFP data %%
tic
HilbertLFP=hilbert(data); % hilbert transformed each column as a dataset (except only one line).
shuffled_phase=HilbertLFP(randperm(length(HilbertLFP)),:);
toc

LegendMarker2='';
%% NEV units
[sz,~]=size(Alltimeframes);
for z=1:sz
    SpikeFreqCell=cell(6,length(NEVunits)+1);
    s_SpikeFreqCell=cell(6,length(NEVunits)+1);
    A=num2cell(NEVchannels);
    B=num2cell(NEVunits);
    SpikeFreqCell(1,:)=['NEVchannels',A];
    SpikeFreqCell(2,:)=['NEVchannels',B];
    SpikeFreqCell(3,1)={'Phase of SU'};
    SpikeFreqCell(4,1)={'SpikeHistCounts'};
    SpikeFreqCell(5,1)={'PolarityCounts'};
    SpikeFreqCell(6,1)={'PolarityIndex'};

    s_SpikeFreqCell(1,:)=['NEVchannels',A];
    s_SpikeFreqCell(2,:)=['NEVchannels',B];
    s_SpikeFreqCell(3,1)={'Phase of SU'};
    s_SpikeFreqCell(4,1)={'SpikeHistCounts'};
    s_SpikeFreqCell(5,1)={'PolarityCounts'};
    s_SpikeFreqCell(6,1)={'PolarityIndex'};
    clear A B;

    Failed_timeframes{z,2}=[];
    timeframe=Alltimeframes(z,:);
    if frametype(z)==1
        frametypetext='move';
    elseif frametype(z)==2
        frametypetext='rest';
    elseif frametype(z)==3
        frametypetext='baseline move';
    elseif frametype(z)==4
        frametypetext='baseline rest';
    end

    dataframe=timeframe*fs;
    AngleEdges=-pi:2*pi/BinNumber:pi; % create Bin edges for polar system
    for h=1:length(NEVunits);  
        SU_Name(h)={[LegendMarker2,' ',sprintf('Chan%d-U%d',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. 
        E=NEV.Data.Spikes.TimeStamp(C(D)); % find spikes that match the sorted channel-unit.
        if filtertype=="decimate_fir"
            E=round(E/(NSx.MetaTags.SamplingFreq/dsfs));
            display('**Downsampled Spikes frames');
        end
        spikesinframe=E(E>dataframe(1)&E<dataframe(2));

        %F=E/3  %%%nope, but how to transform a 30k data into 10k?
    %     disp(size(C));
    %     SUpoints(h,:)=zeros(1,NEV.MetaTags.DataDuration); % plot points of SU in the timeline, step 1: fill with zeros
    %     SUpoints(h,E)=1; % plot points of SU in the timeline
        %A=SpikePhase(:,:rad2deg(angle(HilbertLFP(E))); % fine [-180,180] angle each phase

        SpikePhase=angle(HilbertLFP(spikesinframe,:))+pi/2; % correct the -90 degree shift after hilbert transform
        SpikePhase=wrapToPi(SpikePhase); % wrap angle to [-pi pi] range.
        s_SpikePhase=angle(shuffled_phase(spikesinframe,:))+pi/2; % correct the -90 degree shift after hilbert transform
        s_SpikePhase=wrapToPi(s_SpikePhase); % wrap angle to [-pi pi] range.
        
        if length(spikesinframe)<HistThr
            Failed_timeframes{z,1}=[timeframe(1),timeframe(2)];
            Failed_timeframes{z,2}=vertcat(Failed_timeframes{z,2},[NEVchannels(h),NEVunits(h)]);
        end
   
        SpikeFreqCell{3,h+1}=SpikePhase; %Store raw angle data into cell arrays(because the array size of each SU is different, need to use cell)
        for j=length(LFPChannels):-1:1
            SpikeHistCounts(:,j)=histcounts(SpikePhase(:,j),AngleEdges);
        end
        SpikeFreqCell{4,h+1}=SpikeHistCounts; % save SpikeHistCount to cells
        PolaritySums=sum(exp(i*(SpikePhase))); % absolute vector sum of each angle (complex number)
        PolarityIndex=PolaritySums/length(spikesinframe); % normalized vector sum, should between [0,1]
        SpikeFreqCell{5,h+1}=PolaritySums; 
        SpikeFreqCell{6,h+1}=PolarityIndex; 

        % for shuffled data

        s_SpikeFreqCell{3,h+1}=s_SpikePhase; %Store raw angle data into cell arrays(because the array size of each SU is different, need to use cell)
        for j=length(LFPChannels):-1:1
            s_SpikeHistCounts(:,j)=histcounts(s_SpikePhase(:,j),AngleEdges);
        end
        s_SpikeFreqCell{4,h+1}=s_SpikeHistCounts; % save SpikeHistCount to cells
        s_PolaritySums=sum(exp(i*(s_SpikePhase))); % absolute vector sum of each angle (complex number)
        s_PolarityIndex=s_PolaritySums/length(spikesinframe); % normalized vector sum, should between [0,1]
        s_SpikeFreqCell{5,h+1}=s_PolaritySums; 
        s_SpikeFreqCell{6,h+1}=s_PolarityIndex; 
  
    end
    [rf,~]=size(Failed_timeframes{z,2});

    % Calculate Polarity count and Polarity index
    PolarityProfiles=zeros(5,length(NEVunits)*length(LFPChannels)); 
    T=length(LFPChannels);
    for k=1:1:length(NEVunits);
       PolarityProfiles(1,((k-1)*T+1):k*T)=NEVchannels(k);
       PolarityProfiles(2,((k-1)*T+1):k*T)=NEVunits(k);
       PolarityProfiles(3,((k-1)*T+1):k*T)=LFPChannels;
       PolarityProfiles(4,((k-1)*T+1):k*T)=SpikeFreqCell{5,k+1}; % Polarity counts 
       PolarityProfiles(5,((k-1)*T+1):k*T)=SpikeFreqCell{6,k+1}; % Polarity Index
    end
    PolarityProfile2=abs(PolarityProfiles);
    PolarityProfile2(6,:)=rad2deg(angle(PolarityProfiles(4,:))); % Polarity index: degree [-180,180]
    PolarityProfile2(7,:)=wrapTo360(PolarityProfile2(6,:)); % Polarity index: degree[0 360]
    PolarityProfile2=[names PolarityProfile2];
    
    % for shuffled data
    s_PolarityProfiles=zeros(5,length(NEVunits)*length(LFPChannels)); 
    T=length(LFPChannels);
    for k=1:1:length(NEVunits);
       s_PolarityProfiles(1,((k-1)*T+1):k*T)=NEVchannels(k);
       s_PolarityProfiles(2,((k-1)*T+1):k*T)=NEVunits(k);
       s_PolarityProfiles(3,((k-1)*T+1):k*T)=LFPChannels;
       s_PolarityProfiles(4,((k-1)*T+1):k*T)=s_SpikeFreqCell{5,k+1}; % Polarity counts 
       s_PolarityProfiles(5,((k-1)*T+1):k*T)=s_SpikeFreqCell{6,k+1}; % Polarity Index
    end
    s_PolarityProfile2=abs(s_PolarityProfiles);
    s_PolarityProfile2(6,:)=rad2deg(angle(s_PolarityProfiles(4,:))); % Polarity index: degree [-180,180]
    s_PolarityProfile2(7,:)=wrapTo360(s_PolarityProfile2(6,:)); % Polarity index: degree[0 360]
    s_PolarityProfile2=[names s_PolarityProfile2];
    
    %% Graphical illustration (polarhistogram)
    GetFilePath=[pwd '\']
    SubPath=[sprintf('%d-%d(%s)',timeframe(1),timeframe(2),frametypetext)]; 
    SubPath=[SubPath,'\']
    if kilo==1
        FullFilePath=[GetFilePath,'(kilo)SpikePhaseAnalysisTimeFrames(v2_3b6)\',SubPath];
    elseif kilo==0
        FullFilePath=[GetFilePath,'SpikePhaseAnalysisTimeFrames(v2_3b6)\',SubPath];
    end
    mkdir(FullFilePath)
    failedfilepath=[FullFilePath,'TooFewSpikes\'];
    mkdir(failedfilepath)

    K=0;
    order=1;
    fID=1;
    for h=1:1:length(NEVunits);
        if fID<=rf
            if NEVchannels(h)==Failed_timeframes{z,2}(fID,1)&NEVunits(h)==Failed_timeframes{z,2}(fID,2)
            fID=fID+1;
            SpikePhase=SpikeFreqCell{3,h+1};
              for j=1:1:length(LFPChannels);
                  K=K+1;
                  figure(K);
                  subplot(2,1,1);
                  H(K)=polarhistogram(SpikePhase(:,j),BinNumber);
                  title([sprintf('(NA)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


                  subplot(2,1,2);
                  max_lim=1;
                  x_fake=[0 max_lim 0 -max_lim];
                  y_fake=[max_lim 0 -max_lim 0];
                  h_fake=compass(x_fake,y_fake);
                  hold on;
                  compass(PolarityProfiles(5,order));
                  set(h_fake,'Visible','off');
                  title(sprintf('PI:%f',PolarityProfile2(5,order+1)));

                  hold off; 

                  saveas(K,[failedfilepath,sprintf('(NA)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
                  saveas(K,[failedfilepath,sprintf('(NA)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


                  order=order+1;
              end 
            
        else
            SpikePhase=SpikeFreqCell{3,h+1};
          for j=1:1:length(LFPChannels);
              K=K+1;
              figure(K);
              subplot(2,1,1);
              H(K)=polarhistogram(SpikePhase(:,j),BinNumber);
              title([sprintf('%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


              subplot(2,1,2);
              max_lim=1;
              x_fake=[0 max_lim 0 -max_lim];
              y_fake=[max_lim 0 -max_lim 0];
              h_fake=compass(x_fake,y_fake);
              hold on;
              compass(PolarityProfiles(5,order));
              set(h_fake,'Visible','off');
              title(sprintf('PI:%f',PolarityProfile2(5,order+1)));

              hold off; 

              saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
              saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


              order=order+1;

          end
            end 
        else
          SpikePhase=SpikeFreqCell{3,h+1};
          for j=1:1:length(LFPChannels);
              K=K+1;
              figure(K);
              subplot(2,1,1);
              H(K)=polarhistogram(SpikePhase(:,j),BinNumber);
              title([sprintf('%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


              subplot(2,1,2);
              max_lim=1;
              x_fake=[0 max_lim 0 -max_lim];
              y_fake=[max_lim 0 -max_lim 0];
              h_fake=compass(x_fake,y_fake);
              hold on;
              compass(PolarityProfiles(5,order));
              set(h_fake,'Visible','off');
              title(sprintf('PI:%f',PolarityProfile2(5,order+1)));

              hold off; 

              saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
              saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


              order=order+1;

        end
        end
    end
    for t=1:rf
        A=find(str2double(PolarityProfile2(1,:))==Failed_timeframes{z,2}(t,1));
        B=find(str2double(PolarityProfile2(2,A))==Failed_timeframes{z,2}(t,2));
        PolarityProfile2(:,A(B))=[];
        PolarityProfiles(:,A(B)-1)=[];
        A=find(cell2mat(SpikeFreqCell(1,2:end))==Failed_timeframes{z,2}(t,1))+1;
        B=find(cell2mat(SpikeFreqCell(2,A))==Failed_timeframes{z,2}(t,2));
        SpikeFreqCell(:,A(B))=[];
    end
    filename=sprintf('SpikePhaseAnalysis_%d-%ds_%s.mat',timeframe(1),timeframe(2),frametypetext);
    save([FullFilePath,filename],'SourceFile','SpikeFreqCell','PolarityProfiles','PolarityProfile2');
    
    
    FullFilePath=[FullFilePath,'shuffled\']
    mkdir(FullFilePath)
    failedfilepath=[FullFilePath,'TooFewSpikes\'];
    mkdir(failedfilepath)

    K=0;
    order=1;
    fID=1;
    for h=1:1:length(NEVunits);
        if fID<=rf
            if NEVchannels(h)==Failed_timeframes{z,2}(fID,1)&NEVunits(h)==Failed_timeframes{z,2}(fID,2)
            fID=fID+1;
            s_SpikePhase=s_SpikeFreqCell{3,h+1};
              for j=1:1:length(LFPChannels);
                  K=K+1;
                  figure(K);
                  subplot(2,1,1);
                  H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber);
                  title([sprintf('(NA)(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


                  subplot(2,1,2);
                  max_lim=1;
                  x_fake=[0 max_lim 0 -max_lim];
                  y_fake=[max_lim 0 -max_lim 0];
                  h_fake=compass(x_fake,y_fake);
                  hold on;
                  compass(s_PolarityProfiles(5,order));
                  set(h_fake,'Visible','off');
                  title(sprintf('PI:%f',s_PolarityProfile2(5,order+1)));

                  hold off; 

                  saveas(K,[failedfilepath,sprintf('(NA)(Shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
                  saveas(K,[failedfilepath,sprintf('(NA)(Shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


                  order=order+1;
              end 
            
        else
            s_SpikePhase=s_SpikeFreqCell{3,h+1};
          for j=1:1:length(LFPChannels);
              K=K+1;
              figure(K);
              subplot(2,1,1);
              H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber);
              title([sprintf('(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


              subplot(2,1,2);
              max_lim=1;
              x_fake=[0 max_lim 0 -max_lim];
              y_fake=[max_lim 0 -max_lim 0];
              h_fake=compass(x_fake,y_fake);
              hold on;
              compass(s_PolarityProfiles(5,order));
              set(h_fake,'Visible','off');
              title(sprintf('PI:%f',s_PolarityProfile2(5,order+1)));

              hold off; 

              saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
              saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


              order=order+1;

          end
            end 
        else
          s_SpikePhase=s_SpikeFreqCell{3,h+1};
          for j=1:1:length(LFPChannels);
              K=K+1;
              figure(K);
              subplot(2,1,1);
              H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber);
              title([sprintf('(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]);


              subplot(2,1,2);
              max_lim=1;
              x_fake=[0 max_lim 0 -max_lim];
              y_fake=[max_lim 0 -max_lim 0];
              h_fake=compass(x_fake,y_fake);
              hold on;
              compass(s_PolarityProfiles(5,order));
              set(h_fake,'Visible','off');
              title(sprintf('PI:%f',s_PolarityProfile2(5,order+1)));

              hold off; 

              saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);
              saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]);


              order=order+1;

        end
        end
    end
    for t=1:rf
        A=find(str2double(s_PolarityProfile2(1,:))==Failed_timeframes{z,2}(t,1));
        B=find(str2double(s_PolarityProfile2(2,A))==Failed_timeframes{z,2}(t,2));
        s_PolarityProfile2(:,A(B))=[];
        s_PolarityProfiles(:,A(B)-1)=[];
        A=find(cell2mat(s_SpikeFreqCell(1,2:end))==Failed_timeframes{z,2}(t,1))+1;
        B=find(cell2mat(s_SpikeFreqCell(2,A))==Failed_timeframes{z,2}(t,2));
        s_SpikeFreqCell(:,A(B))=[];
    end
    
    filename=sprintf('(shuffled)SpikePhaseAnalysis_%d-%ds_%s.mat',timeframe(1),timeframe(2),frametypetext);
    save([FullFilePath,filename],'SourceFile','s_SpikeFreqCell','s_PolarityProfiles','s_PolarityProfile2');
    

compilation{2,z}=timeframe;
compilation{3,z}=frametypetext;
compilation{4,z}=PolarityProfile2;
compilation{7,z}=s_PolarityProfile2;
end

if kilo==1
    FullFilePath=[GetFilePath,'(kilo)SpikePhaseAnalysisTimeFrames(v2_3b6)\'];
elseif kilo==0
    FullFilePath=[GetFilePath,'SpikePhaseAnalysisTimeFrames(v2_3b6)\'];
end
save([FullFilePath,'PolarityCompilation.mat'],'SourceFile','compilation');