clear all
kkk1=0;
aa=0;
while aa==0
    clear FFinv FFinv1 clustermitos1 Ensemble signal3 FFmean FFstd Ensemble_peak peakcutoff cmmi cmma
    clear data data2 Ensemble signal2 totpeaks
    clear muf_rw ZZx ZZy L data data2 Maxhist  dminus dplus
    kkk1=kkk1+1;
    prompt={'Proceed with cell ' };
defans={num2str(kkk1)};
fields = {'one'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt, 'Go to', 1, defans,options);
if ~isempty(info)              %see if user hit cancel
   info = cell2struct(info,fields);  
[filename,pathname]= uigetfile(['C:\CVRC\MitochondrialMetabolism\*.mat'],'Choose the _frequ_dist file:');
 if exist([pathname filename])==2 % check whether there is a file
  fid=fopen([pathname filename]); % open file
    position=0;
    while ~feof(fid);
    data=load([pathname filename],'params','mitosincluster');
    file{kkk1}=filename(1:6);
    fclose(fid);break;
   end;
    else
        display('You suck, take an existing file!')
 end
 rw=data.params.rw;
 fs=data.params.fs;
 bb{kkk1}=data.params.bb;
 bbend{kkk1}=data.params.bbend;
 mitosincluster{kkk1}=data.mitosincluster;
 clear data
 prompt={'Name following file: '};
defans={filename(1:10)};
fields = {'name'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt, 'Mito#: ', 1, defans,options);
if ~isempty(info)              %see if user hit cancel
   info = cell2struct(info,fields);
   filename3 = info.name;                               
end
k1=1;

[filename1,pathname1]= uigetfile(pathname,'Choose the _wavefrequs_final:');
 if exist([pathname1 filename1])==2 % check whether there is a file
  fid=fopen([pathname1 filename1]); % open file
    position=0;
    while ~feof(fid);
    data=load([pathname1 filename1]);
    fclose(fid);break;
   end;
    else
        display('You suck, take an existing file!')
 end
[filename1,pathname1]= uigetfile(pathname,'Choose the _mat:');
 if exist([pathname1 filename1])==2 % check whether there is a file
  fid=fopen([pathname1 filename1]); % open file
    position=0;
    while ~feof(fid);
    data2=load([pathname1 filename1]);
    fclose(fid);break;
   end;
    else
        display('You suck, take an existing file!')
 end

 Ntott{kkk1,k1}=length(data.FF);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%take mat-file and determine frequency clusters for individual myocytes%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear ZZ
prompt={'Sampling frequency [sec]','Auto CC cutoff for maxpeak','rwsize','Spacing','Peak cutoff','CC cutoff','Alternative CC cutoff', 'Alternative2 CC cutoff','Intervall around peak [mHz]:','Peak-Intervall cutoff (innercutoff):','Normalized Mito cutoff: '};
defans={num2str((1./fs)*1000),num2str(0.5),num2str(rw),num2str(0.1),num2str(0.05),num2str(0.95),num2str(0.95),num2str(0.95),num2str(0.5),num2str(0.5),num2str(0.05)};
fields = {'fs','aCC','rw','name','name2','name3','name4','name5','roundseg','innercutoff','mitocutoff'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt, 'Define parameters: ', 1,defans, options);
if ~isempty(info)              %see if user hit cancel
    info = cell2struct(info,fields);
    seg=str2num(info.name);
    peakcutoff1=str2num(info.name2);
    ccut=str2num(info.name3);
    cccut=str2num(info.name4);
    cccut_alt2=str2num(info.name5);
    rw=str2num(info.rw);
    auto_CC_cutoff=str2num(info.aCC);
    fs=1/str2num(info.fs).*1000;                %%CAVE: fs in [mHz]!!
    roundseg=str2num(info.roundseg);
    innercutoff=str2num(info.innercutoff);
    mito_cutoff=str2num(info.mitocutoff);
    
end
[IND,num]=size(data.FF);
for k1=1:num
prompt={'Frames to initiation: ','Stop-frames: ','Upper-frequ-cutoff: ','Lower-frequ-cutoff: '};
defans={num2str(bb{kkk1}{k1}),num2str(bbend{kkk1}{k1}),...
    num2str(round2(fs./data.params.upper_cutoff{k1},log10(1/data.params.interpol_seg))),...
    num2str(round2(fs./(data.params.lower_cutoff{k1}),log10(1/data.params.interpol_seg)))};
fields = {'name','name2','name3','name4'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt, 'Frames: ', 1, defans,options);
if ~isempty(info)              %see if user hit cancel
   info = cell2struct(info,fields);
   bb{kkk1}{k1} = str2num(info.name); 
   bbend{kkk1}{k1}=str2num(info.name2);
   frequUcut{k1}=str2num(info.name3);
   frequLcut{k1}=str2num(info.name4);
