function [omat1,omat2,prfl,sigf1,sigf2]=sift2dv1_2_peak(c1cs2,s1cs2,timetab,ncyc,...
    nzf1,nzf2,tailc,orgfid,...
    sw1,sw2,ph1,ph2,wndon,wndend1,...
    blkzero,thre,...
    obs1,obs2,cent1,cent2,sigf1,sigf2,pick)
% *** 2D SIFT cycle (noise and refrow not requried anymore) ***
% in=t1/t2 time, out=sifted t1/t2 time matrix.
% use siftv2 for 2D processing
% use proc1d for 1D processing
% Y. Matsuki, April, 2008
% usage:
% [omat1,omat2,prfl]=sift2dv1(c1cs2,s1cs2,timetab,ncyc,...
%     sigf1,sigf2,nzf1,tailc,orgfid,...
%     sw1,sw2,ph1,ph2,wndon,wndend1,...
%     blkzero,thre,pick)
% ========================================================
% ncyc: number of SIFT cycle
% timetab: vector specifying ACQUIRED time points
% sigf2: signal region in NMR F2
% sigf1: signal region in NMR F1
% nzf1: number of total point after z-filling in t1
% tailc: tail constraint to zero, feb0208. use =0.
% orgfid: clm vector taken from a full matrix
%       to see correlation between rmsd and dynamic range. feb0708
% sw1: spectral width of NMR 1st dim
% ph1: phase [p0 p1] for NMR 1st dim
% wndend: 0-1, point window function dies down to zero for NMR 1st dim.
% blkzero: 1 for block zeroing, 0 for thresholding using thre
% pick: row vector with F2 points used for noise calculation

% feb1709, tue, nzf2 added.
% aug1610, mon, nzf2 in fft, nt2org after ifft is introduced.

% sept0712, modified sift2dv1.m, use fzrpeak.

% nov2813, modified sift2dv1_peak.m 
% noise calculation was obsolete, changed according to sift2dv2_peak.m

% jan0814, plist is deleted from the argument list as it is not used.



% prfl=zeros(ncyc,6);
% mx=zeros(ncyc,1);
% mstd=zeros(ncyc,1);
% mdzsq=zeros(ncyc,1);

tic
    nt1org=length(c1cs2(1,:));
    nt2org=length(c1cs2(:,1));
    tic
    % modified aug1610, ym, to use z-fill in t2.
%     temp1=fft(c1cs2,[],1);
%     temp2=fft(s1cs2,[],1);
    temp1=fft(c1cs2,nzf2,1);
    temp2=fft(s1cs2,nzf2,1);
    toc
%     disp('FT2, done.');

%     temp1=mphasec(temp1,8012.8,[88 0],1);
%     temp2=mphasec(temp2,8012.8,[88 0],1);
    temp1=mphasec(temp1,sw2,ph2,1);
    temp2=mphasec(temp2,sw2,ph2,1);

    [temp1,temp2]=rrcplx(temp1,temp2);
%     disp('swap for FT1, done.');

    cycm1=temp1;
    cycm2=temp2;
    
%     sdymrg=0;
%     dymrg=0.1;
    
% % % obsolete (pick and db01 are generated below)
% % %     brt=brightpart(sigf1,nzf1);

% comment out, apr1008,mtk
%======== make vectors indicating bright and dark portion =======
%     nbf1=length(sigf1(:,1));
%     chk=zeros(nzf1,1);
%     for n=1:nzf1
%         tgl=0;
%         for m=1:nbf1
%             if (n>=sigf1(m,1) && n<=sigf1(m,2))
%                 tgl=1;
%             end
%         end
%     chk(n,1)=tgl;
%     end
% %==========================
% % % % % obsolete   
% % % % % make random coordinate>>>>>>>>
% % % % % npts=500;
% % % % npts=10; % changed mtk, exp dw.
% % % % blen=sum(brt);
% % % % dlen=nzf1-blen;
% % % % ncol=length(pick);
% % % % 
% % % % % ncol=0;
% % % % % nbf2=length(sigf2(:,1));
% % % % % for s=1:nbf2
% % % % %     ncol=ncol+(sigf2(s,2)-sigf2(s,1))+1;
% % % % % end
% % % % 
% % % % if ncol==1
% % % %     % random point coordinate in 1D matrix
% % % %     brand=zeros(npts,1);
% % % %     drand=zeros(npts,1);
% % % % 
% % % %     for p=1:npts
% % % %         brand(p,:)=round(rand(1)*blen);
% % % %         if brand(p,1)==0
% % % %             brand(p,1)=1;
% % % %         end
% % % %     end
% % % %     for p=1:npts
% % % %         drand(p,:)=round(rand(1)*dlen);
% % % %         if drand(p,1)==0
% % % %             drand(p,1)=1;
% % % %         end
% % % %     end
% % % % else
% % % %     % random point coordinate in 2D matrix
% % % %     brand=zeros(npts,2);
% % % %     drand=zeros(npts,2);
% % % %     
% % % %     for p=1:npts
% % % %         brand(p,:)=round(rand(1,2).*[blen,ncol]);
% % % %         if brand(p,1)==0
% % % %             brand(p,1)=1;
% % % %         end
% % % %         if brand(p,2)==0
% % % %             brand(p,2)=1;
% % % %         end
% % % %     end
% % % %     for p=1:npts
% % % %         drand(p,:)=round(rand(1,2).*[dlen,ncol]);
% % % %         if drand(p,1)==0
% % % %             drand(p,1)=1;
% % % %         end
% % % %         if drand(p,2)==0
% % % %             drand(p,2)=1;
% % % %         end
% % % %     end
% % % % end
% % % % % <<<<<<<<<<<<<<


