tfarma_llfu1.m

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

function [l, v]= tfarma_llfu1(theta_ABreim, Pyy, MAR, LAR, MMA, LMA, re_im)
% function [l, v]= tfarma_llfuy(theta_ABreim, Pyy, MAR, LAR, MMA, LMA, re_im)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computation of the NEGATIVE underspread log-likelihood function of
% Gaussian TFARMA(MAR, LAR; MMA, LMA; N; 1/2) models with
% parameters theta_AB. Takes care of re_im! 

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

Amlref= Aml;
Bmlref= Bml;
theta_AB= [param_stack_ml(Aml(:, 2:end)); param_stack_ml(Bml)];
theta_ABreim= [real(param_stack_ml(Aml(:, 2:end))); ...
	       imag(param_stack_ml(Aml(:, 2:end))); ...
	       real(param_stack_ml(Bml)); ...
	       imag(param_stack_ml(Bml))];

%Psi= tf_multiwin(N, 20, 20, 5, 1, 1);
Psi= ones(N);

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
qyy= corr_est(y, y, -1, 1/2);
ayy= fft(qyy);
ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
Ayy= ayy.*conj(Psi);

Ayy= tfarma_ambi(Aml, Bml, N, 1/2);
Pyy= real(nm_to_nk(ml_to_nm(Ayy)));
mesh(abs(Ayy))

%e0= randn(N, 1);
%y= tfarma_gen(e0, Aml, Bml, beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions:
N= length(Pyy);

theta_AB= [theta_ABreim(1:MAR*(2*LAR+1)) + j*theta_ABreim(MAR*(2* LAR+1)+1:2*MAR*(2*LAR+1));
	   theta_ABreim(2*MAR*(2*LAR+1)+1:2*MAR*(2*LAR+1)+(MMA+1)*(2*LMA+1)) ...
	  + j*theta_ABreim(2*MAR*(2*LAR+1)+(MMA+1)*(2*LMA+1)+1:end)];

theta_A= theta_AB(1:MAR*(2*LAR+1));
theta_B= theta_AB(MAR*(2*LAR+1)+1:end);

Aml= [[zeros(LAR, 1); 1; zeros(LAR, 1)] param_destack_ml(theta_A, MAR-1, LAR)];
Bml= param_destack_ml(theta_B, MMA, LMA);

Bnm= param_expand(Bml, N);

l=  2*sum(log(abs(Bnm(:, 1)))) + real(sum(sum(conj(Pyy)./tfarma_wvsp(Aml, Bml, N, 1/2))))/N;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slow, exact reference:
Anm= param_expand(Aml, N);
Bnm= param_expand(Bml, N);
H= tvarma_impr(Anm, Bnm);
R= H*H';
lref= log(real(det(R))) + real(y'*inv(R)*y);
[l lref]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE SCORE
if(nargout==2)
   difAml= -nk_to_nm(nm_to_ml( Pyy.*param_weyl(Aml, N, 1/2)./abs(param_weyl(Bml, N, 1/2)).^2 ));
   scoAml= 2*conj(difAml(N/2+1-LAR:N/2+1+LAR, N/2+2:N/2+1+MAR));

   difBml= +nk_to_nm(nm_to_ml( Pyy.*param_weyl(Bml, N, 1/2).*abs( ...
       param_weyl(Aml, N, 1/2)./param_weyl(Bml, N, 1/2) ).^2./abs(param_weyl(Bml, N, 1/2)).^2 ));
   scoBml= 2*conj(difBml(N/2+1-LMA:N/2+1+LMA, N/2+1:N/2+1+MMA));

   difB0l= -fft(1./Bnm(:, 1));
   scoB0l= 2*conj([difB0l(end-LMA+1:end); difB0l(1:LMA+1)]);

   scoBml(:, 1)= scoBml(:, 1) + scoB0l;

   v= -1/2*[...
       real(param_stack_ml(scoAml)); ...
       imag(param_stack_ml(scoAml)); ...
       real(param_stack_ml(scoBml)); ...
       imag(param_stack_ml(scoBml)); ...
      ];
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%