tfarma_llfu.m

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

function [l, v]= tfarma_llfu(theta_AB, Pyy, MAR, LAR, MMA, LMA)
% function [l, v]= tfarma_llfu(theta_AB, Pyy, MAR, LAR, MMA, LMA)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computation of the NEGATIVE log-likelihood function (underspread
% approximation) of a real-valued Gaussian TFARMA(MAR, LAR; MMA, LMA; N; 1/2)
% models with split-parameter vector theta_AB. 

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

Amlref= Aml;
Bmlref= Bml;
theta_AB= [param_stack_ml(param_split(Aml(:, 2:end))); param_stack_ml(param_split(Bml))];
%theta_AB= [param_stack_ml(param_split(Aml(:, 2:end)))];

Pyy= tfarma_wvsp(Aml, Bml, N, alpha);
figure(1);mesh(Pyy)

[Psi, Mask]= tf_multiwin(N, 4*(MAR+LAR), 4*LAR, 5, 1, 1);
%Psi= ones(N);

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
Ryy= corr_est(y, y, -1, 1/2);
Ayy= nm_to_ml(Ryy);
Pyy= real(nm_to_nk(ml_to_nm(Ayy.*conj(Psi))));
figure(2);mesh(Pyy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha= 1/2;
beta = 1/2;

% Dimensions:
N= size(Pyy, 1);

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

A= param_destack_ml(theta_A, MAR-1, LAR);
B= param_destack_ml(theta_B, MMA  , LMA);

A= [[zeros(LAR, 1); 1; zeros(LAR, 1)] A];

Aml= param_merge(A);
Bml= param_merge(B);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Amlref-Aml
Bmlref-Bml
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Bnm= param_expand(Bml, N);
l=  2*sum(log(abs(Bnm(:, 1)))) + real(sum(sum(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= real(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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE SCORE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargout>1)
   LA= param_weyl(Aml, N, alpha);
   LB= param_weyl(Bml, N, alpha);
   WA= abs(LA).^2;
   WB= abs(LB).^2;
   v= [];
% The TFAR part
   for m0= 1:MAR
      vr= [];
      vi= [];
      km0= ones(N, 1)*exp(2*pi*j/N*(-N/2:N/2-1)*m0);
      vr00= -2*sum(sum(Pyy./WB.*real(LA.*km0)))/N;
      for l0= 1:LAR
         coco= 2*cos(2*pi/N*l0*(0:N-1)).'*ones(1, N);
         sisi= 2*sin(2*pi/N*l0*(0:N-1)).'*ones(1, N);
         vrm0l0= -2*sum(sum(coco.*Pyy./WB.*real(LA.*km0)))/N;
	 vim0l0=  2*sum(sum(sisi.*Pyy./WB.*real(LA.*km0)))/N;
	 vr= [vrm0l0; vr];
	 vi= [vi; vim0l0];
      end;
      v= [v; -vr; -vr00; -vi];
   end;
% The TFMA part for m0= 0
   vr= [];
   vi= [];
   vr00= -2*sum(sum(1./real(Bnm(:, 1)*ones(1, N)) - WA.*Pyy./WB./WB.*real(LB)))/N;
   for l0= 1:LMA
      coco= 2*cos(2*pi/N*l0*(0:N-1)).'*ones(1, N);
      sisi= 2*sin(2*pi/N*l0*(0:N-1)).'*ones(1, N);
      vrm0l0= -2*sum(sum(coco.*(1./real(Bnm(:, 1)*ones(1, N)) - WA.*Pyy./WB./WB.*real(LB))))/N;
      vim0l0=  2*sum(sum(sisi.*(1./real(Bnm(:, 1)*ones(1, N)) - WA.*Pyy./WB./WB.*real(LB))))/N;
      vr= [vrm0l0; vr];
      vi= [vi; vim0l0];
   end;
   v= [v; -vr; -vr00; -vi];
% The TFMA part for m0> 0
   for m0= 1:MMA
      vr= [];
      vi= [];
      km0= ones(N, 1)*exp(2*pi*j/N*(-N/2:N/2-1)*m0);
      vr00= 2*sum(sum(WA.*Pyy./WB./WB.*real(LB.*km0)))/N;
      for l0= 1:LMA
         coco= 2*cos(2*pi/N*l0*(0:N-1)).'*ones(1, N);
         sisi= 2*sin(2*pi/N*l0*(0:N-1)).'*ones(1, N);
         vrm0l0=  2*sum(sum(WA.*Pyy./WB./WB.*coco.*real(LB.*km0)))/N;
	 vim0l0= -2*sum(sum(WA.*Pyy./WB./WB.*sisi.*real(LB.*km0)))/N;
	 vr= [vrm0l0; vr];
	 vi= [vi; vim0l0];
      end;
      v= [v; -vr; -vr00; -vi];
   end;
end;

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

Amlref= Aml;
Bmlref= Bml;
theta_AB= [param_stack_ml(param_split(Aml(:, 2:end))); param_stack_ml(param_split(Bml))];
%theta_AB= [param_stack_ml(param_split(Aml(:, 2:end)))];

Pyy= tfarma_wvsp(Aml, Bml, N, alpha);
figure(1);mesh(Pyy)

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
Ryy= corr_est(y, y, -1, 1/2);
Ayy= nm_to_ml(Ryy);
[Psi, Mask]= tf_multiwin(N, 2*MAR, 2*LAR, 5, 1, 1);
Psi= ones(N);
Pyy= real(nm_to_nk(ml_to_nm(Ayy.*conj(Psi))));
figure(2);mesh(Pyy)


delta= 0.0001;
for ii= 1:length(theta_AB)
   vv= zeros(size(theta_AB));
   vv(ii)= delta;
   [l1, v1]= tfarma_llfu(theta_AB   , Pyy, MAR, LAR, MMA, LMA);
   [l2, v2]= tfarma_llfu(theta_AB+vv, Pyy, MAR, LAR, MMA, LMA);
   [(l2-l1)/delta v1(ii)]
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%