%     %~*~*~*~*~*~* make mask ~*~*~*~*~*~*~
%     spec=zeros(nzf2,nzf1); % fake spectrum
%     list=plist;
%     step=1;
%     refn=1;
%     nt1=nzf1;
%     nt2=nzf2;
%     [sigf1 sigf2]=brightPeak_v3(spec,list,step,refn,...
%         sw1,sw2,nt1,nt2,obs1,obs2,cent1,cent2);
    
%     %~*~*~*~*~*~* make 'pick' ~*~*~*~*~*~*~
%     pick=[];
%     for n=1:length(sigf2(:,1))
%         yes=0;
%         temp=ceil(mean([sigf2(n,1),sigf2(n,2)]));
%         if isempty(pick)
%             pick=[pick;temp];
%         else
%             if sum(pick==temp)
%             else
%                 pick=[pick;temp];
%             end
%         end
%     end
% ooooooooooooooooooooooooooooooooooooooooooooo
% look up "pick" on the fly. mtk, dec1914.
if isempty(pick)
%     disp('pick is empty')

    temp=zeros(nzf2,1);
    for n=1:size(sigf2,1)
        temp(sigf2(n,1):sigf2(n,2),1)=cycm1(sigf2(n,1):sigf2(n,2),1);
    end
    
    [a,b]=max(temp);
    pick=b;
    clear temp
end
% ooooooooooooooooooooooooooooooooooooooooooooo

ncol=length(pick(:,1));
    
    %~*~*~*~*~*~* make a checker matrix ~*~*~*~*~*~*~
    db01=zeros(nzf2,nzf1);
    for m=1:nzf2
        for n=1:nzf1
            bright=0;
            for o=1:length(sigf1(:,1))
                if m>=sigf2(o,1) && m<=sigf2(o,2) && n>=sigf1(o,1) && n<=sigf1(o,2)
                    bright=1;
                end
            end
            if bright==1
                db01(m,n)=1;
            end
        end
    end
    %~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~* 

    
    
    
    
    for cyc=1:ncyc
%        if sdymrg<dymrg
%         sdymrg=dymrg;
%         disp('dynamic range=');
%         disp(sdymrg);

%=================================================
% for checking dynamic range and RMSD
%         tochk1=zeros(ncol,nzf1);
% modified jun1510, ym, to use z-fill in t1.
        tochk1=zeros(ncol,length(cycm1(1,:)));
%         loc=1;
        for r=1:ncol
%             num=sigf2(r,2)-sigf2(r,1)+1;
%             tochk1(loc:num,:)=cycm1(sigf2(r,1):sigf2(r,2),:);
%             loc=loc+num;
            tochk1(r,:)=cycm1(pick(r),:);
        end
        tochk1=tochk1.';
        % #-of-t1-pts by ncol matrix.

%         cycmv(:,1)=cycm1(216,:);
%         cycmv(:,1)=cycm1(refrow,:);

%          tochk1=cycm1(refrow,:);
%          tochk1=sineb(tochk1,90,1.0,2);

    % should be commented out eventually. in real cases no use.
    %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%         if wndon==1