end
end

% seg=0.05;       %take spacing as scale resolution dj (see wavelet)
% peakcutoff=0.05 %%cut off of bigger peaks


for k1=1:num
    ll=frequLcut{k1}:data.params.interpol_seg:frequUcut{k1};
vec=round(ll./seg);
    for a=1:length(data.FF{2,1})
        a

       for i=1:length(data.FF)
           if isempty(data.FF{i,k1})==0
           muf_rw{a,k1,1}(i)=data.FF{i,k1}(a);
           end
       end
   
% vec=data.params.frequvector{k1};
ccc=round(max(muf_rw{a,k1,1}(:))./seg);


if ccc<vec(end)-seg             %%make algorithm faster by not taking all frequ-scale points
    ZZx{a,k1}=vec(1:find(vec==ccc+1)).*seg;
else
    ZZx{a,k1}=vec.*seg;
end
ZZZx{a,k1,kkk1}=ZZx{a,k1};
% ZZx{a,k1}=interp1(x,y,xx);                  %%interpolated x-frequency-scale
ZZy{a,k1}=hist(muf_rw{a,k1,1},ZZx{a,k1});   %%number of corresponding mitos
ZZZy{a,k1,kkk1}=ZZy{a,k1};
end
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%standard deviation and mean %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k1=1:num 
for k=1:IND
     for i=1:length(data.FF{2,1})
         if isempty(data.FF{k,k1})==0
             FFinv{i,k1}(k)=data.FF{k,k1}(i);
         else
             FFinv{i,k1}(k)=0;
         end
     end
end
end
    
for k1=1:num
 for i=1:length(data.FF{2,1})
     c=find(FFinv{i,k1}>0);
     FFinv1{i,k1}=FFinv{i,k1}(c);
     cc=find(FFinv1{i,k1}>ZZx{i,k1}(1) & FFinv1{i,k1}<ZZx{i,k1}(end));
     FFmean{k1}(i)=mean(FFinv1{i,k1}(cc));
     FFstd{k1}(i)=std(FFinv1{i,k1});
 end
end
 
% %--------------------------------------------------------------------------
% 
% %--------------------------------------------------------------------------
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%now classify separated peaks and order them: %%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for k1=1:num
    for a=1:length(data.FF{2,1})
         a
%         peakcutoff{a,k1}=0.1.*;
        for kk=1:length(ZZx{a,k1})                      %%CAVE: cut-off still in scale
            Ensemble.bars{a,kk,k1}=...
                find(muf_rw{a,k1}>=ZZx{a,k1}(kk)-seg./2 & ...
                muf_rw{a,k1}<=ZZx{a,k1}(kk)+seg./2);        %%giving mitos in each bar for each RW-frame
        end
        
        Maxhist{a,k1}=min(find(ZZy{a,k1}(2:end-1)==max(ZZy{a,k1}(2:end-1))))+1;   %find maximum of ZZy frequ distribution
%         peakcutoff{a,k1}=FFmean{k1}(a)+FFstd{k1}(a);
        peakcutoff{a,k1}=0.05.*ZZy{a,k1}(Maxhist{a,k1});
