tfarma.m

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

function [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, crit)
% function [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(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 identify TFARMA(M, L) models. The
% maximum TFAR search space defaults to MMAX, LMAX= N/8, The order
% estimation IC defaults to the AIC. 
% TODO: M/LMAX even

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

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

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



[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, '1.5');
[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX);
[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX);
[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x);

%nargin= 4;

MMAX= 5;
LMAX= 5;
MM= 100;
NRITER= zeros(MMAX, LMAX+1);
for mm= 1:MM
   x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
   [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, 'MDL');
   [MA, LA]= param_dim(AmlEst);
   [MB, LB]= param_dim(BmlEst);
   [mm MA LA MB LB]
   NRITER= NRITER + NrIter(:, :, 1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

if(nargin<=1)
   MMAX= floor(N/8);
   LMAX= floor(N/8);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some constants, etc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 1/2;
beta  = 1/2;
N     = length(x);
lambda= .9800;
rho   = 1-log(12);

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

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

SigmaHat= 9999999999*ones(MMAX, LMAX+1, 3);
NrIter  = zeros(MMAX, LMAX+1, 3);

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFAR Identification (1:Mmax)x(0:LMAX)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFAR(%d, %d) using the "%s"', MMAX, LMAX, crit)
TFAR= {};
TFARunstab= {};
for L= 0:LMAX

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% 127/32/32
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(~L)
      Mmax= N/2-1;
   else
      Mmax= min(floor(N/4/L), MMAX);
   end;%if(~L)
   [L Mmax]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   warning off MATLAB:divideByZero
   Mmax= min(floor(N/4/L), 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   [ATFAR, BTFAR]= tfar_est_tfywu(Ax(N/2+1-3*L:N/2+1+3*L, N/2-Mmax+1:N/2+Mmax+1), N);
   for M= 1:Mmax
      [Aml1, lambdamax, mmTFAR, nrFDIR]= param_stabilize(ATFAR(:, 1:M+1, M), N, lambda);
      [M L mmTFAR];
      B= BTFAR(:, :, M);
      TFAR{M, L+1}= {Aml1 B};
      TFARunstab{M, L+1}= {Aml1 B};
      eprime= tfarma_inv(x, Aml1, B/param_get(B, 0, 0));
      SigmaHat(M, L+1, 1)= eprime'*eprime/N;
      NrIter(M, L+1, 1)= mmTFAR;
   end;
end;%for L= 0:LMAX
NATFAR= (1:MMAX)'*(2*(0:LMAX)+1);
NBTFAR= ones(MMAX, 1)*(2*(0:LMAX)+1);
ICTFAR= log(SigmaHat(:, :, 1)) + PENALTY*(NATFAR+NBTFAR)/N;
[MTFAR, LTFAR] = find(ICTFAR==min(min(ICTFAR)));LTFAR= LTFAR-1;
AmlTFAR= TFAR{MTFAR, LTFAR+1}{1};
BmlTFAR= TFAR{MTFAR, LTFAR+1}{2};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFMA Identification (1:Mmax/2)x(0:LMAX/2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFMA(%d, %d) using the "%s"', MMAX/2, LMAX/2, crit)
TFMA= {};
TFMAunstab= {};
for L= 0:LMAX/2
   warning off MATLAB:divideByZero
   Mmax= min(floor(N/8/L), MMAX/2);
   warning on MATLAB:divideByZero
   for M= 1:Mmax
      BB= tfma_est_dmu(TFAR{2*M, 2*L+1}{1}, TFAR{2*M, 2*L+1}{2}, M, L, N);
      [Bml1, lambdamax, mmTFMA, nrFDIR]= param_stabilize(BB(:, 1:M+1, M+1), N, lambda);
      [M L mmTFMA];
      TFMA{M, L+1}= Bml1;
      TFMAunstab{M, L+1}= BB(:, 1:M+1, M+1);      
      eprime= tfarma_inv(x, 1, Bml1/param_get(Bml1, 0, 0));
      SigmaHat(M, L+1, 2)= eprime'*eprime/N;
      NrIter(M, L+1, 2)= mmTFMA;
   end;
end;
NATFMA= zeros(MMAX/2, 1)*(2*(0:LMAX/2)+1);
NBTFMA= (1:MMAX/2)'*(2*(0:LMAX/2)+1);
ICTFMA= log(SigmaHat(1:Mmax, 1:LMAX/2+1, 2)) + PENALTY*(NATFMA+NBTFMA)/N;
[MTFMA, LTFMA] = find(ICTFMA==min(min(ICTFMA)));LTFMA= LTFMA-1;
AmlTFMA= [zeros(LTFMA, 1); 1; zeros(LTFMA, 1)];
BmlTFMA= TFMA{MTFMA, LTFMA+1};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFARMA Identification (1:Mmax/2)x(0:LMAX/2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Searching up to TFARMA(%d, %d) using the "%s"', MMAX/2, LMAX/2, crit)
TFARMA= {};
TFARMAunstab= {};
for L= 0:LMAX/2
   warning off MATLAB:divideByZero
   Mmax= min(floor(N/8/L), MMAX/2);
   warning on MATLAB:divideByZero
   for M= 1:Mmax
      [AA, BB]= tfarma_est_gkmu(TFAR{2*M, 2*L+1}{1}, TFAR{2*M, 2*L+1}{2}, M, L, M-1, L, N);
      [Bml1, lambdamax, mmTFMA, nrFDIR]= param_stabilize(BB(:, 1:M+1, M+1), N, lambda);
      [M L mmTFMA];
      TFMA{M, L+1}= Bml1;
      TFMAunstab{M, L+1}= BB(:, 1:M+1, M+1);      
      eprime= tfarma_inv(x, 1, Bml1/param_get(Bml1, 0, 0));
      SigmaHat(M, L+1, 2)= eprime'*eprime/N;
      NrIter(M, L+1, 2)= mmTFMA;
   end;
end;
NATFMA= zeros(MMAX/2, 1)*(2*(0:LMAX/2)+1);
NBTFMA= (1:MMAX/2)'*(2*(0:LMAX/2)+1);
ICTFMA= log(SigmaHat(1:Mmax, 1:LMAX/2+1, 2)) + PENALTY*(NATFMA+NBTFMA)/N;
[MTFMA, LTFMA] = find(ICTFMA==min(min(ICTFMA)));LTFMA= LTFMA-1;
AmlTFMA= [zeros(LTFMA, 1); 1; zeros(LTFMA, 1)];
BmlTFMA= TFMA{MTFMA, LTFMA+1};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model Type Selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%P0= tfarma_wvsp(Aml, Bml, N, alpha);
P1= tfarma_wvsp(AmlTFAR, BmlTFAR, N, alpha);
P2= tfarma_wvsp(AmlTFMA, BmlTFMA, N, alpha);
%figure(1);tf_show(P0)
figure(2);tf_show(P1)
figure(3);tf_show(P2)

if(ICTFAR(MTFAR, LTFAR+1)<ICTFMA(MTFMA, LTFMA+1))
   sprintf('TFAR!')
   AmlEst= AmlTFAR;
   BmlEst= BmlTFAR;
else
   sprintf('TFMA!')
   AmlEst= AmlTFMA;
   BmlEst= BmlTFMA;
end;