{+doc
+NOMBRE: umatricesruida
+CREACION: 20180624
+AUTORES: rch
+REGISTRO:
+TIPO: Unidad Pascal.
+PROPOSITO: definición clase TMatricesRuida.
+PROYECTO: SimSEE
+REVISION:
+AUTOR:
+DESCRIPCION:
  Intenta definir clase para calcular las matrices de inyección de ruido
en las transformaciones de ampliación de un CEGH.
   El código se basa en lo desarrollado por Federico Barreto en 2016-2017

-doc}
unit umatricesruidas;


{$mode delphi}

interface

uses
  Classes, SysUtils, matreal;

type
  { TMatricesRuida }

  TMatricesRuida = class
  private

    // Filtro:     X_k+1 = A * X_k + B * R_k
    // Reduccion:  Z_k = M_r * X_k
    // Ampliacion: X_k = M_a * Z_k + B_a * W_a

    // Matrices del filtro
    A, B: TMatR;
    // Matrices de reduccion y ampliacion
    Mr, Ma: TMatR;

    MatrizRuidaMultiRetardos: TMatR;

    //*********** AUXILIARES ********
    // A'
    At,
    // M_a * M_r
    MaMr,
    // (M_a * M_r)'
    MaMr_t,
    // B * B'
    BBt,
    // Estado actual y estado siguiente
    XXt, XsXst: TMatR;

    //***********

    BampEE: TMatR;

    lst: TList;

  public

    (*

    *)
    constructor Create(A, B, MReductora, MAmplificadora: TMatR);

    procedure calcularse;

    // Obtiene la matriz ruida del paso k a partir de la matriz
    // <X_k-1 * X_k-1'>
    // True cuando el metodo convergio
    function Matriz_RUIDA_k(lambda_k: TMatR): boolean;

    function Get_Bamp(kPaso: integer): TMatR;

    function Count: integer;
  end;

const
  ERROR = 1.0e-5;


implementation

uses
  umatriz_ruida;

{ TMatricesRuida }

function TMatricesRuida.Get_Bamp(kPaso: integer): TMatR;
begin

  if A.nf <> A.nc then
    Result := MatrizRuidaMultiRetardos
  else
  begin
    if lst.Count <= kPaso then
      Result := lst.Last
    else
      Result := lst[kPaso];
  end;
end;

function TMatricesRuida.Count: integer;
begin
  Result := lst.Count;
end;

constructor TMatricesRuida.Create(A, B, MReductora, MAmplificadora: TMatR);
var
  Bt: TMatR;
begin

  lst := TList.Create;

  self.A := A;
  self.B := B;
  self.Ma := MAmplificadora;
  self.Mr := MReductora;


  At := TMatR.Create_Clone(A);
  At.Transponer;
  if A.nf <> A.nc then
  begin
  {$IFDEF RUIDA_MULTI_RETARDOS}
    MatrizRuidaMultiRetardos := Matriz_RUIDA_MultiRetardos(A, B, Ma, Mr);
  {$ELSE}
    MatrizRuidaMultiRetardos := nil;
  {$ENDIF}
  end;

  BampEE := Matriz_RUIDA_EE(A, B, Ma, Mr);

  //XXt inicialmente es una matriz nula
  XXt := TMatR.Create_Init(A.nf, A.nf);
  XXt.Ceros;

  if Mr = nil then
  begin
    MaMr := nil;
    MaMr_t := nil;
  end
  else
  begin
    MaMr := TMatR.Create_Init(Ma.nf, Mr.nc);
    MaMr.Mult(Ma, Mr);
    MaMr_t := TMatR.Create_Clone(MaMr);
    MaMr_t.Transponer;
  end;

  Bt := TMatR.Create_Clone(B);
  Bt.Transponer;
  BBt := TMatR.Create_Init(B.nf, B.nf);
  BBt.Mult(B, Bt);
  Bt.Free;

end;

procedure TMatricesRuida.calcularse;
begin
  i:=0;
  repeat
    i:=i+1;
    lambda:=GetLambda(i);
    conv:=Matriz_RUIDA_k(lambda);

    *****

  until conv or (max(globs.calcNPasosSim,globs.calcNPasosOpt)=i);


end;