%         clear p ZZy_n
%         p{k1}(1)=1;    %set correlation coeff 0
%         p{k1}(2)=0;
%         ZZy_n{a,k1}=ZZy{a,k1};
%         Maxhistt2=min(find(ZZy_n{a,k1}(2:end-1)==min(ZZy_n{a,k1}(2:end-1))))+1;  %find minimal point in Distr. as start for while loop  
% %         
%         if a<rw/2
%             Maxhist{a,k1}=min(find(ZZy_n{a,k1}(2:end)==max(ZZy_n{a,k1}(2:end-1))))+1;
%             
%         elseif a>=rw/2 & a<=length(data.gpmean{2,1})-rw/2
%             z=0;
%         for k=1:2
%                 if p{k1}(2)<=auto_CC_cutoff      %%start looking for first peak, that has Auto-CC > auto_CC_cutoff (0.75)
%                      p{k1}(1)=p{k1}(2);
%                     z=z+1;
%                     ZZy_n{a,k1}(Maxhistt2)=0;
%                     Maxhistt=Maxhistt2;
%                     Maxhistt2=min(find(ZZy_n{a,k1}(2:end-1)==max(ZZy_n{a,k1}(2:end-1))))+1;           %%determine max-frequ
%                     if length(Ensemble.bars{a,Maxhistt2,k1})<=innercutoff.*length(Ensemble.bars{a,Maxhistt,k1})
%                         p{k1}(2)=0;
%                     else
%                     c=Ensemble.bars{a,Maxhistt2,k1};            
%                     ll=(max(a-rw./2+1,1):min(length(data.gpmean{2,1}),a+rw./2));    %set window rw for inner CC
%                     Z{a}=zeros(1,length(ll));
%                     for kk=1:length(c)
%                         Z{a}=Z{a}+data.gpmean{c(kk)+1,k1}(ll);        %%take mean of tmre signal of 
%                     end                                                     %% tmRe signal, the "+1" stems from the fact that the first entry in gpmean is 0 (compare wavelet_analysis_v2)
%                     Z{a}=Z{a}./length(c);
%                     LUz=median(Z{a});
%                     R{a,k1}=zeros(2,2);
%                     for kk=1:length(c)
%                             LUm=median(data.gpmean{c(kk)+1,k1}(ll));
%                         CORRCOEF{kk,k1}=...                                    %%do CC
%                         corrcoef(Z{a}-LUz,data.gpmean{c(kk)+1,k1}(ll)-LUm);%%get corr coeff and check in while loop if it satisfies the auto_CC_cutoff
%                         R{a,k1}=R{a,k1}+CORRCOEF{kk,k1};
%                     end
%                     R{a,k1}=R{a,k1}./length(c);
%                     p{k1}(2)=R{a,k1}(2,1); 
%                     end
%                 end
%         end
%          
%             if z==2 & p{k1}(2)<=p{k1}(1) 
%             Maxhist{a,k1}=Maxhistt;
%             elseif z==2 & p{k1}(2)>p{k1}(1) & ...
%                     length(Ensemble.bars{a,Maxhistt,k1})>=innercutoff.*length(Ensemble.bars{a,Maxhistt,k1})
%             Maxhist{a,k1}=Maxhistt2;
%             else
%             Maxhist{a,k1}=Maxhistt2;
%             end
%         
%         elseif a>length(data.gpmean{2,1})-rw/2
%             Maxhist{a,k1}=min(find(ZZy_n{a,k1}(2:end-1)==max(ZZy_n{a,k1}(2:end-1))))+1;
%         end


%             for kk=max(Maxhist{a,k1}-(roundseg)./seg,2):1:min(Maxhist{a,k1}+(roundseg)./seg,length(ZZx{a,k1})-1)        %%check in 1mHz interval around maxpeak if there are similar high peaks and take them in
%                 if length(Ensemble.bars{a,kk,k1})>=innercutoff.*length(Ensemble.bars{a,Maxhist{a,k1},k1})
%                     test(kk-max(Maxhist{a,k1}-(roundseg)./seg,2)+1)=1;
%                 else
%                     test(kk-max(Maxhist{a,k1}-(roundseg)./seg,2)+1)=0;
%                 end
%             end
%             ccmin=min(min(find(test>0))+max(Maxhist{a,k1}-(roundseg)./seg,2)-1,Maxhist{a,k1});
%             ccmax=max(max(find(test>0))+max(Maxhist{a,k1}-(roundseg)./seg,2)-1,Maxhist{a,k1});
            
            
            
            ccmin=max(Maxhist{a,k1}-(roundseg)./seg,2);
            ccmax=min(Maxhist{a,k1}+roundseg./seg,length(ZZx{a,k1})-1);
        
                
        Ensemble_peak{a,k1}(1)=0;           %%throw away cut-off frequ
        for kk=2:length(ZZx{a,k1})-1
            if ZZy{a,k1}(kk)>=peakcutoff{a,k1} %.*ZZy{a,k1}(Maxhist{a,k1}) 
                 Ensemble_peak{a,k1}(kk)=2;                 %%determ if bar significant                        
            elseif ZZy{a,k1}(kk)<peakcutoff{a,k1}  & ZZy{a,k1}(kk)>0   %peakcutoff.*ZZy{a,k1}(Maxhist{a,k1})
                Ensemble_peak{a,k1}(kk)=1;
            else
                Ensemble_peak{a,k1}(kk)=0;
            end
        end
        Ensemble_peak{a,k1}(end+1)=0;
        c10=find(Ensemble_peak{a,k1}>0);
        c11=Ensemble_peak{a,k1}(c10);
        
        kk=max(find(c10<Maxhist{a,k1}));
        while c11(kk)==2 & kk>1
            kk=kk-1;
        end
        c1=max(2,c10(kk+1));
        
        
        kk=min(find(c10>Maxhist{a,k1}));
        while c11(kk)==2 & kk<length(c11)
            kk=kk+1;
        end
        c2=min(length(ZZx{a,k1})-1,c10(kk-1));
        Ensemble_peak{a,k1}(c1:c2)=2;
        clear c10
        c10=find(Ensemble_peak{a,k1}==2);
        c12=find(Ensemble_peak{a,k1}==1);
        Ensemble_peak{a,k1}(c12)=0;
        Ensemble_peak{a,k1}(c10)=1;
        
        
        [Ensemble.peak_L{a,k1},Ensemble.peak_no{a,k1}]=bwlabel(Ensemble_peak{a,k1});        %%every peak with corrseponding bars. 
       
    end    
        clear ZZy_n
