tfarma_fit.m

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

function [MDL, BIC, AIC, GIC, MIC, mdl, bic, aic, gic, mic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda)
% function [MDL, BIC, AIC, GIC, mdl, bic, aic, gic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Searches TFARMA(M, L, M-1, L) models for M= 1:Mmax, L= 0:Lmax. 

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

Mmax= max([2*MAR 2]);
Lmax= max([2*LAR 2]);
Psi= ones(N);%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);

Ax= fft(corr_est(x, x, -1, alpha));
Ax= [ Ax(N/2+1:N, :); Ax(1:N/2, :)].*conj(Psi);

E= zeros(Mmax, Lmax+1);
Eprime= zeros(Mmax, Lmax+1);

for L= 0:Lmax
   for M= 1:Mmax
      AA= tfarma_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2+(M-1)-M+1:N/2+(M-1)+M+1), M-1, N);

% Inverse TFAR Filtering
      [Aml1, lambdamax, mm]= param_stabilize(AA(:, 1:M+1, M), N, lambda);
      uM= tfarma_inv(x, Aml1, 1);
      Au= fft(corr_est(uM, uM, -1, alpha));
      Au= [ Au(N/2+1:N, :); Au(1:N/2, :) ].*conj(Psi);
      [B, REG, PMIN]= tfma_est_cepsb(Au, M-1, L);
      [Bml1, lambdamax, mm]= param_stabilize(B, N, lambda);
	    
% Inverse TFMA Filtering
      eprime= tfarma_inv(uM, 1, Bml1/param_get(Bml1, 0, 0));
      e     = tfarma_inv(uM, 1, Bml1);
      E(M, L+1)     = e'*e/N;
      Eprime(M, L+1)= eprime'*eprime/N;
   end;
end;


e= log(E);
eprime= log(Eprime);

rho= 1-log(12);

MDL= eprime + (log(N+1)+rho)  *(NA+NB+1/2)/N;
BIC= eprime + log(N)          *(NA+NB    )/N;
AIC= eprime + 2               *(NA+NB    )/N;
GIC= eprime + (2+rho)         *(NA+NB    )/N;
MIC= eprime + (1+rho+log(N)/2)*(NA+NB    )/N;

mdl= e + (log(N+1)+rho)  *(NA+NB+1/2)/N;
bic= e + log(N)          *(NA+NB    )/N;
aic= e + 2               *(NA+NB    )/N;
gic= e + (2+rho)         *(NA+NB    )/N;
mic= e + (1+rho+log(N)/2)*(NA+NB    )/N;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
for n= 8:8
   N= 2^n;
   re_im= 'i';
   mo_no= 'n';
   alpha= 1/2;
   beta = 1/2;
   lambda= .9800;
   MM= 20;

   for M= 1:1
      for L= 2:2
         MAR  =   M;
         LAR  =   L;
         MMA= MAR-1;
         LMA  = LAR;
         tfpm_file_gen;
