tfarma20.m

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

function [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, crit, longord)
% function [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, crit, longord)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% An automatic procedure to identify TFAR(M, L), TFMA(M, L), or
% monic TFARMA(M, L) models. The maximum search space defaults to
% MMAX, LMAX= N/16, The order estimation IC defaults to the
% AIC. The long order factor defaults to 2. The signal is NOT 
% normalized to energy N! The TFARMA IC is decreased by log(2)!
% EstInfo= {SigHat, InfoCr, NrIter, NrPara, EstInn, AmlAR, ... 
%          BmlAR, AmlMA, BmlMA, AmlAA, BmlAA, PENALTY, Ax};

% TODO: 
% monic TFARMA

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MAR  =   3;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;

%Psi= ones(N);
%Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);

x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
longord= 3;
MMAX= max(4, longord*max(MAR, MMA));
LMAX= longord*max(LAR, LMA);
crit= 'AIC';

P0= tfarma_wvsp(Aml, Bml, N, alpha);
figure(1);tf_show(P0);drawnow
%[AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, '1.5', 2);
[AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, '1.5');
[AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX);
[AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX);
[AmlEst, BmlEst, EstInfo]= tfarma20(x);
%nargin= 4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The nargin sory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargin<=4)
   longord= 2;
end;

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

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

if(nargin<=1)
   MMAX= floor(N/16);
   LMAX= 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;

MMAX= min(MMAX, floor(N/16));
LMAX= min(LMAX, floor(N/16));