end



for k1=1:num                                    %%determine the mean signal Ensemble.signal for each max peak
    for a=1:length(data.FF{2,1})
        a
        ll=(max(a-rw./2+1,1):min(length(data.gpmean{2,1}),a+rw./2));
     
        for k=1:Ensemble.peak_no{a,k1}
            c=find(Ensemble.peak_L{a,k1}==k);           %%for each coherent cluster peak find corresponding histogram bars
            Ensemble.signal{a,k1,k}=zeros(1,length(ll));
            z=0;
            for kk=1:length(c)
                cc=Ensemble.bars{a,c(kk),k1};       %%find mitos in the bars for whole cluster peak
                if Maxhist{a,k1}==c(kk)
                    Ensemble.maxpeak{a,k1}=k;       %%find maximum cluster peak
                end
                for kkk=1:length(cc)
                    z=z+1;
%                 Ensemble.signal{a,k1,k}=Ensemble.signal{a,k1,k}+smooth(data.gpmean{cc(kkk)+1}(a:a+127),80,'rloess')';
                Ensemble.signal{a,k1,k}=Ensemble.signal{a,k1,k}+data.gpmean{cc(kkk)+1,k1}(ll);  %%get mean signal for all cluster peaks
                end
            end
            
            Ensemble.signal{a,k1,k}=Ensemble.signal{a,k1,k}./z;     %divide by Number of Mitos in cluster peak
            L{a,k1,k}=z;                                            %%number of mitos in clusterpeak k
        end
        
        for k=1:length(ZZx{a,k1})
            Ensemble.signal2{a,k1,k}=zeros(1,length(ll));
            z=0;
            for kk=1:length(Ensemble.bars{a,k,k1})
                Ensemble.signal2{a,k1,k}=Ensemble.signal2{a,k1,k}+data.gpmean{Ensemble.bars{a,k,k1}(kk)+1,k1}(ll);
                z=z+1;
            end
            if z~=0
            Ensemble.signal2{a,k1,k}=Ensemble.signal2{a,k1,k}./z;
            end
        end
        
        [maxtab,mintab]=peakdet(ZZy{a,k1},0.5);                 %%find max and min of histogram
        ss=interp1(maxtab(:,1),maxtab(:,2),1:length(ZZx{a,k1}));    %%interpolate between max over range of histogram
        Diff=diff(ss);                                              %%find difference function
        pCl=1;
        MM=min(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1})); 
        k=max(MM-2,2);                  %beginning of maximum cluster peak
        MaP=Ensemble.maxpeak{a,k1};                             
        while pCl>=ccut & k>1
            k=max(MM-2,2);
            while k>1 & Diff(k)>=0 
                k=k-1;
            end
            Minbreakleft=k;                             %%find first Min. of interpolated maximum function to the left
            while Ensemble_peak{a,k1}(k)==0 & k>=2
                k=k-1;
            end                                         %%find first peak after Min to the left that is above cutoff (i.e. belonging to a cluster peak)
            eins=k;
            k=Minbreakleft;
            while  k>1 & Diff(k)<=0                       %%find first Max after Min to the left
                k=k-1;
            end
            if ZZy{a,k1}(k)>ZZy{a,k1}(k+1)
                zwei=k;
            else
                zwei=k+1;
            end
            
            
          
            if eins==1 & zwei==1
                Nextsignifpeakleft=2;
                pCl=0;
                Peaknol=1;
                Minbreakleft=2;
            else
                Nextsignifpeakleft=max(eins,zwei);           %%take Max. of left max and first peak above cutoff and check peak

