vtfar.m

用于模拟时变非平稳的ARMA过程

function [AMLEst, B0LEst, EstInfo]= vtfar(X, MMAX, LMAX, crit)
% function [AMLEst, B0LEst, EstInfo]= vtfar(X[, MMAX, LMAX, crit])
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% An automatic procedure to fit a VTFAR(M, L) model to a multivariate 
% nonstationary signal x. The maximum search space defaults to 
%   MMAX= min(8, floor(N/16)), 
%   LMAX= min(4, floor(N/16)), 
% The order estimation IC <crit> defaults to the AIC. 
% EstInfo= {SigHat, InfoCr, NrPara, PENALTY, AX};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   2;
LAR  =   1;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
F    =   4;
filename= sprintf('vtfar%02d%02d%02d%04d', MAR, LAR, F, N)
load(filename);
%-------------
alpha= 1/2;
beta = 1/2;
AML0= AML;
B0L0= B0L;

E0= randn(N, F);
X= vtfarma_gen(E0, AML, B0L, F, alpha);
MMAX= MAR+1;
LMAX= LAR+1;
crit= 'AIC';


% THE CHANNEL
clear;tfpm
alpha= 1/2;
load updown_h_rusk.mat
figure(1);clf;subplot(2, 1, 1);mesh(abs(h_I2))

N= 512;
J= 5;

MMIN= 1;
MMAX= 6;
LMIN= 0;
LMAX= 3;
BMIN= 1;
BMAX= 3;

crit= 'AIC';

for offset= 0:200:1800
   X= 1000*h_I2(offset+1:offset+N, 29:29+J);
   subplot(2, 1, 2);mesh(abs(X/1000));drawnow
   [AML, B0L]= vtfar20(X, MMIN, MMAX, LMIN, LMAX, BMIN, BMAX, crit);
   [Mest, Lest, Best]
   filename= sprintf('ident%02d%02d%02d-%04d', J, BMIN, BMAX, offset)
   save(filename, 'AML', 'B0L', 'InfoCr', 'SigHat', 'Mest', 'Lest', 'Best');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some constants, etc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 1/2;      % TF shift parameter
beta  = 1/2;
N     = size(X, 1);
lambda= .9800;    % maximum root magnitude
rho   = 1-log(12);% for the MDL

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The nargin sory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargin<=3)
   crit= 'AIC';
end;

if(nargin<=2)
   LMAX= min(8, floor(N/16));
end;

if(nargin<=1)
   MMAX= min(8, floor(N/16));
   LMAX= min(4, floor(N/16));
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The penalty
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch(lower(crit))
 case 'mdl'
  PENALTY= log(N+1)+rho;
 case 'aic'
  PENALTY= 2;
 otherwise
  PENALTY= str2num(crit);
end;

[N, F]= size(X);
alpha= 1/2;

% Compute the EXAFs
%X= X/sqrt(energy(X));
[AX, ambi]= vambi_est(X, MMAX, LMAX);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:F
   for jp= 1:F
      [j jp]
      figure(1);clf;tf_show(reshape(AX(j, jp, :, :), N, N))
      figure(2);clf;mesh(abs(reshape(AX(j, jp, :, :), N, N)))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Prepare
InfoCr= 9999999*ones(MMAX, LMAX+1);
NrPara= InfoCr;
SigHat= InfoCr;
for m= 1:MMAX
   for l= 0:LMAX
      NrPara(m, l+1)= ((m+1)*(2*l+1)-1)*F^2;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VTFAR Identification 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to VTFAR(%d, %d) using the "%s"', MMAX, LMAX, crit)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for L= 0:LMAX
   for M= 1:MMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%for L= 1:1
%   for M= 2:2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%      [M L]

      [AML, B0L]= vtfar_parest(ambi, M, L, N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[AML0 AML]
sum(sum(sum(sum(abs(AML-AML0)))))/sum(sum(sum(sum(abs(AML0)))))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the MoNiC Model!!!!!!!!!!!!!
      b0L= zeros(size(B0L));
      bbb= inv(B0L(:, :, L+1, 1));
      for l= -L:L
         b0L(:, :, L+1+l, 1)= B0L(:, :, L+1+l, 1)*bbb;
      end;

      Ehat= vtfarma_inv(X, AML, b0L, F, alpha);
      V= Ehat'*Ehat/N;
      SigHat(M, L+1)= log(real(det(V)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for M= 1:MMAX;
end;%for L= 0:LMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

InfoCr= SigHat + PENALTY*NrPara/N;
[Mest, Lest]= find(InfoCr(:, :)==min(min(InfoCr(:, :))));
Lest= Lest-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[MAR LAR Mest Lest]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[AMLEst, B0LEst]= vtfar_parest(ambi, Mest, Lest, N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Further Detailled Info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargout==3)
   EstInfo= {SigHat, InfoCr, NrPara, PENALTY, AX};
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   2;
LAR  =   1;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
F    =   4;
filename= sprintf('vtfar%02d%02d%02d%04d', MAR, LAR, F, N)
load(filename);
%-------------
alpha= 1/2;
beta = 1/2;
AML0= AML;
B0L0= B0L;

MMAX= MAR+1;
LMAX= LAR+1;
crit= 'AIC';


MM= 100;
MLest= [];
for mm= 1:MM
   X= vtfarma_gen(randn(N, F), AML, B0L, F, alpha);
   [AMLEst, B0LEst]= vtfar(X);
   [Mest, Lest, F]= param_dim(AMLEst);
   [mm Mest Lest MAR LAR]
   MLest= [MLest; [Mest Lest]];
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%