vtfar20.m

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

function [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMIN, MMAX, LMIN, LMAX, BMIN, BMAX, crit)
% function [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMIN, MMMAX, LMIN, LMAX, BMIN, BMAX, crit)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% An automatic procedure to identify VTFAR(M, L) models. The
% maximum search space defaults to MMAX, LMAX= N/16, The order
% estimation IC defaults to the AIC. The signal is NOT normalized
% to energy N! 
% EstInfo= {SigHat, InfoCr, NrIter, NrPara, EstInn, PENALTY, AMBI};
%
% TODO: 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N=  64;
M0=  3;
L0=  2;
J0=  1;
B0=  1;

MMIN= 1;
MMAX= M0+1;
LMIN= 0;
LMAX= L0+1;
BMIN= 1;
BMAX= J0;
crit= 'AIC';

filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)
load(filename);
AML0= AML;
B0L0= B0L;


%AML0= eye(J0);B0L0= eye(J0);
E0= randn(N, J0);
X= vtfarma_gen(E0, AML0, B0L0, B0, 1/2);
L= L0;
M= M0;

[AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMAX, LMIN, LMAX, BMIN, BMAX, crit)




y= tfarma_gen(E0, reshape(AML, 2*L+1, M+1), reshape(B0L, 2*L+1, 1), 1/2);
norm(X-y)/norm(X)





% 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, Mest, Lest, Best, InfoCr, Varian]= 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', 'Varian', 'Mest', 'Lest', 'Best');
end;

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

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

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

% Compute the EXAFs
X= X/sqrt(energy(X));
[AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      [j jp]
      figure(1);clf;tf_show(reshape(AMBI(j, jp, :, :), N, N))
      figure(2);clf;mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))
      pause
   end;
end;

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

% Prepare
InfoCr= 9999999*ones(MMAX, LMAX+1, J);
NrPara= InfoCr;
Varian= InfoCr;
for m= 1:MMAX
   for l= 0:LMAX
      for b= 1:J
         NrPara(m, l+1, b)= ((m+1)*(.7*2*l+1)-1)*number_band(J, b);
      end;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for L= LMIN:LMAX
   for M= MMIN: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;

      for B= BMIN:BMAX
%         clear AMBI
%	 pack
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         clear Linv
         Linv= vtfarma_Linv(AML, b0L, N, B, alpha);
         V= zeros(J);
         for tau= 1:J
            for taup= 1:J
               for t= 1:J
                  for tp= 1:J
%                     [M L B tau taup t tp]
                     V(tau, taup)= V(tau, taup) + sum(sum(...
                     reshape(Linv(tau, t, :, :), N, N).*...
                     reshape(conj(Linv(taup, tp, :, :)), N, N).*...
                     nm_to_nk(ml_to_nm(reshape(AMBI(t, tp, :, :), N, N)))));
                  end;
               end;
            end;
         end;
         V= V/N^2;
         e= real(eig(V));
         e= e.*(e>0);
	 if(prod(e))
	    Varian(M, L+1, B)= log(prod(e));
	 else
	    Varian(M, L+1, B)= 99999999;
	 end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         Ehat= vtfarma_inv(X, AML, b0L, B, 1/2);
         V= Ehat'*Ehat/N;
	 Varian(M, L+1, B)= log(real(det(V)));
      end;%for B= 0:J-1

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

InfoCr= Varian + PENALTY*NrPara/N;
MinCr= max(max(max(InfoCr)));
Best= 0;
for b= 1:J
   if(min(min(InfoCr(:, :, b)))<MinCr)
      [Mest, Lest]= find(InfoCr(:, :, b)==min(min(InfoCr(:, :, b))));
      MinCr= min(min(InfoCr(:, :, b)));
      Best= b;
   end;
end;
Lest= Lest-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M0 L0 B0 Mest Lest Best]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[AML, B0L]= vtfar_parest(ambi, Mest, Lest, N);
for j= 1:J
   for jp= 1:J
      if(abs(j-jp)>=Best)
         AML(j, jp, :, :)= zeros(2*Lest+1, Mest+1);
         B0L(j, jp, :, :)= zeros(2*Lest+1, 1);
      end;
   end;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N=  256;
M0= 2;
L0= 1;
J0= 4;
B0= 2;

MMAX= M0+2;
LMAX= L0+1;
crit= 'AIC';

filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)
load(filename);
AML0= AML;
B0L0= B0L;

MM= 100;
MLBest= [];
for mm= 1:MM
   X= vtfarma_gen(randn(N, J0), AML0, B0L0, B0, 1/2);
   [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMAX, 0, LMAX, 1, J0, crit);
   [mm Mest Lest Best M0 L0 B0]
   MLBest= [MLBest; [Mest Lest Best]];
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%