%             clear Peaknol
            Peaknol=-1;
            for k=1:Ensemble.peak_no{a,k1}
                c=find(Ensemble.peak_L{a,k1}==k);
                for kk=1:length(c)
                    if c(kk)==Nextsignifpeakleft
                        Peaknol=k;                      %%find corresponding cluster peak
                    end
                end
            end

            EE=zeros(1,length(ll));
            zz=0;
            if Peaknol~=-1 & Peaknol<MaP
                for k=Peaknol+1:MaP
                    EE=EE+Ensemble.signal{a,k1,k}.*L{a,k1,k};
                    zz=zz+L{a,k1,k};                            %%get mean signal of corresponding cluster peak
                end
            else 
                EE=Ensemble.signal2{a,k1,Nextsignifpeakleft}./length(Ensemble.bars{a,Nextsignifpeakleft,k1});
            end
                
            if MaP>Peaknol & Peaknol~=-1
                EE=EE./zz;
                Sigmax=EE;
                Sigleft=Ensemble.signal{a,k1,Peaknol};
                LUmax=median(Sigmax);
                LUl=median(Sigleft);
                CORRCOEFl=corrcoef(Sigmax-LUmax,Sigleft-LUl);
                pCl=CORRCOEFl(2,1);                             %%check if corresponding cluster peak is highly correlated to maximum cluster peak
                MM=min(find(Ensemble.peak_L{a,k1}==Peaknol));
                MaP=Peaknol;
            elseif Peaknol==-1
                EEmax=zeros(1,length(ll));
                z=0;
                            for k=Minbreakleft:Maxhist{a,k1}
                                EEmax=EEmax+Ensemble.signal2{a,k1,k};
                                    z=z+1;
                            end
                        EEmax=EEmax./z;
                Sigmax=EEmax;
                Sigleft=EE;
                LUmax=median(Sigmax);
                LUl=median(Sigleft);
                CORRCOEFl=corrcoef(Sigmax-LUmax,Sigleft-LUl);
                pCl=CORRCOEFl(2,1);                             %%check if corresponding cluster peak is highly correlated to maximum cluster peak
                c=find(ZZy{a,k1}>0);
                k=Nextsignifpeakleft-1;
                while k>1 & ZZy{a,k1}(k)==0 
                    k=k-1;
                end
                MM=k;
%                 MaP=Nextsignifpeakleft;
            else
                pCl=0;
                Peaknol=Peaknol-1;
            end
            end
        end
        if Peaknol==1 & Minbreakleft<=min(find(Ensemble.peak_L{a,k1}==1))       %%determine start-bar in Histogram by also considering border conditions
            cl=Minbreakleft;         
        elseif Peaknol==Ensemble.peak_no{a,k1}
            cl=min(find(Ensemble.peak_L{a,k1}==Peaknol));
         k=1;
         while cl>=mintab(k,1)
            k=k+1;
         end
         cl=max(2,mintab(k-1,1));
        elseif Peaknol==-1
            cl=Minbreakleft;
        else
         Peaknol=Peaknol+1;
         cl=min(Minbreakleft,min(find(Ensemble.peak_L{a,k1}==Peaknol)));
         k=1;
         if cl>mintab(1,1)
         while cl>=mintab(k,1)
            k=k+1;
         end
         if k==1
             cl=2;
         else
         cl=max(2,mintab(k-1,1));
         end
         end
        end  
        
        
        %%%%%%%%%%%%%RIght hand side%%%%%%%%%%%%%%%%%%%%%
         pCr=1;                                                                     %%now start the same for the right hand side of the histogram
        MM=max(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1})); 
        MaP=Ensemble.maxpeak{a,k1};
        k=min(MM+2,length(ZZx{a,k1})-1);
        while pCr>=ccut & k<length(ZZx{a,k1})
            k=min(MM+2,length(ZZx{a,k1})-1);
            while k<length(ZZx{a,k1}) & Diff(k)<=0 
                k=k+1;
            end
            Minbreakright=k;
            while k<=length(ZZx{a,k1})-1 & Ensemble_peak{a,k1}(k)==0 
                k=k+1;
            end
            eins=k;
            k=Minbreakright;
            while  k<length(ZZx{a,k1})-1 & Diff(k)>=0                    %%find first Max after Min to the left
                k=k+1;
            end
            if k==length(ZZx{a,k1})-1
                zwei=k+1;
            else
                if ZZy{a,k1}(k+1)>ZZy{a,k1}(k)
                zwei=k+1;
                else
                    zwei=k;
                end
            end
            
            if eins==length(ZZx{a,k1}) & zwei==length(ZZx{a,k1})
            Nextsignifpeakright=length(ZZx{a,k1})-1;
            Peaknor=Ensemble.peak_no{a,k1};
            pCr=0;
            Minbreakright=length(ZZx{a,k1})-1;
            else
          
                Nextsignifpeakright=min(eins,zwei);
                Peaknor=-1;
            for k=1:Ensemble.peak_no{a,k1}
                c=find(Ensemble.peak_L{a,k1}==k);
                for kk=1:length(c)
                    if c(kk)==Nextsignifpeakright
                        Peaknor=k;
                    end
                end
            end

            
            EE=zeros(1,length(ll));
            zz=0;
            if Peaknor>MaP & Peaknor~=-1
            for k=MaP:Peaknor-1
            EE=EE+Ensemble.signal{a,k1,k}.*L{a,k1,k};
            zz=zz+L{a,k1,k};
            end
            else
              EE=Ensemble.signal2{a,k1,Nextsignifpeakright}./length(Ensemble.bars{a,Nextsignifpeakright,k1});
            end
            if Peaknor>MaP & Peaknor~=-1
            EE=EE./zz;
            Sigmax=EE;
            Sigright=Ensemble.signal{a,k1,Peaknor};
            LUmax=median(Sigmax);
            LUr=median(Sigright);
            CORRCOEFr=corrcoef(Sigmax-LUmax,Sigright-LUr);
            pCr=CORRCOEFr(2,1); 
            MM=max(find(Ensemble.peak_L{a,k1}==Peaknor));
            MaP=Peaknor;
            elseif Peaknor==-1
                EEmax=zeros(1,length(ll));
                z=0;
                    for k=Maxhist{a,k1}:Minbreakright
                        EEmax=EEmax+Ensemble.signal2{a,k1,k};
                        z=z+1;
                    end
                    EEmax=EEmax./z;
                Sigmax=EEmax;
                Sigright=EE;
                LUmax=median(Sigmax);
                LUr=median(Sigright);
                CORRCOEFl=corrcoef(Sigmax-LUmax,Sigright-LUr);
                pCr=CORRCOEFl(2,1);                             %%check if corresponding cluster peak is highly correlated to maximum cluster peak
                c=find(ZZy{a,k1}>0);
                k=Nextsignifpeakright+1;
                while k<length(ZZx{a,k1}) & ZZy{a,k1}(k)==0 
                    k=k+1;
                end
                MM=k;