%-------------------------

         Mmax= max([2*M 2]);
         Lmax= max([2*L 2]);
         Psi= ones(N);
         [Psi, Mask, v2]= tf_multiwin(N, 2*Mmax+2, 2*Lmax+2, 5, 0, 0);
         v2
         MDLhist= zeros(Mmax, Lmax+1);
         BIChist= zeros(Mmax, Lmax+1);
         AIChist= zeros(Mmax, Lmax+1);
         GIChist= zeros(Mmax, Lmax+1);
         MIChist= zeros(Mmax, Lmax+1);
         mdlhist= zeros(Mmax, Lmax+1);
         bichist= zeros(Mmax, Lmax+1);
         aichist= zeros(Mmax, Lmax+1);
         gichist= zeros(Mmax, Lmax+1);
         michist= zeros(Mmax, Lmax+1);
         for mm= 1:MM
            x= tfarma_gen(randn(N, 1), Aml, Bml/param_get(Bml, 0, 0), 1/2);
            [MDL, BIC, AIC, GIC, MIC, mdl, bic, aic, gic, mic]= tfarma_fit(x, Mmax, Lmax, Psi, lambda);
            [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
            [MBIC, LBIC] = find(BIC==min(min(BIC)));LBIC= LBIC-1;
            [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
            [MGIC, LGIC] = find(GIC==min(min(GIC)));LGIC= LGIC-1;
            [MMIC, LMIC] = find(MIC==min(min(MIC)));LMIC= LMIC-1;
            [Mmdl, Lmdl] = find(mdl==min(min(mdl)));Lmdl= Lmdl-1;
            [Mbic, Lbic] = find(bic==min(min(bic)));Lbic= Lbic-1;
            [Maic, Laic] = find(aic==min(min(aic)));Laic= Laic-1;
            [Mgic, Lgic] = find(gic==min(min(gic)));Lgic= Lgic-1;
            [Mmic, Lmic] = find(mic==min(min(mic)));Lmic= Lmic-1;
            [N MAR LAR mm MMDL, LMDL, MBIC, LBIC, MAIC, LAIC, MGIC, LGIC, MMIC, LMIC]
            [N MAR LAR mm Mmdl, Lmdl, Mbic, Lbic, Maic, Laic, Mgic, Lgic, Mmic, Lmic]
            MDLhist(MMDL, LMDL+1)= MDLhist(MMDL, LMDL+1) + 1;
            BIChist(MBIC, LBIC+1)= BIChist(MBIC, LBIC+1) + 1;
            AIChist(MAIC, LAIC+1)= AIChist(MAIC, LAIC+1) + 1;
            GIChist(MGIC, LGIC+1)= GIChist(MGIC, LGIC+1) + 1;
            MIChist(MMIC, LMIC+1)= MIChist(MMIC, LMIC+1) + 1;
            mdlhist(Mmdl, Lmdl+1)= mdlhist(Mmdl, Lmdl+1) + 1;
            bichist(Mbic, Lbic+1)= bichist(Mbic, Lbic+1) + 1;
            aichist(Maic, Laic+1)= aichist(Maic, Laic+1) + 1;
            gichist(Mgic, Lgic+1)= gichist(Mgic, Lgic+1) + 1;
            michist(Mmic, Lmic+1)= michist(Mmic, Lmic+1) + 1;
	    
   figure(1);subplot(2, 2, 1);imagesc(MDLhist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(1);subplot(2, 2, 2);imagesc(BIChist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(1);subplot(2, 2, 3);imagesc(AIChist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(1);subplot(2, 2, 4);imagesc(GIChist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(2);subplot(2, 2, 1);imagesc(mdlhist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(2);subplot(2, 2, 2);imagesc(bichist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(2);subplot(2, 2, 3);imagesc(aichist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
   figure(2);subplot(2, 2, 4);imagesc(gichist)
   set(gca, 'XTick', 1:Lmax+1)
   set(gca, 'XTickLabel', 0:Lmax)
            drawnow
         end;
         exp= 'ord';
         resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no)
         save(resultname, 'lambda', 'Mmax', 'Lmax', 'MM', 'MDLhist', 'BIChist', 'AIChist', 'GIChist', 'MIChist', 'mdlhist', 'bichist', 'aichist', 'gichist', 'michist')
      end;
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
re_im= 'i';
mo_no= 'n';
alpha= 1/2;
beta = 1/2;

M= 3
L= 2
MSEn= [];
meanMn= [];
meanLn= [];
for n= 6:9
   N= 2^n;
   exp= 'ord';
   resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no)
   load(resultname)

   MDLM= sum(MDLhist.').';MDLL= sum(MDLhist);
   BICM= sum(BIChist.').';BICL= sum(BIChist);
   AICM= sum(AIChist.').';AICL= sum(AIChist);
   GICM= sum(GIChist.').';GICL= sum(GIChist);
   MICM= sum(MIChist.').';MICL= sum(MIChist);
   mdlM= sum(mdlhist.').';mdlL= sum(mdlhist);
   bicM= sum(bichist.').';bicL= sum(bichist);
   aicM= sum(aichist.').';aicL= sum(aichist);
   gicM= sum(gichist.').';gicL= sum(gichist);
   micM= sum(michist.').';micL= sum(michist);

figure(n);clf;
subplot(4, 2, 1);bar(MDLM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 2);bar(MDLL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 3);bar(BICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 4);bar(BICL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 5);bar(AICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 6);bar(AICL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 7);bar(GICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 8);bar(GICL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

mean= [
    sum((1:Mmax)'.*MDLM)/MM;
    sum((1:Mmax)'.*BICM)/MM;
    sum((1:Mmax)'.*AICM)/MM;
    sum((1:Mmax)'.*GICM)/MM;
    sum((1:Mmax)'.*MICM)/MM;
    sum((1:Mmax)'.*mdlM)/MM;
    sum((1:Mmax)'.*bicM)/MM;
    sum((1:Mmax)'.*aicM)/MM;
    sum((1:Mmax)'.*gicM)/MM;
    sum((1:Mmax)'.*micM)/MM;
    ];
meanMn= [meanMn mean];

mean= [
    sum((0:Lmax).*MDLL)/MM;
    sum((0:Lmax).*BICL)/MM;
    sum((0:Lmax).*AICL)/MM;
    sum((0:Lmax).*GICL)/MM;
    sum((0:Lmax).*MICL)/MM;
    sum((0:Lmax).*mdlL)/MM;
    sum((0:Lmax).*bicL)/MM;
    sum((0:Lmax).*aicL)/MM;
    sum((0:Lmax).*gicL)/MM;
    sum((0:Lmax).*micL)/MM;
    ];
meanLn= [meanLn mean];

mse= [
       10*log10((sum(((1:Mmax)'-M).^2.*MDLM) + sum(((0:Lmax)-L).^2.*MDLL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*BICM) + sum(((0:Lmax)-L).^2.*BICL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*AICM) + sum(((0:Lmax)-L).^2.*AICL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*GICM) + sum(((0:Lmax)-L).^2.*GICL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*MICM) + sum(((0:Lmax)-L).^2.*MICL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*mdlM) + sum(((0:Lmax)-L).^2.*mdlL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*bicM) + sum(((0:Lmax)-L).^2.*bicL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*aicM) + sum(((0:Lmax)-L).^2.*aicL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*gicM) + sum(((0:Lmax)-L).^2.*gicL))/MM/(M^2+L^2));
       10*log10((sum(((1:Mmax)'-M).^2.*micM) + sum(((0:Lmax)-L).^2.*micL))/MM/(M^2+L^2));
       ];
   MSEn= [MSEn mse];
   drawnow
end

N= 256;
L= 2;
MSEm= [];
meanMm= [];
meanLm= [];
for M= 2:5
   exp= 'ord';
   resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no)
   load(resultname)

   MDLM= sum(MDLhist.').';MDLL= sum(MDLhist);
   BICM= sum(BIChist.').';BICL= sum(BIChist);
   AICM= sum(AIChist.').';AICL= sum(AIChist);
   GICM= sum(GIChist.').';GICL= sum(GIChist);
   MICM= sum(MIChist.').';MICL= sum(MIChist);
   mdlM= sum(mdlhist.').';mdlL= sum(mdlhist);
   bicM= sum(bichist.').';bicL= sum(bichist);
   aicM= sum(aichist.').';aicL= sum(aichist);
   gicM= sum(gichist.').';gicL= sum(gichist);
   micM= sum(michist.').';micL= sum(michist);

figure(M);clf;
subplot(4, 2, 1);bar(MDLM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 2);bar(MDLL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 3);bar(BICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 4);bar(BICL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 5);bar(AICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 6);bar(AICL, 1);axis([.5 Lmax+1.5 0 MM])
set(gca, 'XTick', 1:Lmax+1)
set(gca, 'XTickLabel', 0:Lmax)

subplot(4, 2, 7);bar(GICM, 1);axis([.5 Mmax+0.5 0 MM])
subplot(4, 2, 8);bar(GICL, 1);axis([.5 Lmax+1.5 0 MM])