tfar_fit2.m

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

function [E, Eprime]= tfar_fit2(x, Mmax, Lmax, Psi, lambda)
% function [E, Eprime]= tfar_fit2(x, Mmax, Lmax, Psi, lambda)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Searches TFAR(M, 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  =   0;
LMA  = LAR;
re_im= 'i';
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);
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);

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

E= 99999*ones(Mmax, Lmax+1, Lmax+1);
Eprime= 99999*ones(Mmax, Lmax+1, Lmax+1);

for L= 0:Lmax
   M= 1;
      [AA, BB]= tfar_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2-M+1:N/2+M+1), N);
      [Aml1, lambdamax, mm, nrFDIR]= param_stabilize(AA(:, 1:M+1, M), N, lambda);
      Bml1= BB(:, :, M);
      e= tfarma_inv(x, Aml1, Bml1);
      eprime= tfarma_inv(x, Aml1, Bml1/param_get(Bml1, 0, 0));
      E(M, L+1, 1)     = e'*e/N;
      Eprime(M, L+1, 1)= eprime'*eprime/N;
   
   for M= 2:Mmax
      [AA, BB]= tfar_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2-M+1:N/2+M+1), N);
      for LL= 0:L
% Inverse TFAR Filtering
         Mask= ones(2*L+1, M+1);
         Mask(1:L-LL, end)= zeros(L-LL, 1);
         Mask(L+2+LL:end, end)= zeros(L-LL, 1);
%         figure(1);plot(abs(param_tdir(AA(:, 1:M+1, M), N)))
%         figure(2);plot(abs(param_tdir(AA(:, 1:M+1, M).*Mask, N)))
         [Aml1, lambdamax, mm, nrFDIR]= param_stabilize(AA(:, 1:M+1, M).*Mask, N, lambda);
         Bml1= BB(:, :, M);
         e= tfarma_inv(x, Aml1, Bml1);
         eprime= tfarma_inv(x, Aml1, Bml1/param_get(Bml1, 0, 0));
         E(M, L+1, L-LL+1)     = e'*e/N;
	 [M, L, LL];
         Eprime(M, L+1, L-LL+1)= eprime'*eprime/N;
%	 drawnow
      end;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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= 100;

   for M= 2:2
      for L= 3:3
         MAR  =   M;
         LAR  =   L;
         MMA  =   0;
         LMA  = LAR;
         tfpm_file_gen;
%-------------------------
         Bml= Bml/param_get(Bml, 0, 0);
         Mmax= max([2*M 2]);
         Lmax= max([2*L 2]);
         [Psi, Mask, v2]= tf_multiwin(N, 2*Mmax+2, 2*Lmax+2, 5, 0, 0);
%         Psi= ones(N);
         EE= zeros(Mmax, Lmax+1, Lmax+1, MM);
         EEprime= zeros(Mmax, Lmax+1, Lmax+1, MM);
         for mm= 1:MM
            x= tfarma_gen(randn(N, 1), Aml, Bml/param_get(Bml, 0, 0), 1/2);
            [E, Eprime]= tfar_fit2(x, Mmax, Lmax, Psi, lambda);
	    EE(:, :, :, mm)= E;
	    EEprime(:, :, :, mm)= Eprime;
            [N M L mm]
	 end;
         exp= 'ord2A';
         resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no)
         save(resultname, 'lambda', 'Mmax', 'Lmax', 'MM', 'EE', 'EEprime')
      end;
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%