%                 MaP=Nextsignifpeakright;
            else
                pCr=0;
                Peaknor=Peaknor+1;
            end
            end
        end
        if Peaknor==Ensemble.peak_no{a,k1} & Minbreakright>=max(find(Ensemble.peak_L{a,k1}==Ensemble.peak_no{a,k1}))
            cr=Minbreakright; 
        elseif Peaknor==1
            cr=max(find(Ensemble.peak_L{a,k1}==Peaknor));
         k=length(mintab(:,1));
         while cr<=mintab(k,1)
            k=k-1;
         end
         cr=min(length(mintab(:,1)),mintab(k+1,1));
        elseif Peaknor==-1
            cr=Minbreakright;
        else
         Peaknor=Peaknor-1;
%          if Peaknor==Ensemble.maxpeak{a,k1}
             
         cr=max(find(Ensemble.peak_L{a,k1}==Peaknor));
         k=length(mintab(:,1));
         if cr<mintab(end,1)
         
         while cr>=mintab(k,1) & k>2
            k=k-1;
         end
         if k==length(mintab(:,1))
             cr=length(ZZx{a,k1})-1;
         else
         cr=max(Minbreakright,min(length(mintab(:,1)),mintab(k+1,1)));
         end
         end
        end   
            
            
        
    cmmi{k1}(a)=cl;
    cmma{k1}(a)=cr;
        
        
        
        
    end
end




clear rightpeaks leftpeaks clusterbars  CORRCOEFl CORRCOEFr Ensemble.signal
% 
%%------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Alternative clusterbar selection by only taking those in general with a
%%%CC>0.8%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k1=1:num
    for a=1:length(data.FF{2,1})
        a
        ll=(max(a-rw./2+1,1):min(length(data.gpmean{2,1}),a+rw./2));
        
         %%now the second alternative version: take signals from each
            %%ind. mito:
            for k=1:length(data.FF)
                if isempty(data.FF{k,k1})==0
                signal3{k,k1}=data.gpmean{k+1,k1}(ll);
                end
            end
            
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%now do the CC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             cc1(1)=min(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1}));
%             cc1(2)=max(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1}));
            
            
%         clustermaxmin=min(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1}));
%         clustermaxmax=max(find(Ensemble.peak_L{a,k1}==Ensemble.maxpeak{a,k1}));
%         clustermax=clustermaxmin:1:clustermaxmax;
        clustermax=cmmi{k1}(a):cmma{k1}(a);
%         clear clustermaxmitos
        clustermaxmitos{a,k1}=Ensemble.bars{a,clustermax(1),k1};
        if length(clustermax)>=2
          for kk=2:length(clustermax)
              clustermaxmitos{a,k1}=[clustermaxmitos{a,k1},Ensemble.bars{a,clustermax(kk),k1}];
          end
        end
        
 
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           %%%this using the second alternative version%%%%%%%%%%%%%%%%
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
           for k=1:length(data.FF)
               if isempty(data.FF{k,k1})==0
               if isempty(find(clustermaxmitos{a,k1}==k))==1
                   clear Sigmax
                   Sigmax=Ensemble.signal{a,k1,Ensemble.maxpeak{a,k1}};
                   LUmax=median(Sigmax);
                   Sigtot=signal3{k,k1};
                   LUtot=median(Sigtot);
                   CORRCOEFtot=...
                       corrcoef(Sigmax-LUmax,...
                       Sigtot-LUtot);
                   pCtot=CORRCOEFtot(2,1);
                if pCtot>=cccut_alt2
                totpeaks{a,k1}(k)=1;
                else
                totpeaks{a,k1}(k)=0;
                end
               else
                   totpeaks{a,k1}(k)=1;
               end
               end
           end
        clear cctot
        cctot=find(totpeaks{a,k1}==1);
        

