tfma_fit_aic.m

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

function B= tfma_fit_aic(x, Mmax, Lmax, Psi)
% function B= tfma_fit_aic(x, Mmax, Lmax, Psi)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
load('~/matlab/data/Motordaten/bmw1000.mat');
MM= 100;
N0= 1468;
N= 256;
offset= 64;
alpha= 1/2;

x= c1(offset:2:offset+2*N-1, 1)-mean(c1(offset:2:offset+2*N-1, 1));
plot(x)
Mmax= 7;
Lmax= 9;

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

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

%Psi= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);
%Psi= Psi(:, N/2-Mmax+1:N/2+Mmax+1);

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);

for LMA= 1:Lmax
   for MMA= 1:Mmax
      B= tfma_est_cepsb(Aalpha, MMA, LMA);
      e= real(tfarma_inv(x, 1, B));
      figure(1);clf;plot(e)
      figure(2);clf;tf_show(tfarma_wvsp(1, B, N, alpha));
      drawnow
      [MMA LMA var(e)];
      E(MMA, LMA)= var(e);
   end;
end;
AIC= N*log(E) + 4*(1:MMA)'*(2:1+LMA);
[MMA, LMA]= find(AIC==min(min(AIC)));

B= tfma_est_cepsb(Aalpha, MMA, LMA);