tfar_est_mlu.m

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

function [A, B, exitflag, output]= tfar_est_mlu(Ayy, Aml0, Bml0, opts)
% function [A, B, exitflag, output]= tfar_est_mlu(Ayy, Aml0, Bml0, opts)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Underspread ML estimator for real-valued TFAR(M, L; N; 1/2)
% models. Log-likelihood in tfar_llfu.m. 

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

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, :)];
Psi= tf_multiwin(N, MAR, LAR, 5, 1, 1);

Ayy= ayy.*conj(Psi);
ayy= Ayy(N/2+1-3*LAR:N/2+1+3*LAR, N/2+1-MAR:N/2+1+MAR);
[AA, BB]= tfar_est_tfywu(ayy, N);
Aml0= AA(:, :, end);
Bml0= BB(:, :, end);

%opts= optimset('MaxIter', N, 'MaxFunEvals', N, 'GradObj', 'on', 'DerivativeCheck', 'off', 'Display', 'iter', 'Diagnostics', 'on');
opts= optimset('GradObj', 'on', 'Display', 'iter', 'Diagnostics', 'on');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions:
[MAR, LAR]= param_dim(Aml0);

theta_AB= [param_stack_ml(param_split(Aml0(:, 2:end))); param_stack_ml(param_split(Bml0))];
%theta_AB= [param_stack_ml(param_split(Aml0(:, 2:end)))];

% RANDOM INIT. TRY OUT!!!!!!!!!!!!!!
%theta_AB= randn(size(theta_AB))

Pyy= real(nm_to_nk(ml_to_nm(Ayy)));

% Call fminunc:
[X, lX, exitflag, output]= fminunc('tfar_llfu', theta_AB, opts, Pyy, MAR, LAR);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call fmincon:NIET JUIST!!!
[X, lX, exitflag, output]= fmincon('tfar_llfu', theta_AB, [], [], [], [], -3*ones(size(theta_AB)), 3*ones(size(theta_AB)), [], opts, Pyy, MAR, LAR);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AB= param_destack_ml(X, MAR, LAR);
%AB= param_destack_ml(X, MAR-1, LAR);
A= [[zeros(LAR, 1); 1; zeros(LAR, 1)] AB(:, 1:end-1)];
%A= [[zeros(LAR, 1); 1; zeros(LAR, 1)] AB(:, 1:end)];
B= AB(:, end);
%B= [zeros(LAR, 1); 1; zeros(LAR, 1)];

A= param_merge(A);
B= param_merge(B);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm([Aml0 Bml0], [Aml Bml], 1)
param_norm([A    B   ], [Aml Bml], 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    =  64;
MAR  =   2;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;
%-------------
beta = 1/2;
mm= 1;
MM= 10;
opts= optimset('GradObj', 'on', 'Display', 'iter', 'Diagnostics', 'on');

AMLyw= zeros(2*LAR+1, MAR+2, MM);
AMLml= zeros(2*LAR+1, MAR+2, MM);
%AMLyw= zeros(2*LAR+1, MAR+1, MM);
%AMLml= zeros(2*LAR+1, MAR+1, MM);

Iter= [];
NFun= [];
EFLG= [];
Psi= tf_multiwin(N, MAR, LAR, 5, 1, 1);
for mm= 1:MM
   mm
   y= real(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= ayy(N/2+1-3*LAR:N/2+1+3*LAR, N/2+1-MAR:N/2+1+MAR);
   [AA, BB]= tfar_est_tfywu(ayy, N);
   Aml0= AA(:, :, end);
   Bml0= BB(:, :, end);
   
   [A, B, exitflag, output]= tfar_est_mlu(Ayy, Aml0, Bml0, opts);
   Iter= [Iter; output.iterations];
   NFun= [NFun; output.funcCount];
   EFLG= [EFLG; exitflag];
   figure(1);clf;plot(Iter);drawnow
   figure(2);clf;plot(NFun);drawnow
   AMLyw(:, :, mm)= [Aml0 Bml0];
   AMLml(:, :, mm)= [A B];
%   AMLyw(:, :, mm)= Aml0;
%   AMLml(:, :, mm)= A;
end;
[mean(Iter) mean(NFun)]
[Myw, Vyw, Byw]= param_mse(AMLyw, [Aml Bml])
[Mml, Vml, Bml]= param_mse(AMLml, [Aml Bml])
%[Myw, Vyw, Byw]= param_mse(AMLyw, Aml)
%[Mml, Vml, Bml]= param_mse(AMLml, Aml)
[norm(Myw) norm(Mml)]
mean(AMLyw, 3)
mean(AMLml, 3)
sum(EFLG)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%