unit DblSpectr;

interface

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

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

 TDblSpectrNumericStates=record
  Position:integer;
  Num,
  SumW,
  WSum,
  SumSqr1,SumSqr2,
  Sum1,Sum2,
  NY1,NY2,
  NN1,
  MT1,MT2,
  f0:double;
  end;

 TDblSpectrState=record
  Data,X,Y,GXY:PComplexArr;
  W,REZ:TExtArray;
  NumStates:TDblSpectrNumericStates;
  end;
 PDblSpectrState=^TDblSpectrState;

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

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

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

implementation

function TStatDblSpctrConvert(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 PDblSpectrProp(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 TStatDblSpctr;
 label L0;
 var K,N,SW,M1,M2,D1,D2,CNY1,CNY2,CNN1,CM1,CM2,a1,a2,b1,b2,U:double;
  i:integer;
begin
 Result:=0;
 with Vars^,Vars^.NumStates,PDblSpectrProp(Prop.Arr)^ do
 case Action of
  f_GetConvertFuncAdr: Result:=integer(@TStatDblSpctrConvert);
  f_GetCount     :begin
   CU.arr^[1]:=CU.arr^[0];
   CY.arr^[0]:=Size^ div 2 -1;
   CY.arr^[1]:=CY.arr^[0]
   end;
  f_InitMem      :begin
   ReallocMem(Data,Size^*SOfCplx);
   ReallocMem(X,Size^*SOfCplx);
   ReallocMem(Y,Size^*SOfCplx);
   ReallocMem(GXY,Size^*SOfCplx);
   W.ChangeCount(Size^);
   REZ.ChangeCount(CY.arr^[1]);
   end;
  f_InitTime     :time:=at;
  f_InitState    :begin
   Position:=0;
   Num:=0;
   WSum:=0;
   SumW:=0;
   Sum1:=0;Sum2:=0;
   SumSqr1:=0;SumSqr2:=0;
   MT1:=0;MT2:=0;
   NY1:=0;NY2:=0;
   NN1:=0;
   for i:=0 to Size^-1 do begin
    GXY^[i]:=NULL;
    K:=winfunc(i,Size^,Win^);
    W.arr^[i]:=K;
    WSum:=WSum+K*K
    end;
   AY.Ptr(0).FillArray(0);
   AY.Ptr(1).FillArray(0);
   goto L0
   end;
  f_Create       :begin
    Vars:=New(PDblSpectrState);
    with PDblSpectrState(Vars)^ do begin
      W:=TExtArray.Create(1);
      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);
   W.Free;
   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+1)*f0;
   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
    N:=Num+Size^;
    SW:=SumW+WSum;
    M1:=Sum1;M2:=Sum2;
    D1:=SumSqr1;D2:=SumSqr2;
    CNY1:=NY1;CNY2:=NY2;
    CNN1:=NN1;
    CM1:=MT1;CM2:=MT2;
      //Предварительная подготовка данных
    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));
     D1:=D1+U*U;
     CM1:=CM1+U;
     Data^[i].Re:=U*W.arr^[i];
     U:=(Data^[i].Im-(b2*(Num+i+1)+a2));
     D2:=D2+U*U;
     CM2:=CM2+U;
     Data^[i].Im:=U*W.arr^[i]
     end;
      //Вычисление взаимного спектра
    DoubleFFT(Data,X,Y,Size^);
     for i:=0 to Size^-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
      Num:=N;
      SumW:=SW;
      Sum1:=M1;Sum2:=M2;
      SumSqr1:=D1;SumSqr2:=D2;
      NY1:=CNY1;NY2:=CNY2;
      NN1:=CNN1;
      MT1:=CM1;MT2:=CM2;
      Move(X^,GXY^,Size^*SOfCplx);
     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^=2 then for i:=0 to Size^ div 2 -2 do begin
      AY.Ptr(0).arr^[i]:=f0*i;
      if X^[i].Re<>0 then AY.Ptr(1).arr^[i]:=arctan(X^[i].Im/X^[i].Re)
       else AY.Ptr(1).arr^[i]:=0.5*pi*sign_(X^[i].Im)
    end else begin
     if OutMode^=1 then K:=K*(N-1)/sqrt((D1-sqr(CM1)/N)*(D2-sqr(CM2)/N));
     for i:=0 to AY.Ptr(0).Count-1 do begin
      AY.Ptr(0).arr^[i]:=f0*(i+1);
      AY.Ptr(1).arr^[i]:=K*sqrt(sqr(X^[i+1].Re)+sqr(X^[i+1].Im));
     end;
     Move(AY.Ptr(1).arr^,REZ.arr^,REZ.Count*SOfR);
    end;
    end;
   if (Position=Size^-1) then Position:=0 else inc(Position);
   time:=time+Tau^;
   end;
  end; //END Case
end;   //END Procedure

function TStatDblSpctrRst;
begin
 Result:=0;
 with PDblSpectrProp(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^,Size^*SOfCplx);
                BlockRead(F,GXY^,Size^*SOfCplx);
                BlockRead(F,W.arr^,W.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,GXY^,Size^*SOfCplx);
                BlockWrite(F,W.arr^,W.Count*SOfR);
                BlockWrite(F,REZ.arr^,REZ.Count*SOfR);
               except
                Result:=er_WriteFile
               end;
 end
end;

end.