SigHat= exp(99)*ones(MMAX+1, LMAX+1, 3);
NrIter= zeros(MMAX+1, LMAX+1, 3);
NrPara= NrIter;
InfoCr= NrIter;
EstInn= zeros(MMAX+1, LMAX+1, 3, N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonparametric signal statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ax= fft(corr_est(x, x, longord*MMAX, 1/2));
Ax= [Ax(N/2+1:N, :); Ax(1:N/2, :)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFAR Identification 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFAR(%d, %d) using the "%s"', MMAX, LMAX, crit)
AR= {};
ARunstab= {};
for L= 0:longord*LMAX

   warning off MATLAB:divideByZero
   Mmax= min(max(longord*MMAX, floor(longord*N/16/L)), longord*MMAX);
   warning on MATLAB:divideByZero

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TAPER DESIGN (CANCEL???)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   [Psi, Mask, v2]= tf_multiwin(N, 3*Mmax, 4*L, 0, 2, 1);
   v2
   subplot(1, N/4+1, L+1)
   tf_show(Psi);drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   [AAR, BAR]= tfar_est_tfywu(Ax(N/2+1-3*L:N/2+1+3*L, N/2-longord*Mmax+1:N/2+longord*Mmax+1), N);
   [AAR, BAR]= tfar_est_tfywu(Ax(N/2+1-3*L:N/2+1+3*L, :), N);
   for M= 1:Mmax
      if( (M<=MMAX)&(L<=LMAX) )
	 [Aml1, lambdamax, mmAR, nrFDIR]= param_stabilize(AAR(:, 1:M+1, M), N, lambda);
         [M L mmAR];
         B= BAR(:, :, M);
         AR{M+1, L+1}      = {Aml1 B};
         eprime= tfarma_inv(x, Aml1, B/param_get(B, 0, 0));
         figure(99);subplot(2, 1, 1);plot(real(eprime));subplot(2, 1, 2);plot(imag(eprime));title(sprintf('TFAR(%d, %d); %d', M, L, mmAR));drawnow
         SigHat(M+1, L+1, 1)= eprime'*eprime/N;
         NrIter(M+1, L+1, 1)= mmAR;
	 EstInn(M+1, L+1, 1, :)= eprime;
      end;
      ARunstab{M+1, L+1}= {AAR(:, 1:M+1, M) BAR(:, :, M)};
   end;
end;%for L= 0:LMAX
NrPara(:, :, 1)= (1:MMAX+1)'*(2*(0:LMAX)+1) - 1;
InfoCr(:, :, 1)= log(SigHat(:, :, 1)) + PENALTY*NrPara(:, :, 1)/N;
[MAREst, LAREst] = find(InfoCr(:, :, 1)==min(min(InfoCr(:, :, 1))));
MAREst= MAREst-1;LAREst= LAREst-1;
AmlAR= AR{MAREst+1, LAREst+1}{1};
BmlAR= AR{MAREst+1, LAREst+1}{2};





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFMA Identification 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFMA(%d, %d) using the "%s" and alpha= %d', MMAX, LMAX, crit, longord)
MA= {};
MAunstab= {};
for L= 0:LMAX

%   warning off MATLAB:divideByZero
%   Mmax= min(floor(N/8/L), MMAX);
%   warning on MATLAB:divideByZero

   for M= 1:MMAX
      BB= tfma_est_dmu(ARunstab{longord*M+1, longord*L+1}{1}, ARunstab{longord*M+1, longord*L+1}{2}, M, L, N);
      [Bml1, lambdamax, mmMA, nrFDIR]= param_stabilize(BB(:, 1:M+1, M+1), N, lambda);
      [M L mmMA];
      MA{M+1, L+1}      = {1 Bml1};
      MAunstab{M+1, L+1}= {1 BB(:, 1:M+1, M+1)};      
      eprime= tfarma_inv(x, 1, Bml1/param_get(Bml1, 0, 0));
      figure(99);subplot(2, 1, 1);plot(real(eprime));subplot(2, 1, 2);plot(imag(eprime));title(sprintf('TFMA(%d, %d); %d', M, L, mmMA));drawnow
      SigHat(M+1, L+1, 2)= eprime'*eprime/N;
      NrIter(M+1, L+1, 2)= mmMA;
      EstInn(M+1, L+1, 2, :)= eprime;
   end;
end;
NrPara(:, :, 2)= (1:MMAX+1)'*(2*(0:LMAX)+1) - 1;
InfoCr(:, :, 2)= log(SigHat(:, :, 2)) + PENALTY*NrPara(:, :, 2)/N;
[MMAEst, LMAEst] = find(InfoCr(:, :, 2)==min(min(InfoCr(:, :, 2))));
MMAEst= MMAEst-1;LMAEst= LMAEst-1;
AmlMA= MA{MMAEst+1, LMAEst+1}{1};
BmlMA= MA{MMAEst+1, LMAEst+1}{2};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFARMA Identification 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFARMA(%d, %d) using the "%s" and alpha= %d', MMAX, LMAX, crit, longord)
ARMA= {};
AAunstab= {};
for L= 0:LMAX

   warning off MATLAB:divideByZero
   Mmax= min(floor(N/8/L), MMAX);
   warning on MATLAB:divideByZero

   for M= 2:MMAX
      [A, B]= tfarma_est_gkmu(ARunstab{longord*M+1, longord*L+1}{1}, ARunstab{longord*M+1, longord*L+1}{2}, M, L, M-1, L, N);
      [Aml1, lambdamax, mmAA1, nrFDIR]= param_stabilize(A, N, lambda);
      [Bml1, lambdamax, mmAA2, nrFDIR]= param_stabilize(B, N, lambda);
      [M L mmAA1 mmAA2];
      AA{M+1, L+1}      = {Aml1 Bml1};
      AAunstab{M+1, L+1}= {A B};
      eprime= tfarma_inv(x, Aml1, Bml1/param_get(Bml1, 0, 0));
      figure(99);subplot(2, 1, 1);plot(real(eprime));subplot(2, 1, 2);plot(imag(eprime));title(sprintf('TFARMA(%d, %d); %d', M, L, mmAA1+mmAA2));drawnow
      SigHat(M+1, L+1, 3)= eprime'*eprime/N;
      NrIter(M+1, L+1, 3)= mmAA1+mmAA2;
      EstInn(M+1, L+1, 3, :)= eprime;
   end;
end;
NrPara(:, :, 3)= 2*(0:MMAX)'*(2*(0:LMAX)+1) - 1 ;
InfoCr(:, :, 3)= log(SigHat(:, :, 3)) + PENALTY*NrPara(:, :, 3)/N;
[MAAEst, LAAEst] = find(InfoCr(:, :, 3)==min(min(InfoCr(:, :, 3))));
MAAEst= MAAEst-1;LAAEst= LAAEst-1;
AmlAA= AA{MAAEst+1, LAAEst+1}{1};
BmlAA= AA{MAAEst+1, LAAEst+1}{2};

%THAT'S IT:???
%InfoCr(:, :, 3)= InfoCr(:, :, 3) - log(1.5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model Type Selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%P0= tfarma_wvsp(Aml, Bml, N, alpha);
%PAR= tfarma_wvsp(AmlAR, BmlAR, N, alpha);
%PMA= tfarma_wvsp(AmlMA, BmlMA, N, alpha);
%PAA= tfarma_wvsp(AmlAA, BmlAA, N, alpha);
%figure(1);tf_show(P0)
%figure(90);tf_show(PAR)
%figure(91);tf_show(PMA)
%figure(92);tf_show(PAA)

if(InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MMAEst+1, LMAEst+1, 2) & ...
   InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MAAEst+1, LAAEst+1, 3)-log(2))
   sprintf('TFAR(%d, %d)', MAREst, LAREst)
   AmlEst= AmlAR;
   BmlEst= BmlAR;
else
   if(InfoCr(MMAEst+1, LMAEst+1, 2)<InfoCr(MAAEst+1, LAAEst+1, 3)-log(2))
      sprintf('TFMA(%d, %d)', MMAEst, LMAEst)
      AmlEst= AmlMA;
      BmlEst= BmlMA;
   else
      sprintf('TFARMA(%d, %d)', MAAEst, LAAEst)
      AmlEst= AmlAA;
      BmlEst= BmlAA ;     
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   sprintf('TFAR(%d, %d)', MAREst, LAREst)
   AmlEst= AmlAR;
   BmlEst= BmlAR;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Further Detailled Info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EstInfo= {SigHat, InfoCr, NrIter, NrPara, EstInn, AmlAR, BmlAR, AmlMA, BmlMA, AmlAA, BmlAA, PENALTY, Ax};
EstInfo= {SigHat, InfoCr, NrIter, NrPara, EstInn, PENALTY, Ax};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MAR  =   0;
LAR  =   2;
MMA  =   3;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;

MMAX= 6;
LMAX= 6;
MM= 20;
STATUS= {};
for mm= 1:MM
   x= tfarma_gen(randn(N, 1), Aml, Bml, beta);

   [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, 'MDL', 2);
   STATUS{1, 1, mm}= EstInfo;

   [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, 'MDL', 3);
   STATUS{1, 2, mm}= EstInfo;

   [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, 'AIC', 2);
   STATUS{2, 1, mm}= EstInfo;

   [AmlEst, BmlEst, EstInfo]= tfarma20(x, MMAX, LMAX, 'AIC', 3);
   STATUS{2, 2, mm}= EstInfo;
end;
filename= sprintf('tfarma20_%04d%d%d%d%d%s%s', N, MAR, LAR, MMA, LMA, re_im, mo_no)
save(filename, 'STATUS')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
mm= 1;
MM= 20;

filename= sprintf('tfarma20_%04d%d%d%d%d%s%s', N, MAR, LAR, MMA, LMA, re_im, mo_no)
load(filename)

TFAR= 0;
TFMA= 0;
TFAA= 0;
MAR= 0;LAR= 0;
MMA= 0;LMA= 0;
MAA= 0;LAA= 0;
for mm= 1:MM
   InfoCr= STATUS{2, 2, mm}{2};
   
   [MAREst, LAREst] = find(InfoCr(:, :, 1)==min(min(InfoCr(:, :, 1))));
   MAREst= MAREst-1;LAREst= LAREst-1;
   [MMAEst, LMAEst] = find(InfoCr(:, :, 2)==min(min(InfoCr(:, :, 2))));
   MMAEst= MMAEst-1;LMAEst= LMAEst-1;
   [MAAEst, LAAEst] = find(InfoCr(:, :, 3)==min(min(InfoCr(:, :, 3))));
   MAAEst= MAAEst-1;LAAEst= LAAEst-1;

   if(InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MMAEst+1, LMAEst+1, 2) & ...
      InfoCr(MAREst+1, LAREst+1, 1)<InfoCr(MAAEst+1, LAAEst+1, 3))
      sprintf('TFAR(%d, %d)', MAREst, LAREst)
      MAR= MAR+MAREst;LAR= LAR+LAREst;
      TFAR= TFAR+1;
   else
      if(InfoCr(MMAEst+1, LMAEst+1, 2)<InfoCr(MAAEst+1, LAAEst+1, 3))
         sprintf('TFMA(%d, %d)', MMAEst, LMAEst)
	 MMA= MMA+MMAEst;LMA= LMA+LMAEst;
         TFMA= TFMA+1;
      else
         sprintf('TFARMA(%d, %d)', MAAEst, LAAEst)
	 MAA= MAA+MAAEst;LAA= LAA+LAAEst;
         TFAA= TFAA+1;
      end;
   end;
end;
[100*TFAR/MM MAR/TFAR LAR/TFAR]
[100*TFMA/MM MMA/TFMA LMA/TFMA]
[100*TFAA/MM MAA/TFAA LAA/TFAA]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%