tfar.m

来自「用于模拟时变非平稳的ARMA过程」· M 代码 · 共 260 行

function [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, crit)
% function [AmlEst, B0lEst, EstInfo]= tfar(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 TFAR(M, L) model to a non- 
% stationary 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)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   3;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;

x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
MMAX= 2*MAR;
LMAX= 2*LAR;
crit= 'AIC';

P0= tfarma_wvsp(Aml, Bml, N, alpha);
figure(1);tf_show(P0);drawnow
%nargin= 1;

% Some call examples 
[AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, '1.5');
[AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX);
[AmlEst, B0lEst, EstInfo]= tfar(x, MMAX);
[AmlEst, B0lEst]= tfar(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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<=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;

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

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TAPER DESIGN (NOT ACTIVE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   [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+1-Mmax:N/2+1+Mmax), N);
   for M= 1:Mmax
      B= BAR(:, :, M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if(0)% VARIANCE COMPUTATION BY INVERSE FILTERING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [Aml1, lambdamax, mmAR, nrFDIR]= param_stabilize(AAR(:, 1:M+1, M), N, lambda);
      [M L mmAR];
      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)= eprime'*eprime/N;
      NrIter(M+1, L+1)= mmAR;
      EstInn(M+1, L+1, :)= eprime;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% VARIANCE COMPUTATION BY INNER PRODUCT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      AR{M+1, L+1}      = {AAR(:, 1:M+1, M) B};
      Lalpha= tfarma_weyl(AAR(:, 1:M+1, M), B/param_get(B, 0, 0), N, alpha);
      SigHat(M+1, L+1)= real(sum(sum(conj(Rx)./abs(Lalpha)^2)))/N;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;
end;%for L= 0:LMAX
NrPara(:, :)= (1:MMAX+1)'*(2*(0:LMAX)+1) - 1;
InfoCr(:, :)= log(SigHat(:, :)) + PENALTY*NrPara(:, :)/N;
[MAREst, LAREst] = find(InfoCr(:, :)==min(min(InfoCr(:, :))));
MAREst= MAREst-1;LAREst= LAREst-1;
AmlEst= param_stabilize(AR{MAREst+1, LAREst+1}{1}, N, lambda);
B0lEst= AR{MAREst+1, LAREst+1}{2};

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   3;
LAR  =   2;
MMA  =   0;
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
   mm
   x= tfarma_gen(randn(N, 1), Aml, Bml);

   [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, 'MDL');
   STATUS{1, mm}= EstInfo;

   [AmlEst, B0lEst, EstInfo]= tfar(x, MMAX, LMAX, 'AIC');
   STATUS{2, mm}= EstInfo;
end;
filename= sprintf('tfar_%04d%d%d%s%s', N, MAR, LAR, 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  =   3;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%