--- spectrum.pas	Mon Jun 26 12:53:55 2000
+++ Spectrum.pas	Mon Jun 26 02:04:34 2000
@@ -1227,8 +1227,8 @@
   DX,X:TXValue;
   Gsum,Y:TYValue;
 begin
-  if count<16 then exit;
-  DX:=points^[1].X-points^[0].X;
+  if (count<16)or(S=0) then exit;

+  DX:=(points^[count-1].X-points^[0].X)/(count-1);

   ig:=0; Gsum:=0; X:=0;
   repeat
     Y:=1/sqrt(2*pi)/S*exp(-(X*X)/(S*S));
@@ -1286,8 +1286,7 @@
     BaseMaxDiff .. maximal absolute value of second order derivative
                    of course, if zero, automatic tuning is applied
 }
-
-procedure XYClassBase(points:PExpPoints;count:integer;isbase:PBoolArray;
+procedure XYClassBase2Der(points:PExpPoints;count:integer;isbase:PBoolArray;

                       BaseMinInterval,BaseMaxDiff:single);
 var
   diffp: PExpPoints;
@@ -1300,8 +1299,8 @@
   {Compute filtered second order derivative}
   GetMem(diffp,count*sizeof(TExpPoint));
   Move(points^,diffp^,count*sizeof(TExpPoint));
-  XYFilterGaus(diffp,count,BaseMinInterval/2);
   XYDiff2Ord(diffp,count);
+  XYFilterGaus(diffp,count,BaseMinInterval/2);

   {Find appropriate value when BaseMaxDiff is zero}
   if BaseMaxDiff=0 then
   begin
@@ -1356,6 +1355,117 @@
   FreeMem(diffp,count*sizeof(TExpPoint));
 end;
 
+{ Classify points as base or candidates for peaks

+  and store result in isbase

+    BaseMinInterval .. minimal length of couse satisfiing conditions

+    BaseMaxSn .. maximal standart deviation of baseline course,

+                 if zero, automatic tuning is applied

+}

+procedure XYClassBaseStdDev(points:PExpPoints;count:integer;isbase:PBoolArray;

+                            BaseMinInterval,BaseMaxStdDev:single);

+var

+  disp: PExpPoints;

+  i, j, ib, iblock: integer;

+  LastX, DX: TXValue;

+  npoints: integer;

+  SumY,SumY2:double;

+  SumYNxt,SumY2Nxt:double;

+  ADisp:single;

+  BaseMaxDisp:single;

+begin

+

+  {Consider constant sampling frequency -- It is not ptimal}

+  DX:=(points^[count-1].X-points^[0].X)/(count-1);

+  if DX=0 then exit;

+  npoints:=round(BaseMinInterval/DX);

+  if npoints<10 then npoints:=10;

+  if count<2*npoints then exit;

+

+  {Compute filtered standard deviation}

+  GetMem(disp,count*sizeof(TExpPoint));

+  Move(points^,disp^,count*sizeof(TExpPoint));

+  XYFilterGaus(disp,count,BaseMinInterval/4);

+

+  {Find appropriate value when BaseMaxStdDev is zero}

+  if BaseMaxStdDev=0 then

+  begin

+    BaseMaxStdDev:=1e-5;

+  end;

+  BaseMaxDisp:=Sqr(BaseMaxStdDev)*Sqr(npoints);

+

+  {Compute disperse for moving window

+   Computation must be restarted after short block,

+   because float format can accumulate precission errors

+   A+B+C-A-B-C <> 0 }

+  SumYNxt:=0;SumY2Nxt:=0;

+  i:=0;

+  while i<npoints do

+  begin

+    SumYNxt:=SumYNxt+disp^[i].Y; SumY2Nxt:=SumY2Nxt+Sqr(disp^[i].Y);

+    inc(i)

+  end;

+  j:=0;

+  while i<count do

+  begin

+    iblock:=i+npoints;

+    if iblock>count then iblock:=count;

+    SumY:=SumYNxt;SumY2:=SumY2Nxt;

+    SumYNxt:=0;SumY2Nxt:=0;

+    repeat

+      ADisp:=SumY2*npoints-Sqr(SumY);

+      SumY:=SumY-disp^[j].Y; SumY2:=SumY2-Sqr(disp^[j].Y);

+      SumY:=SumY+disp^[i].Y; SumY2:=SumY2+Sqr(disp^[i].Y);

+      SumYNxt:=SumYNxt+disp^[i].Y; SumY2Nxt:=SumY2Nxt+Sqr(disp^[i].Y);

+      disp^[j].Y:=ADisp;

+      inc(i);inc(j);

+    until i>=iblock;

+  end;

+  npoints:=npoints div 2;

+  {Center compuded disperse around corresponding X}

+  dec(i);dec(j);

+  while i>count-npoints do

+  begin

+    disp^[i].Y:=disp^[j].Y;

+    dec(i);

+  end;

+  while j>=0 do

+  begin

+    disp^[i].Y:=disp^[j].Y;

+    dec(i);dec(j);

+  end;

+  while i>=0 do

+  begin

+    disp^[i].Y:=disp^[i+1].Y;

+    dec(i);

+  end;

+

+  {Classify points}

+  LastX:=disp^[0].X; ib:=0;

+  for i:=0 to count-1 do

+  begin

+    if disp^[i].Y<=BaseMaxDisp then

+    begin

+      if (disp^[i].X-LastX)>=BaseMinInterval then

+      begin

+        while ib<=i do

+        begin

+          isbase^[ib]:=true;

+          inc(ib);

+        end;

+      end;

+    end else begin

+      LastX:=disp^[i].X;

+      while ib<=i do

+      begin

+        isbase^[ib]:=false;

+        inc(ib);

+      end;

+    end;

+  end;

+  FreeMem(disp,count*sizeof(TExpPoint));

+end;

+

+{ Find index of nearest [X,Y] point for given X }

 function XYFindXIndex(points:PExpPoints;count:integer;
                       Xval:TXValue):integer;
 var
@@ -1375,6 +1485,7 @@
   Result:=ib;
 end;
 
+{ Autodetect baseline }

 function TSpectrum.AutoDetectForBase: boolean;
 var
   X, X1: TXValue;
@@ -1398,7 +1509,8 @@
   GetMem(isbase, cpt * sizeof(boolean));
   bmi := AcqData.ULM.BaseMinInterval;
   bmd := AcqData.ULM.BaseMaxDiff;
-  XYClassBase(pa, cpt, isbase, bmi, bmd);
+  XYClassBaseStdDev(pa, cpt, isbase, bmi, bmd);

+  {XYClassBase2Der(pa, cpt, isbase, bmi, bmd);}

 
   {Create base description}
   ipt:=0;
@@ -1455,8 +1567,9 @@
   baseln: TULBRObj; { in ULBRType.PAS }
   points: TMemoryStream;
   pa, tmppa: PExpPoints;{ulantype}
-  bmi,slope: single;
-  PeakMinWidth:single;
+  slope: single;

+  MinPeakWidth:single;

+  MinPeakHeight:single;

   RejectLevel:single;
 
  procedure StoreNewPeak(PTX:TXValue);
@@ -1468,9 +1581,12 @@
    PHHHit:    boolean;
    PHHLimX:   TXValue;
  begin
+   pti:=XYFindXIndex(tmppa, dpt, PTX);

+   {Reject small peaks under MinPeakHeight}

+   if abs(tmppa^[pti].Y)<MinPeakHeight then exit;

+   {Create new peak object}

    newpeak := TULPRObj(Peaks.Add(ULPRID));
    newpeak.X:= PTX;
-   pti:=XYFindXIndex(tmppa, dpt, PTX);
    PHHY:=tmppa^[pti].Y/2;
    { Find Left Peak End }
    PHHLimX:=0;
@@ -1490,7 +1606,7 @@
      if tmppa^[i].Y<=MinY then begin
        MinX:=tmppa^[i].X; MinY:=tmppa^[i].Y;
      end
-     else if tmppa^[i].X<MinX-2*PeakMinWidth then break;
+     else if tmppa^[i].X<MinX-2*MinPeakWidth then break;

    end;
    newpeak.X1:=MinX;
    { Find Right Peak End }
@@ -1510,12 +1626,12 @@
      if tmppa^[i].Y<=MinY then begin
        MinX:=tmppa^[i].X; MinY:=tmppa^[i].Y;
      end
-     else if tmppa^[i].X>MinX+2*PeakMinWidth then break;
+     else if tmppa^[i].X>MinX+2*MinPeakWidth then break;

    end;
    newpeak.X2:=MinX;
  end;
 
- procedure ProcessOneReg; {imput is tmppa,dpt,PeakMinWidth}
+ procedure ProcessOneReg; {imput is tmppa,dpt,MinPeakWidth}

  var
    i : integer;
    PTB,PTE: TXValue;    {Peak Top Begin/End}
@@ -1534,20 +1650,20 @@
      begin
        if AY>PTY then begin PTB:=tmppa^[i].X; PTE:=PTB; PTY:=AY; end
        else if AY=PTY then begin PTE:=tmppa^[i].X; end
-       else if tmppa^[i].X>=PTB+PeakMinWidth then
+       else if tmppa^[i].X>=PTB+MinPeakWidth then

        begin
 	 if PTY>=RejectLevel then
          begin
 	   PTX:=(PTB+PTE)/2; {Retention Time of New Peak}
-	   if (PTX>tmppa^[0].X+PeakMinWidth/3) and
-              (PTX<tmppa^[dpt-1].X-PeakMinWidth/3) then
+	   if (PTX>tmppa^[0].X+MinPeakWidth/3) and

+              (PTX<tmppa^[dpt-1].X-MinPeakWidth/3) then

              StoreNewPeak(PTX);
 	 end;
 	 FindPT:=False;
        end;
      end else begin
        if AY<PTY then begin PTY:=AY; PTE:=tmppa^[i].X; end
-       else if tmppa^[i].X>=PTE+PeakMinWidth then FindPT:=True;
+       else if tmppa^[i].X>=PTE+MinPeakWidth then FindPT:=True;

      end;
      inc(i);
    until i>=dpt;
@@ -1560,10 +1676,10 @@
     points.CopyFrom(AcqData.ULAD.Data, AcqData.ULAD.Data.Size);
     cpt := AcqData.ULAD.Data.Size div sizeof(TExpPoint);
     pa := PExpPoints(points.Memory);
-    bmi := AcqData.ULM.BaseMinInterval;
+    MinPeakWidth:=AcqData.ULM.MinPeakWidth;

+    MinPeakHeight:=AcqData.ULM.MinPeakHeight;

     Peaks.Clear;
     {Detect peaks above found base ranges}
-    PeakMinWidth:=bmi;
     RejectLevel:=0;
     ibl:=0;
     while ibl<BaseLine.ChildCount do
@@ -1583,7 +1699,7 @@
           tmppa^[ipt3].Y:=tmppa^[ipt3].Y-
              (tmppa^[ipt3].X-baseln.X1)*slope-baseln.Y1;
         end;
-        XYFilterGaus(tmppa,dpt,PeakMinWidth/2); {Not needed for filtered data}
+        XYFilterGaus(tmppa,dpt,MinPeakWidth/2); {Not needed for filtered data}

         ProcessOneReg;
         FreeMem(tmppa, dpt * sizeof(TExpPoint));
       end;
