tfarma1.m

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

function [MARhat, LARhat, MMAhat, LMAhat]= tfarma1(x)
% function [MARhat, LARhat, MMAhat, LMAhat]= tfarma1(x)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% An automatic procedure to identify (non)stationary (TF)(AR)(MA)
% models. 
%TODO: Write LOGGING!!!

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

x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
figure(1);clf;mesh(tfarma_wvsp(Aml, Bml, N, alpha));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha= 1/2;
beta = 1/2;
N= length(x);
lambda= .9800;
N= length(x);
Lmax= N/8;

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

%Ax= tfarma_ambi(Aml, Bml, N, alpha);

figure(2);clf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TFAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigmahatTFAR= 9999999999*ones(N/2, Lmax+1);
nriteratTFAR= zeros(N/2, Lmax+1);
lambdamxTFAR= zeros(N/2, Lmax+1);
NTFAR= (1:(N/2))'*(2*(0:Lmax)+1)-1;
TFAR= {};
for L= 0:Lmax
   if(~L)
      Mmax= N/2-1;
   else
      Mmax= min(floor(N/4/L), N/8);
   end;%if(~L)
   [L Mmax];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
      M;
      if(L)
	 [Amlstab, lambdamax, mm]= param_stabilize(ATFAR(:, 1:M+1, M), N, lambda, 1);
         ATFAR(:, 1:M+1, M)= Amlstab;
         nriteratTFAR(M+1, L+1)= mm;
         lambdamxTFAR(M+1, L+1)= lambdamax;
      end;
      e= tfarma_inv(x, ATFAR(:, 1:M+1, M), BTFAR(:, 1, M)/param_get(BTFAR(:, 1, M), 0, 0));
      sigmahatTFAR(M+1, L+1)= e'*e/N;
%      var(tfarma_inv(x, ATFAR(:, 1:M+1, M), BTFAR(:, 1, M)));
%      eprime= tfarma_inv(x, Aml1, Bml1/param_get(Bml1, 0, 0));
%      E(M, L+1, 1)     = e'*e/N;
%      Eprime(M, L+1, 1)= eprime'*eprime/N;
   end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% PLOT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   if(L)
      for M= 1:Mmax
         subplot(Lmax, N/4, N/4*(L-1)+M)
         tf_show(tfarma_wvsp(ATFAR(:, 1:M+1, M), BTFAR(:, 1, M), N, alpha))
         drawnow;
      end;
   end;

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

%   TFAR= {TFAR {ATFAR, BTFAR}};
end;%for L= 0:Lmax

SigmahatTFAR= log(sigmahatTFAR);
AIC_TFAR= SigmahatTFAR + 2*NTFAR/N;
MDL_TFAR= SigmahatTFAR + (log(N)-1.5)*NTFAR/N;

%AIC_TFAR==min(min(AIC_TFAR));
%MDL_TFAR==min(min(MDL_TFAR));

[MARhat, LARhat] = find(MDL_TFAR==min(min(MDL_TFAR)));
LARhat= LARhat-1;
MARhat= MARhat-1;
MMAhat= 0;
LMAhat= 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   4;
LAR  =   3;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
lambda= .9800;
Bml= Bml/param_get(Bml, 0, 0);
MM= 100;

for mm= 1:MM
   x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
   [MARhat, LARhat, MMAhat, LMAhat]= tfarma1(x);
   [MARhat, LARhat, MMAhat, LMAhat]
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%