tfar_fit.m

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

function [TFFPE, TFAIC]= tfar_fit(x, Mmax, Lmax, Psi, lambda)
% function [TFFPE, TFAIC]= tfar_fit(x, Mmax, Lmax, Psi, lambda))
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Searches M= 1:Mmax, L= 0:Lmax. 

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

Mmax= 2*MAR;
Lmax= 2*LAR;
Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);
x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha= 1/2;
beta = 1/2;
N= length(x);

Qalpha= corr_est(x, x, -1, alpha);
Aalpha= fft(Qalpha);
Aalpha= [ Aalpha(N/2+1:N, :); Aalpha(1:N/2, :)].*conj(Psi);
E= zeros(Mmax, Lmax+1);
NA= (1:Mmax)'*(2*(0:Lmax)+1) + ones(Mmax, 1)*2*(0:Lmax);
for L= 0:Lmax
   [AA, BB]= tfar_est_tfywu(Aalpha(N/2-3*L+1:N/2+3*L+1, N/2-Mmax+1:N/2+Mmax+1), N);

% Inverse Filtering
   for M= 1:Mmax
      [Aml, lambdamax, mm]= param_stabilize(AA(:, 1:M+1, M), N, lambda);
%      [M L mm]
      eM= tfarma_inv(x, Aml, BB(:, :, M));
%      figure(99);plot(real(eM));drawnow
      rM= sum(eM'*eM)/N;
      E(M, L+1)= rM;
   end;
end;

e= log(E);

TFFPE= e+log((N+NA+1)./(N-NA-1));
TFAIC= e+2*NA/N;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
MM= 100;
mm= 1
Mmax= 2*MAR+2;
Lmax= 2*LAR+2;
[Psi, Mask, v2]= tf_multiwin(N, 2*Mmax, 2*Lmax, 5, 0, 0);
v2
tf_show(Psi)
AIChist= zeros(Mmax, Lmax+1);
FPEhist= zeros(Mmax, Lmax+1);
for mm= 1:MM
   x= tfarma_gen(randn(N, 1), Aml, Bml, 1/2);
   [TFFPE, TFAIC]= tfar_fit(x, Mmax, Lmax, Psi, .98);
   [MFPE, LFPE] = find(TFFPE==min(min(TFFPE)));
   [MAIC, LAIC] = find(TFAIC==min(min(TFAIC)));
   [MFPE, LFPE, MAIC, LAIC]
   AIChist(MAIC, LAIC)= AIChist(MAIC, LAIC) + 1;
   FPEhist(MFPE, LFPE)= FPEhist(MFPE, LFPE) + 1;
   figure(1);imagesc(AIChist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(2);imagesc(FPEhist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   drawnow
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%