tfarma_fit2.m

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

function [E, Eprime]= tfarma_fit2(x, Mmax, Lmax, Psi, lambda)
% function [E, Eprime]= tfarma_fit2(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    = 256;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  =   2;
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);
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);

NA= (1:Mmax)'*(2*(0:Lmax)+1);
NB= NA - 1;

for L= 0:Lmax
   M= 1;
      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);
      [Aml1, lambdamax, mm, nrFDIR]= param_stabilize(AA(:, 1:M+1, M), N, lambda);
      Aml1= Aml1./param_get(Aml1, 0, 0);
      Aml1(1:L, 1)= zeros(L, 1);
      Aml1(L+2:end, 1)= zeros(L, 1);
      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, nrFDIR]= 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, 1)     = e'*e/N;
      Eprime(M, L+1, 1)= eprime'*eprime/N;
   
   for M= 2: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);
      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);
         Aml1= Aml1./param_get(Aml1, 0, 0);
         Aml1(1:L, 1)= zeros(L, 1);
         Aml1(L+2:end, 1)= zeros(L, 1);
%         figure(3);plot(abs(param_tdir(Aml1, N)))
         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, nrFDIR]= param_stabilize(B.*Mask(:, 2:end), 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, L-LL+1)     = e'*e/N;
         Eprime(M, L+1, L-LL+1)= eprime'*eprime/N;
%	 Aml1
%	 Bml1
%	 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= 5:5
      for L= 2:2
         MAR  =   M;
         LAR  =   L;
         MMA= MAR-1;
         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]= tfarma_fit2(x, Mmax, Lmax, Psi, lambda);
	    EE(:, :, :, mm)= E;
	    EEprime(:, :, :, mm)= Eprime;
            [N MAR LAR mm]
	 end;
         exp= 'ord2C';
         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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%