%             cycfid=sineb(cycm1(226,:).',90,wndend1,1);

% jun2708, fri
% fid rmsd should be taken at all signal-containing F1 slices.
            cycfid=sineb(cycm1.',90,wndend1,1);
%         end
%             rmsd=fidrmsd(orgfid,cycfid,[1]);

    %\\\\\\\\\\\\
    % commented out, feb1909, thu, matsuki.
    % in a de novo case, there is no full reference FID.
    %\\\\\\\\\\\\
%     rmsd=fidrmsd(orgfid,cycfid,pick);
    % <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    
    %         disp('rmsd =');
%         disp(abs(rmsd));

%          tochk1=fft(tochk1,nzf1,2);
%          tochk1=fft(tochk1(refrow,:),nzf1,2);
        if wndon==1
% jun2408,tue,mtk
% to maximize s/n by sift, check exp instead of sine.
            tochk1=sineb(tochk1,90,wndend1,1);
%             tochk1=expo(tochk1,2,sw1,1);
        end
         tochk1=fft(tochk1,nzf1,1);
        % nzf1 by ncol matrix.

        
        %          tochk1=mphasec(tochk1,4000,[0,0],1);
         tochk1=mphasec(tochk1,sw1,ph1,1);
         
        
        p=1;
        for m=1:length(pick(:,1))
            for n=1:nzf1
                if db01(pick(m),n)==0 % jul2513, mtk.
%                 if db01(m,n)==0
                    darknoise(p,1)=tochk1(n,m);
                    p=p+1;
                end
            end
        end
        
       mx=max(max(real(tochk1)));
       ndrk=median(abs(darknoise(:,1)));
       drdrk=mx/ndrk;
        prfl(cyc,1)=drdrk;     % dyn range at dark area
        prfl(cyc,2)=mx;     % max signal intensity
        prfl(cyc,3)=ndrk;     % noise level at dark area
         
         
         
% % % % %         [dymrg,o1,o2,o3,o4]=drange(tochk1,brt,noise);
% % % %         [drbrt,drdrk,mx,nbrt,ndrk]=dranmedian2(tochk1,brt,brand,drand);
% % % % %         [dymrg,o1,o2,o3]=drange(tochk1,sigf1,sigf2,refrow,noisst,noised);
% % % % %         [dymrg,o1,o2,o3]=drange(cycm1,sigf1,sigf2,refrow,noisst,noised);
% % % % %         mx(cyc,:)=o1;
% % % % %         mstd(cyc,:)=o2;
% % % % %         mdzsq(cyc,:)=o3;
% % % %         prfl(cyc,1)=drbrt;  % dyn range at bright area
% % % %         prfl(cyc,2)=drdrk;     % dyn range at dark area
% % % %         prfl(cyc,3)=mx;     % max signal intensity
% % % %         prfl(cyc,4)=nbrt;     % noise level at bright area
% % % %         prfl(cyc,5)=ndrk;     % noise level at dark area

% jun2708, fri
% fid rmsd taken at all signal-containing slices. so average it.
%         prfl(cyc,6)=rmsd;   % FID RMSD

    %\\\\\\\\\\\\
    % commented out, feb1909, thu, matsuki.
    % in a de novo case, there is no full reference FID.
    %\\\\\\\\\\\\
%     prfl(cyc,6)=mean(rmsd);   % FID RMSD
%=================================================

%         tic
        cycm1=fft(cycm1,nzf1,2);
        cycm2=fft(cycm2,nzf1,2);
%         toc
%         disp('FT1, done.');

% off, apr0108
% select block zeroing algorithm or thresholding.
    if blkzero==1
%         cycm1=fzr(cycm1,sigf1,sigf2);
%         cycm2=fzr(cycm2,sigf1,sigf2);
        cycm1=fzrPeak(cycm1,sigf1,sigf2);
        cycm2=fzrPeak(cycm2,sigf1,sigf2);
    else
        cycm1=threhard(cycm1,sw1,ph1,thre);
        cycm2=threhard(cycm2,sw1,ph1,thre);
    end
   

%         tic
        cycm1=ifft(cycm1,[],2);
        cycm2=ifft(cycm2,[],2);
%         toc
%         disp('IFT1, done.');

        cycm1=cycm1(:,1:nt1org);
        cycm2=cycm2(:,1:nt1org);

%         cycm1=rplc(temp1,cycm1,timetab,tailc);
%         cycm2=rplc(temp2,cycm2,timetab,tailc);
%         cycm1=rplc(temp1,cycm1,timetab,tailc,sigf2);
%         cycm2=rplc(temp2,cycm2,timetab,tailc,sigf2);
        cycm1=rplc(temp1,cycm1,timetab,tailc,[1 nzf2]);
        cycm2=rplc(temp2,cycm2,timetab,tailc,[1 nzf2]);
%         disp('replace, done.');

%         disp('cycle done is')
        disp(cyc);
%         else
%             disp('break has occurd at cyc=')
%             disp(cyc);    
%             disp('dynamic range=');
%             disp(dymrg);
%             disp('rmsd =');
%             disp(rmsd);
%             break
%        end
    end
    
    [cycm1,cycm2]=rrcplx(cycm1,cycm2);
%     disp('swap for FT2, done.');
    
%     cycm1=mphasec(cycm1,8012.8,[-88 0],1);
%     cycm2=mphasec(cycm2,8012.8,[-88 0],1);
    cycm1=mphasec(cycm1,sw2,-1*ph2,1);
    cycm2=mphasec(cycm2,sw2,-1*ph2,1);

    omat1=ifft(cycm1,[],1);
    omat2=ifft(cycm2,[],1);
%     disp('IFT2, done.');

omat1=omat1(1:nt2org,:);
omat2=omat2(1:nt2org,:);

toc
disp('All, done.');


