function [fetal_peaks_round1,Ftest_round1,sample2,find1,fetal_ind_cleaned1,qflag]=Fetal_Peakdetection2(FECG,sample,fs,factorcheat,test_flag)
qflag=1;
blocksize=2000;
[nleads,nsamples]=size(FECG);
nblocks=nsamples/blocksize;
bin=400;


factor=factorcheat;


[Ftest_round1,good_sig1]=AddSig(FECG,factor,sample,fs);
fpeaks1 = PeakDetection3(Ftest_round1,fs,sample,.1,4);
find1=find(fpeaks1==1);
if length(find1)<30
    clear find1 fpeaks1 mn1 r1
    sig_power1=sum(good_sig1);
    [mpow,mpowi]=max(sig_power1);
    test_sample=Ftest_round1((mpowi-1)*blocksize+1:mpowi*blocksize);
    [mt,mti]=max(abs(test_sample(26:blocksize-25)));
    mti=mti+25;
    fpeaks1 = PeakDetection3(Ftest_round1,fs,test_sample(mti-25:mti+24),.1,4);
    find1=find(fpeaks1==1);
end
% % % % assert(length(find1)>30,'sample does not match with Ftest_round1 quite well. Very few peaks found.');
if length(find1) < 50
    qflag=0;
    fetal_peaks_round1=[];
    Ftest_round1=[];
    sample2=0;
    find1=[];
    fetal_ind_cleaned1=[];
    return;
end


%checking for untrustworthy peaks
dfetal=diff(find1);
fetal_ind_test1=dfetal>=600;
fetal_ind_test2= fetal_ind_test1+[0 fetal_ind_test1(1:end-1)]+[fetal_ind_test1(2:end) 0];
fetal_ind_test3=filter([0.5 0 0.5],1,dfetal);
fetal_ind_test3=[fetal_ind_test3(3) fetal_ind_test3(3:end) fetal_ind_test3(end)];
ind_set=find(abs(fetal_ind_test3-dfetal).*~fetal_ind_test2>100)+1;

dd=find(dfetal<max(200,mean(dfetal(~fetal_ind_test1))/1.5));
ind_set=union(ind_set,union(dd,dd+1));


fetal_ind_cleaned1=setdiff(find1,find1(ind_set));

if ~test_flag

    fetal_peaks_round1=FillMissedPeaks3(fetal_ind_cleaned1,fs);

    fetal_peaks_round1=sort(fetal_peaks_round1);

    clear dfetal fetal_ind_test1 fetal_ind_test2 fetal_ind_test3 ind_set dd


    finalpeaks1=zeros(1,60000);
    finalpeaks1(fetal_peaks_round1)=1;
    phase = PhaseCalculation(finalpeaks1);


    nanind=find(sum(good_sig1)==0);
    if length(nanind)>=1
        FF=Ftest_round1(~isnan(Ftest_round1));
        pp=phase(~isnan(Ftest_round1));
    %     pp=pphase(1:(nanind(1)-1)*blocksize);
    else
        FF=Ftest_round1;
        pp=phase;
    end

    Fmean = MeanECGExtraction(FF,pp,bin,1);

    [Fsigmax FI] = max(abs(Fmean));
    sample2=Fmean(FI-30:FI+29);
%     sample2=Fmean(FI-30:FI+99);
else
    fetal_peaks_round1=0;
    FF=Ftest_round1(~isnan(Ftest_round1));
    ll=PeakDetection3(FF,fs,sample,0.1,4);
    pp=PhaseCalculation(ll);
    Fmean=MeanECGExtraction(FF,pp,bin,1);
    [dummy, FI] = max(abs(Fmean));
    sample2=Fmean(FI-30:FI+29);
%     sample2=Fmean(FI-30:FI+99);
end


end




function [Ftest, good_sig]=AddSig(FECG,factor,sample,fs)

blocksize=2000;
leadset=factor~=0;
[nleads,nsamples]=size(FECG);
nblocks=nsamples/blocksize;
pinterval=25:41;
interval=21:38;
F=0:150;
Ftest=zeros(1,nsamples);
maxFtest=zeros(1,nblocks);
psd_peak=zeros(1,nblocks);
good_sig=zeros(4,nblocks);


for j=1:nblocks
    for i=1:4    
        [S(:,i),F,T,P(:,i)] = spectrogram(FECG(i,(j-1)*blocksize+1:j*blocksize),blocksize,[],F,fs);
    end

    peakf=PsdPeak2(P,pinterval,leadset);
    psd_peak(j)=peakf;
    Fqindex=logical(zeros(1,4));
    Tqindex=logical(zeros(1,4));
    qindex=logical(zeros(1,4));
    Pqindex=logical(zeros(1,4));

    for i=1:4
%         [dummy FI] = max(abs(FECG(i,(j-1)*blocksize+1:j*blocksize)));
        if leadset(i)
            peaktest=PeakDetection3(FECG(i,(j-1)*blocksize+1:j*blocksize),fs,factor(i)*sample,.01,4);
            peakdiff=diff(find(peaktest==1));
            numpeak=length(peakdiff);
            Tqindex(i)=sum(peakdiff>250 & peakdiff<600) >= 2/3 * numpeak & numpeak>=2 & numpeak <=5 ;
            Pqindex(i)= all(peakdiff>250 & peakdiff<600) & sum(peakdiff)>=0.65*blocksize & all(abs(peakdiff-sum(peakdiff)/numpeak)<100);%& (numpeak==3 | numpeak==4);

        end
    end
    if peakf>0
        Larea=sum(P(1:interval(1)-1,:));
        Rarea=sum(P(interval(end)+1:end,:));
        area_test= Larea./Rarea>=1.2 | sum(P(interval,:))./sum(P)>0.4;% & Rarea./sum(P)<0.275);
        peak_testL=sum(P(peakf-5:peakf+5,:))./Larea>=0.21;
        peak_testR=sum(P(peakf-5:peakf+5,:))./Rarea>=0.21;
        Fqindex=area_test & peak_testL & peak_testR;
        set=setdiff(peakf-2:peakf+2,peakf);
    else
        Larea(Pqindex)=sum(P(1:interval(1)-1,Pqindex));
        Rarea(Pqindex)=sum(P(interval(end)+1:end,Pqindex));
        area_test(Pqindex)= Larea(Pqindex)./Rarea(Pqindex)>=1.2 | sum(P(interval,Pqindex))./sum(P(:,Pqindex))>0.4;

        peak_testL(Pqindex)=sum(P(interval,Pqindex))./Larea(Pqindex)>=0.21;
        peak_testR(Pqindex)=sum(P(interval,Pqindex))./Rarea(Pqindex)>=0.21;
        Fqindex(Pqindex)=area_test(Pqindex) & peak_testL(Pqindex) & peak_testR(Pqindex);
    end
    if peakf>0
        qindex=Tqindex & Fqindex | Pqindex;
    else
        qindex(Pqindex)=Fqindex(Pqindex);
    end
  
    if any(qindex)
        Ftest((j-1)*blocksize+1:j*blocksize)= factor(qindex)*FECG(qindex,(j-1)*blocksize+1:j*blocksize)/sum(qindex);
        good_sig(:,j)=qindex;
        maxFtest(j)=max(Ftest((j-1)*blocksize+1:j*blocksize));
    else
        Ftest((j-1)*blocksize+1:j*blocksize)=NaN;
    end
    
    clear S T P Larea Rarea;
end
end