unit spectrum;

interface

uses fourier,GlType,GlProc,MathObj,WinFuncs,Classes,SysUtils, MVTU_TLB;

type
 //Cвойства спектра
 TSpectrProp=record
  Size:^integer;          //Размер серии
  CalcMode:^integer;      //Способ расчёта- 0-по всем сериям ,1-по каждой серии
  Tau:^double;            //Период квантования
  DelTrend:^integer;      //Удаление линейного тренда- 0-нет , 1- да
  Win:^integer;           //Тип окна данных
  OutMode:^integer;       //Способ вывода-  0-абсолютная плотность ,1-нормированная
 end;
 PSpectrProp=^TSpectrProp;

 TSpectrNumericStates=record
  Position:integer;
  Num,
  SumSqr,
  Sum,
  SW,
  WSum,
  f0,
  NY,
  N2,
  MT:double;
  end;

 TSpectrState=record
  Data:     PComplexArr;
  W,Pwr,REZ:TExtArray;
  NumStates:TSpectrNumericStates;
  end;

 PSpectrState=^TSpectrState;

  //RUN-функция вычисления спектральной плотности.
function TStatSpectr(at,adt:RealType;var time,dtime:RealType;var AU,AY:TPtrExt;var AX,ADX:array of RealType;
  var CU,CY : TIntArray;var Prop:TPtrArray;var Vars:PSpectrState;Action:Integer):Integer;export;

function TStatSpectrRst(at,adt:RealType;var AU,AY:TPtrExt;var AX,ADX:PExtArr;
 var Prop:TPtrArray;var Vars : Pointer;Action:Integer;var F:File):Integer;export;

implementation

function TStatSpectrConvert(Action:  integer;
                          Vars:    Pointer;
                          Prop:    TPtrArray;
                          PropStr: TStringList;
                          BlockId: integer;
                          var MVTU: IMVTU_Server):integer;
  const RecName:PChar = 'Спектральная плотность';
  var
      Res: integer;
begin
  case Action of
    //Возвращаем ссылку на имя записи в базе МВТУ-4
    cnv_GetRecName: Result:=integer(RecName);
    //Конвертируем и устанавливаем параметры при помощи интерфейса МВТУ
    cnv_Convert: with PSpectrProp(Prop.Arr)^ do try
                   //Присваиваем параметры для блока
                   MVTU.SetBlockProp(BlockId,'size',PChar(PropStr[0]),Res);
                   MVTU.SetBlockProp(BlockId,'calcmode',PChar(IntToStr(CalcMode^)),Res);
                   MVTU.SetBlockProp(BlockId,'tau',PChar(PropStr[2]),Res);
                   MVTU.SetBlockProp(BlockId,'deltrend',PChar(IntToStr(DelTrend^)),Res);
                   MVTU.SetBlockProp(BlockId,'win',PChar(IntToStr(Win^)),Res);
                   MVTU.SetBlockProp(BlockId,'outmode',PChar(IntToStr(OutMode^)),Res);
                 finally
                 end;
  end;
end;

function TStatSpectr;
 label L0;
 var M,N,D,U,SumW,K,CNY,CN2,b,CM,a:double;
  i:integer;
