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)%   michael.jachan@tuwien.ac.at 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-AmlBmlref-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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%