tfar_fit_aic.m

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

function [A, B]= tfar_fit_aic(x, Mmax, Lmax, Psi)
% function [AA, BB]= tfar_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, 3*Mmax, 4*Lmax, 0, 2, 1);
Psi= Psi(:, N/2-Mmax+1:N/2+Mmax+1);

Qalpha= corr_est(x, x, Mmax, alpha);
Aalpha= fft(Qalpha);
Aalpha= [ Aalpha(N/2+1:N, :); Aalpha(1:N/2, :)].*conj(Psi);
E= zeros(Mmax, Lmax);

for LAR= 1:Lmax
   Ayy= Aalpha(N/2-3*LAR+1:N/2+3*LAR+1, :);
   [AA, BB]= tfar_est_tfywu(Ayy, N);
   for MAR= 1:Mmax
      e= real(tfarma_inv(x, AA(:, 1:1+MAR, MAR), BB(:, :, MAR)));
%      figure(1);clf;plot(e)
%      figure(2);clf;tf_show(tfarma_wvsp(AA(:, 1:1+MAR, MAR), BB(:, :, MAR), N, alpha));
      drawnow
      [MAR LAR var(e)];
      E(MAR, LAR)= var(e);
   end;
end;
AIC= N*log(E) + 4*(1:MAR)'*(2:1+LAR);
[MAR, LAR]= find(AIC==min(min(AIC)));

Qalpha= corr_est(x, x, MAR, alpha);
Aalpha= fft(Qalpha);
Aalpha= [ Aalpha(N/2+1:N, :); Aalpha(1:N/2, :)];
Ayy= Aalpha(N/2-3*LAR+1:N/2+3*LAR+1, :);
[AA, BB]= tfar_est_tfywu(Ayy, N);
A= AA(:, 1:MAR+1, MAR);
B= BB(:, :, MAR);