begin
 Result:=0;
 with Vars^,Vars^.NumStates,PSpectrProp(Prop.Arr)^ do
 case Action of
  f_GetConvertFuncAdr: Result:=integer(@TStatSpectrConvert);
  f_GetDeriCount,
  f_GetInit      :Result:=0;
  f_GetStateCount:Result:=0;
  f_GetCount     :begin
   CY.arr^[0]:=Size^ div 2 -1;
   CY.arr^[1]:=Size^ div 2 -1;
   end;
  f_InitMem      :begin
   W.ChangeCount(Size^);
   PWR.ChangeCount(CY.arr^[1]);
   REZ.ChangeCount(CY.arr^[1]);
   ReallocMem(Data,SOfCplx*Size^);
   end;
  f_InitTime     :time:=at;
  f_InitState    :begin
   Position:=0;
   Num:=0;
   Sum:=0;
   SumSqr:=0;
   SW:=0;
   WSum:=0;
   NY:=0;
   N2:=0;
   MT:=0;
   for i:=0 to Size^-1 do begin
    U:=winfunc(i,Size^,Win^);
    W.arr^[i]:=U;
    WSum:=WSum+U*U;
    end;
   PWR.FillArray(0);
   AY.Ptr(0).FillArray(0);
   AY.Ptr(1).FillArray(0);
   goto L0;
   end;
  f_Create       :begin
   Vars:=New(PSpectrState);
   with PSpectrState(Vars)^ do begin
    W:=TExtArray.Create(1);
    PWR:=TExtArray.Create(1);
    REZ:=TExtArray.Create(1);
    GetMem(Data,128*SOfCplx);
   end
   end;
  f_Free         :begin
   W.Free;
   PWR.Free;
   REZ.Free;
   Dispose(Data);
   Dispose(Vars);
   end;
  f_RestoreOuts: begin
    Move(REZ.arr^,AY.Ptr(0).arr^,REZ.Count*SOfR);
    for i:=0 to REZ.Count-1 do AY.Ptr(0).arr^[i]:=(i+1)*f0;
   end;
  f_GoodStep     : if time-at<=adt/2 then begin
   L0:Data^[Position].Re:=AU.Ptr(0).arr^[0];
      Data^[Position].Im:=0;
   if Position=Size^-1 then begin
     M:=Sum;D:=SumSqr;N:=Num+Size^;CNY:=NY;CN2:=N2;CM:=MT;
     for i:=0 to Size^-1 do begin
      U:=Data[i].Re;
      M:=M+U;
      K:=Num+i+1;
      CNY:=CNY+K*U;
      CN2:=CN2+sqr(K);
     end;
     SumW:=SW+WSum;
     if DelTrend^=1 then begin
      b:=(CNY-0.5*M*N)/(CN2-N*sqr(0.5*N));
      a:=M/N-b*0.5*N
     end else begin
      b:=0;
      a:=M/N
     end;
     for i:=0 to Size^-1 do begin
      U:=(Data[i].Re-(b*(Num+i+1)+a));
      D:=D+U*U;
      CM:=CM+U;
      Data[i].Re:=U*W.arr^[i]
     end;
     FFT(Data,Size^,sgn_Fw);
     for i:=0 to AY.Ptr(0).Count-1 do
      AY.Ptr(1).arr^[i]:=Pwr.arr^[i]+sqr(Data^[i+1].Re)+sqr(Data^[i+1].Im);
     if (CalcMode^=0) then begin
      Sum:=M;SumSqr:=D;Num:=N;SW:=SumW;NY:=CNY;N2:=CN2;MT:=CM;
      Move(AY.Ptr(1).arr^,Pwr.arr^,AY.Ptr(1).Count*SOfR);
     end;
     if Tau^=0 then begin
       K:=2*adt/SumW;
       f0:=1/(Size^*adt)
      end else begin
       K:=2*Tau^/SumW;
       f0:=1/(Size^*Tau^)
      end;
     if OutMode^=1 then K:=K*(N-1)/(D-sqr(CM)/N);
     for i:=0 to AY.Ptr(0).Count-1 do begin
      AY.Ptr(0).arr^[i]:=(i+1)*f0;
      AY.Ptr(1).arr^[i]:=AY.Ptr(1).arr^[i]*K;
     end;
     Move(AY.Ptr(1).arr^,REZ.arr^,REZ.Count*SOfR)
    end;
   if (Position=Size^-1) then Position:=0 else inc(Position);
   time:=time+Tau^;
   end;
  end; //END Case
end;   //END Procedure

function TStatSpectrRst;
begin
 Result:=0;
 with PSpectrProp(Prop.arr)^,PSpectrState(Vars)^ do
 case Action of
 f_ReadRez   : try
                BlockRead(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_ReadFile
               end;
 f_WriteRez  : try
                BlockWrite(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_WriteFile
               end;
 f_ReadRst   : try
                BlockRead(F,NumStates,SizeOf(NumStates));
                BlockRead(F,Data^,Size^*SOfCplx);
                BlockRead(F,W.arr^,W.Count*SOfR);
                BlockRead(F,PWR.Arr^,Pwr.Count*SOfR);
                BlockRead(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_ReadFile
               end;
 f_WriteRst  : try
                BlockWrite(F,NumStates,SizeOf(NumStates));
                BlockWrite(F,Data^,Size^*SOfCplx);
                BlockWrite(F,W.arr^,W.Count*SOfR);
                BlockWrite(F,PWR.Arr^,Pwr.Count*SOfR);
                BlockWrite(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_WriteFile
               end;
 end
end;

end.