unit DblCorrel;

interface
uses glType,glProc,Mathobj,Fourier,Classes,SysUtils, MVTU_TLB;

type
  //Cвойства функции корреляции
 TCorelProp=record
  Size:^integer;          //Размер серии
  CalcMode:^integer;      //Способ расчёта- 0-по всем сериям ,1-по каждой серии
  Tau:^double;            //Период квантования
  DelTrend:^integer;      //Удаление линейного тренда 0- не удалять, 1-удалять
  end;

 PCorelProp=^TCorelProp;

 TDblCorelNumericStates=record
  Size2,
  Position:integer;
  SerNum,
  T,
  RXX,
  RYY,
  Sum1,Sum2,
  NY1,NY2,
  NN1:double;
  end;

 TCorelState=record
  Data,X,Y,GXY:PComplexArr;
  REZ:TExtArray;
  NumStates:TDblCorelNumericStates;
  end;
 PCorelState=^TCorelState;

      //RUN-функция вычисления взаимной корреляции
      //методом быстрого преобразования Фурье
function TStatCorel(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:PCorelState;Action:Integer):Integer;export;

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

const
 NULL:TComplex=(Re:0;Im:0);

implementation

function TStatCorelConvert(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 PCorelProp(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);
                 finally
                 end;
  end;
end;

function TStatCorel;
 label L0;
 var K,XX,SN,Num,YY,M1,M2,CNY1,CNY2,CNN1,U,a1,b1,a2,b2,N:double;
  i:integer;
begin
 Result:=0;
 with Vars^,PCorelProp(Prop.Arr)^,Vars^.NumStates do
 case Action of
  f_GetConvertFuncAdr: Result:=integer(@TStatCorelConvert);
  f_GetCount     :begin
   CU.arr^[1]:=CU.arr^[0];
   CY.arr^[0]:=Size^;
   CY.arr^[1]:=Size^;
   end;
  f_InitMem      :begin
   Size2:=2*Size^;
   ReallocMem(Data,Size2*SOfCplx);
   ReallocMem(X,Size2*SOfCplx);
   ReallocMem(Y,Size2*SOfCplx);
   ReallocMem(GXY,Size2*SOfCplx);
   REZ.ChangeCount(CY.arr^[1]);
   end;
  f_InitTime     :time:=at;
  f_InitState    :begin
   Position:=0;
   SerNum:=0;
   RXX:=0;
   RYY:=0;
   Sum1:=0;Sum2:=0;
   NY1:=0;NY2:=0;
   NN1:=0;
   for i:=0 to Size2-1 do GXY^[i]:=NULL;
   AY.Ptr(0).FillArray(0);
   AY.Ptr(1).FillArray(0);
   goto L0;
   end;
  f_Create       :begin
   Vars:=New(PCorelState);
   with PCorelState(Vars)^ do begin
    REZ:=TExtArray.Create(1);
    GetMem(Data,128*SOfCplx);
    GetMem(X,128*SOfCplx);
    GetMem(Y,128*SOfCplx);
    GetMem(GXY,128*SOfCplx)
    end;
   end;
  f_Free         :begin
   Dispose(Data);
   Dispose(X);
   Dispose(Y);
   Dispose(GXY);
   REZ.Free;
   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*T;
   end;
  f_GoodStep     : if time-at<=adt/2 then begin
   L0:Data^[Position].Re:=AU.Ptr(0).arr^[0];
      Data^[Position].Im:=AU.Ptr(1).arr^[0];
   if Position=Size^-1 then begin
    SN:=SerNum+1;
    Num:=SerNum*Size^;
    N:=SN*Size^;
    XX:=RXX;YY:=RYY;
    M1:=Sum1;M2:=Sum2;
    CNY1:=NY1;CNY2:=NY2;
    CNN1:=NN1;
    for i:=0 to Size^-1 do begin
     K:=Num+i+1;
     U:=Data^[i].Re;
     M1:=M1+U;
     CNY1:=CNY1+K*U;
     CNN1:=CNN1+sqr(K);
     U:=Data^[i].Im;
     M2:=M2+U;
     CNY2:=CNY2+K*U;
     end;
    if DelTrend^=1 then begin
     b1:=(CNY1-0.5*M1*N)/(CNN1-N*sqr(0.5*N));
     a1:=M1/N-b1*0.5*N;
     b2:=(CNY2-0.5*M2*N)/(CNN1-N*sqr(0.5*N));
     a2:=M2/N-b2*0.5*N;
    end else begin
     b1:=0;
     b2:=0;
     a1:=M1/N;
     a2:=M2/N
    end;
    for i:=0 to Size^-1 do begin
     U:=(Data^[i].Re-(b1*(Num+i+1)+a1));
     XX:=XX+U*U;
     Data^[i].Re:=U;
     U:=(Data^[i].Im-(b2*(Num+i+1)+a2));
     YY:=YY+U*U;
     Data^[i].Im:=U;
     end;
    for i:=Size^ to Size2-1 do Data^[i]:=NULL;
    DoubleFFT(Data,X,Y,Size2);
    if Tau^= 0 then T:=adt else T:=Tau^;
     for i:=0 to Size2-1 do begin
      K:=GXY^[i].Re+X^[i].Re*Y^[i].Re+X^[i].Im*Y^[i].Im;
      X^[i].Im:=GXY^[i].Im+X^[i].Re*Y^[i].Im-Y^[i].Re*X^[i].Im;
      X^[i].Re:=K
      end;
    if CalcMode^=0 then begin
      SerNum:=SN;
      RXX:=XX;
      RYY:=YY;
      Sum1:=M1;Sum2:=M2;
      NY1:=CNY1;NY2:=CNY2;
      NN1:=CNN1;
      Move(X^,GXY^,Size2*SOfCplx);
     end;
    FFT(X,Size2,sgn_Rew);
    for i:=0 to Size^-1 do AY.Ptr(0).arr^[i]:=T*i;
    if (XX>0) and (YY>0) then begin
      K:=0.5/(sqrt(XX*YY)*(Size^-1));
      for i:=0 to Size^-1 do AY.Ptr(1).arr^[i]:=K*X^[i].Re
    end else AY.Ptr(1).FillArray(1);
    end;
   Move(AY.Ptr(1).arr^,REZ.arr^,REZ.Count*SOfR);
   if (Position=Size^-1) then Position:=0 else inc(Position);
   time:=time+Tau^;
   end;
  end; //END Case
end;   //END Procedure

function TStatCorelRst;
begin
 Result:=0;
 with PCorelProp(Prop.arr)^,Vars^,Vars^.NumStates 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^,Size2*SOfCplx);
                BlockRead(F,GXY^,Size2*SOfCplx);
                BlockRead(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_ReadFile
               end;
 f_WriteRst  : try
                BlockWrite(F,NumStates,SizeOf(NumStates));
                BlockWrite(F,Data^,Size2*SOfCplx);
                BlockWrite(F,GXY^,Size2*SOfCplx);
                BlockWrite(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_WriteFile
               end;
 end
end;


end.