%--------------------------------------------------------------------------
          %%%%determine max and min frequencies within the range of 1 STD
          %%%%for all mitos that are >=cccut_alt2 correlated to
          %%%%the max peak
%           clear dd
%           for i=1:length(cctot)
%               dd(i)=data.FF{cctot(i),k1}(a);
%           end
%           Ensemble.maxmean{k1}(a)=mean(dd(:));
%           stdd=std(dd(:));
% %           maxx{a,k1}=max(dd(:));
% %           minn{a,k1}=min(dd(:));
%           Ensemble.maxfrequ{k1}(a)=max(ZZx{a,k1}(2),Ensemble.maxmean{k1}(a)-stdd);
%           Ensemble.minfrequ{k1}(a)=min(ZZx{a,k1}(end-2),Ensemble.maxmean{k1}(a)+stdd);
          dplus{a,k1}=cctot;
          for i=1:length(cctot)
              if FFinv{a,k1}(cctot(i))==min(ZZx{a,k1}(:)) | FFinv{a,k1}(cctot(i))==max(ZZx{a,k1}(:))
                  cctot(i)=0;
              end
          end
        dminus{a,k1}=cctot(find(cctot>0));
        
        if isempty(dminus{a,k1})==1
%           ctot{a,1,k1}=0;
%           clustertot=0;
            Ensemble.clustermitos_alternative2{a,k1}=0;
        else
             Ensemble.clustermitos_alternative2{a,k1}=dminus{a,k1};
        end 
              
          
          
                
    end
end
clear clustermitos1
clustermitos1=Ensemble.clustermitos_alternative2;
clustermitos1=dplus;
clear CORRCOEFl CORRCOEFr CORRCOEFtot pCr pCl pCtot signal2 signal3
% clear rightpeaks leftpeaks totpeaks
version=3;
clear mitosinclusterr
for k1=1:num
    for a=1:length(data.FF{2,1})
        a
%         mito_no{version,k1}(a)=length(find(Bwcluster{a,k1}==1))./length(find(data.bww1{k1}));
        mitosinclusterr{version,k1}(a)=length(clustermitos1{a,k1})./length(data.FF);
        
    end
end

     cmmii{k1,kkk1}=cmmi{k1};
     cmmaa{k1,kkk1}=cmma{k1};
     

 clear signal3 signal signal2 Ensemble
else
    aa=1;
end
end
kkk1=kkk1-1;
L=kkk1;


%%%CAVE: ONLY CONSIDERING k1 = 1!!


lCUT=0;                 %lower frequ cut-off
uCUT=100;                %upper frequ cut-off
for kkk1=1:L
    kkk1
    for k1=1:1
        z1{kkk1,k1}=0;  
        Meanclusterhisty{k1,kkk1}=zeros(1,length(lCUT:seg:uCUT));
        Nsingle{k1,kkk1}=0;
        for a=1:length(cmmii{k1,kkk1})
            
            if a>=bb{kkk1}{k1} & a<=bbend{kkk1}{k1}
                
            clear BBvec
            BBvec=max(lCUT,ZZZx{a,k1,kkk1}(1)):seg:min(ZZZx{a,k1,kkk1}(end),uCUT);      %define frequ vector being the basis for all myocytes
            Clusterhisty{a,k1,kkk1}=zeros(1,length(lCUT:seg:uCUT));                     %predefine histogram y-scale
            for b=1:length(BBvec)
                if cmmaa{k1,kkk1}(a)>length(BBvec)                                  %%check that upper frequ in single histogram is within the limits of frequ vector
                    cmmaa{k1,kkk1}(a)=length(BBvec);
                end
                if cmmii{k1,kkk1}(a)>cmmaa{k1,kkk1}(a)                              %%check that upper frequ 
                    cmmii{k1,kkk1}(a)=cmmaa{k1,kkk1}(a);
                end
                if BBvec(b)>=ZZZx{a,k1,kkk1}(1) & BBvec(b)<=ZZZx{a,k1,kkk1}(end)
            if BBvec(b)<ZZZx{a,k1,kkk1}(cmmii{k1,kkk1}(a)) | BBvec(b)>ZZZx{a,k1,kkk1}(cmmaa{k1,kkk1}(a))
                    Clusterhisty{a,k1,kkk1}(round((BBvec(b)-lCUT)*(1/seg))+1)=0;
            else
                    Clusterhisty{a,k1,kkk1}(round((BBvec(b)-lCUT)*(1/seg))+1)=...
                        ZZZy{a,k1,kkk1}(round((BBvec(b)-ZZZx{a,k1,kkk1}(1))*(1/seg))+1);
            end
                end
            end
            Nsingle{k1,kkk1}=Nsingle{k1,kkk1}+mitosincluster{kkk1}{3,k1}(a).*Ntott{kkk1,k1};
            Meanclusterhisty{k1,kkk1}=Meanclusterhisty{k1,kkk1}+...
                Clusterhisty{a,k1,kkk1}./(mitosincluster{kkk1}{3,k1}(a).*Ntott{kkk1,k1});
            z1{kkk1,k1}=z1{kkk1,k1}+1;
            end
        end
        Meanclusterhisty{k1,kkk1}=Meanclusterhisty{k1,kkk1}./z1{kkk1,k1};
    end