{

  Filtro:       X_k+1 = A * X_k + B * R_k               Ec. 1
  Ampliacion:   X_k+1 = M_a * Z_k+1 + B_a * W_a         Ec. 2

  Aplicando la reduccion al filtro tenemos:

  M_r * X_k+1 = M_r * A * X_k + M_r * B * R_k
  Z_k+1 = M_r * A * X_k + M_r * B * R_k sustituyendo en la Ec. 2

  X_k+1 = M_a * [ M_r * A * X_k + M_r * B * R_k ] + B_a * W_a

  X_k+1 = M_a * M_r * A * X_k + M_a * M_r * B * R_k + B_a * W_a

  Multiplicando por el transpuesto y tomando valor esperado y considerando que
  X_k, R_k y W_k no estan correlacionados:

  <X_k+1 * X_k+1'> = M_a * M_r * A * <X_k * X_k'> * A' * M_r' * M_a' +
                     M_a * M_r * B * B' * M_r' * M_a' +                    Ec. 3
                     B_a * B_a'


  Multiplacando la Ec.1 por su transpuesto y tomando valor esperado tenemos:

  <X_k+1 * X_k+1'> = A * <X_k * X_k'> * A' + B * B'                        Ec. 4


  Despejando B_a * B_a' de Ec. 3 y Ec. 4:

  B_a * B_a' = A * <X_k * X_k'> * A' + B * B' -
               M_a * M_r * A * <X_k * X_k'> * A' * M_r' * M_a' -
               M_a * M_r * B * B' * M_r' * M_a'
}

function TMatricesRuida.Matriz_RUIDA_k(lambda_k: TMatR): boolean;
var
  XXtAt, AXXtAt, BBtLambda, LambdaBBtLambda, H: matreal.TMatR;
  MrH, Fa, Bamp_k_menos_uno: TMatR;

  Bamp_k: TMatR;
  dif: TMatR;
  j, k: integer;

  dimRes: integer;

begin

  XXtAt := TMatR.Create_Init(XXt.nf, At.nc);
  XXtAt.Mult(XXt, At);
  AXXtAt := TMatR.Create_Init(A.nf, XXtAt.nc);
  AXXtAt.Mult(A, XXtAt);

  XXtAt.Free;

  BBtLambda := TMatR.Create_Init(BBt.nf, lambda_k.nc);
  BBtLambda.Mult(BBt, lambda_k);

  LambdaBBtLambda := TMatR.Create_Init(lambda_k.nf, BBtLambda.nc);
  LambdaBBtLambda.Mult(lambda_k, BBtLambda);

  BBtLambda.Free;

  XsXst := TMatR.Create_Init(AXXtAt.nf, AXXtAt.nc);
  XsXst.Suma(AXXtAt, LambdaBBtLambda);

  AXXtAt.Free;
  LambdaBBtLambda.Free;

  //XsXst.WriteArchiXLT('XsXst.xlt');

  H := XXt.RaizPorPotenciaIterada(dimRes, False);

  if (H <> nil) and (dimRes > 0) then
  begin
    MrH := TMatR.Create_Init(Mr.nf, H.nc);
    MrH.Mult(Mr, H);

    MrH.CalcBSEN(Fa);

    Fa.OrtoNormal;
    Fa.Transponer;

    Bamp_k := TMatR.Create_Init(H.nf, Fa.nc);
    Bamp_k.Mult(H, Fa);

    H.Free;
    MrH.Free;
    Fa.Free;

    //Bamp_k.WriteArchiXLT('Ba.xlt');
  end
  else
    Bamp_k := nil;

  Bamp_k_menos_uno := TMatR(lst.Last);
  if (Bamp_k_menos_uno <> nil) and (Bamp_k <> nil) and (lst.Count > 1) then
  begin
    dif := TMatR.Create_Init(Bamp_k.nf, Bamp_k.nc);
    for k := 1 to Bamp_k.nf do
      for j := 1 to Bamp_k.nc do
        dif.pon_e(k, j, Bamp_k.e(k, j) - Bamp_k_menos_uno.e(k, j));

    Result := dif.NormMaxAbs <= ERROR;
  end
  else
    // No convergio
    Result := False;

  lst.Add(Bamp_k);

  XXt.Free;
  XXt := XsXst;

end;

end.



