tfarma_est_gkmu.m

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

function [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N)
% function [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% MONIC TFARMA SYSTEM!!!!!
% MC, LC >= ??

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

Bml= .2*Bml;

KC=  4;
MC= 2*MAR-1;
LC= 2*LAR;

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
Ayy= Aalpha(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1);

%y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
%qyy= corr_est(y, y, MC, 1/2);
%ayy= fft(qyy);
%ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
%Ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :);

% Estimate the long TFAR model
[CA, CB]= tfar_est_tfywu(Ayy, N);
Cml= CA(:, :, end);
D0l= CB(:, :, end);

%Cml= ones(2*LC+1, 1)*(0:MC) + j*(-LC:LC)'*ones(1, MC+1);D0l= -1;

P0= tfarma_wvsp(Aml, Bml, N, 1/2);
figure(1);tf_show(P0)
P1= tfarma_wvsp(Cml, D0l, N, 1/2);
figure(2);tf_show(P1)
figure(3);tf_show(abs(P0-P1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha= 1/2;

% Dimensions
[MC, LC]= param_dim(Cml);
LL   = 2*LMA+1;

% Fill the Toeplitz/block-Toeplitz generator rPM
rPM= [];
for mu= MMA-1:-1:-MMA+1
   rPM= [rPM; fliplr(Cml(LC-2*LMA+1:LC+2*LMA+1, MAR+mu+1).')];
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toeptoep(rPM)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RHS= [];
for m= MAR+1:MAR+MMA
   RHS= [RHS; -Cml(LC-LMA+1:LC+LMA+1, m+1)*param_get(D0l, 0, 0)];
end;

% Solve the system
BBstack= invert_waxkailath(rPM, RHS);

B= [[zeros(LMA, 1); param_get(D0l, 0, 0); zeros(LMA, 1)] reshape(BBstack(:, :, end), LL, MMA)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reshape to BB
BB= zeros(2*LMA+1, 1+MMA, MMA);
for m= 1:MMA
   BB(:, :, m)= [[zeros(LMA, 1); param_get(D0l, 0, 0); zeros(LMA, 1)] reshape(BBstack(1:LL*m, 1, m), LL, m) zeros(LL, MMA-m)];
end;
param_norm(BB(:, :, end), Bml, 1)
param_norm(B, Bml, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The TFAR part
B1= [ [zeros(LC-LMA, 1+MMA); B/param_get(D0l, 0, 0); zeros(LC-LMA, 1+MMA)] zeros(2*LC+1, MC-MMA) ];
A= zeros(2*LAR+1, 1+MAR);
A(LAR+1, 1)= 1;
for m= 1:MAR
   for l= -LAR:LAR
      aml= Cml(LC+1+l, 1+m);
      for mprime= 1:m
         for lprime= -LAR:LAR
            aml= aml + Cml(LC+1+l-lprime, 1+m-mprime)*...
                       B1(LC+1+lprime, 1+mprime);
         end;
      end;
      A(LAR+1+l, 1+m)= aml;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm(A, Aml, 0)
param_norm(B, Bml, 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
MAR  =   3;
LAR  =   2;
MMA  =   2;
LMA  =   2;
N    = 256;
re_im= 'i';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
mm   =   1;

KC=  2;
LC= KC*(LAR);
MC= KC*(MAR);

MM= 1000;
%Psi= tf_multiwin(N, MC, LC, 5, 1, 1);
Psi= ones(N);
psi= Psi(N/2+1-3*LC:N/2+1+3*LC, N/2+1-MC:N/2+1+MC);

AA= zeros(size(Aml));
BB= zeros(size(Bml));
for mm= 1:MM
   y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
   qyy= corr_est(y, y, MC, 1/2);
   ayy= fft(qyy);
   ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
   Ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :).*conj(psi);
% Estimate the long TFAR model
   [CA, CB]= tfar_est_tfywu(Ayy, N);
   Cml= CA(:, :, end);
   D0l= CB(:, :, end);
   [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N);
   [mm param_norm(A, Aml, 0) param_norm(B, Bml, 0)]
   AA= AA+A;
   BB= BB+B;
end;
AA/MM
Aml
BB/MM
Bml
[param_norm(AA/MM, Aml, 0) param_norm(BB/MM, Bml, 0)]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%