end

Meanmeanclusterhisty=zeros(1,length(lCUT:seg:uCUT));
Meanmeanclusterhistx=lCUT:seg:uCUT;
Meanmeanclusterhistx=Meanmeanclusterhistx';
Ntot=0;
A(:,1)=Meanmeanclusterhistx;
for kkk1=1:L
    XX=Meanclusterhisty{1,kkk1}         %%./max(Meanclusterhisty{1,kkk1}(:))
    Meanmeanclusterhisty=Meanmeanclusterhisty+XX;
    Ntot=Ntot+Nsingle{1,kkk1};
    A(:,kkk1+1)=XX;
end

Meanmeanclusterhisty=Meanmeanclusterhisty';
figure; plot(Meanmeanclusterhisty)
%     
YY=Meanmeanclusterhisty;
XXX=Meanmeanclusterhistx;
YY = YY./max(YY);
ff=fittype('gauss1');
fopts=fitoptions('gauss1');
c = find(YY~=0);
xx = XXX(c);
yy = YY(c);
[cfun,gof,output]...
    =fit(xx,yy,ff);
% =fit(Meanmeanclusterhistx,Meanmeanclusterhisty,ff);
figure; 
hold on
plot(XXX,YY)
plot(XXX,cfun(XXX),'r')
GAUSS1 = cfun(XXX);
limits.Fmax=max(find(cfun(XXX)>=0.1.*max(cfun(XXX))));
limits.Fmin=min(find(cfun(XXX)>=0.1.*max(cfun(XXX))));
limits.Ftotmax=find(cfun(XXX)==max(cfun(XXX)));

prompt={['Save files as Overall_frequs']};
defans={pathname};
fields = {'name'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt,'File',1, defans,options);
if ~isempty(info)              %see if user hit cancel
    info = cell2struct(info,fields);
    ii=info.name;
save([ii,'Overall_frequs_WH' ],...
   'ZZZx','ZZZy','cmmii','cmmaa','Meanclusterhisty',...
   'cfun','gof','output','limits',...
   'Meanmeanclusterhisty','Meanmeanclusterhistx','file','bb','bbend','Nsingle','Ntot');
end  
prompt={['Save files as Overall_frequs_origin']};
defans={ii};
fields = {'name'};
options.Resize='on';
options.WindowStyle='normal';
info = inputdlg(prompt,'File',1, defans,options);
if ~isempty(info)              %see if user hit cancel
    info = cell2struct(info,fields);
    ii=info.name;
save([ii,'Overall_frequs_WH_origin' ],...
   'A','YY','GAUSS1');
end  





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% processing from 03/27/12
% YY = YY./max(YY);
% for k=1:length(XXX)
%     xwD(k) = YY(k)*XXX(k);
%     wD(k) = YY(k);
% end
% MD = sum(xwD)./sum(wD);   %%% mean of Dist
% MV = xwD./sum(wD);
% std(MV)
% %% PYR: 18.2924
% %% 
% sqrt((1./length(xwD))*sum((MV-MD).^2))
% a=0;
% for k=1:length(XXX)
%     if wD(k) > 0
%         a=a+1;
%         Dy(a) = YY(k);
%         Dx(a) = XXX(k);
%     end
% end
% 
% [cfun,gof,output]...
%     =fit(Dx',Dy',ff);
% GAUSS = cfun(XXX);
% figure; hold on; plot(XXX,YY); plot(XXX,GAUSS,'r')
% 
% 
% limits.Fmax=max(find(cfun(XXX)>=0.1.*max(cfun(XXX))));
% limits.Fmin=min(find(cfun(XXX)>=0.1.*max(cfun(XXX))));
% limits.Ftotmax=find(cfun(XXX)==max(cfun(XXX)));
% 
%         save([path,'Overall_frequs_lactate_origin' ],...
%    'A','YY','GAUSS');
% 
% save([path,'Overall_frequs_lactate' ],...
%    'ZZZx','ZZZy','cmmii','cmmaa','Meanclusterhisty',...
%    'cfun','gof','output','limits',...
%    'Meanmeanclusterhisty','Meanmeanclusterhistx','file','bb','bbend','Nsingle','Ntot');