vtfar_parest.m

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

function [AML, B0L]= vtfar_parest(ambi, M, L, N)
% function [AML, B0L]= vtfar_parest(ambi, M, L, N)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% get rid of pospart!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N= 256;
M0=  3;
L0=  2;
J0=  2;
B0= J0;

MMAX= 2*M0;
LMAX= 2*L0;

filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)
load(filename);
AML0= AML;
B0L0= B0L;

%E= randn(N, J0);X= vtfarma_gen(E, AML0, B0L0, B0, 1/2);[AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);

M= M0;
L= L0;
AMBI= vtfarma_ambi(AML0, B0L0, N, alpha);ambi= AMBI(:, :, N/2+1-3*LMAX:N/2+1+3*LMAX, N/2+1-MMAX:N/2+1+MMAX);
for j= 1:J0
   for jp= 1:J0
      [j jp]
      figure(2);tf_show(reshape(AMBI(j, jp, :, :), N, N))
      figure(4);mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))
      pause
   end;
end;



AA= zeros(J0);
m= 0;
l= 0;
for mp= 1:M
   for lp= -L:L
      AA= AA + AML0(:, :, L+1+lp, 1+mp)*ambi(:, :, 3*LMAX+1+l-lp, MMAX+1+m-mp);
   end;
end;
AA/N
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[J, J, LMAX, MMAX]= size(ambi);
LMAX= (LMAX-1)/6;
MMAX= (MMAX-1)/2;
alpha= 1/2;

% Fill the TBTB Matrix
A= zeros(J*M*(2*L+1));
for m= 1:M
   for mp= 1:M
      for l= -L:L
         for lp= -L:L
            A((m -1)*(2*L+1)*J+(L+l )*J+1:(m -1)*(2*L+1)*J+(L+l +1)*J, ...
              (mp-1)*(2*L+1)*J+(L+lp)*J+1:(mp-1)*(2*L+1)*J+(L+lp+1)*J)...
                = ambi(:, :, 3*LMAX+1+lp-l, MMAX+1+mp-m);
         end;
      end;
   end;
end;

% Fill the RHS vector
a= zeros(J*M*(2*L+1), J);
for m= 1:M
   for l= -L:L
      a((m-1)*(2*L+1)*J+(L+l)*J+1:(m-1)*(2*L+1)*J+(L+l+1)*J, :)= ambi(:, :, 3*LMAX+1-l, MMAX+1-m);
   end;
end;



% The fast inversion;-)
theta= -conj(inv(A)*a);



% Remember the parameters
AML= zeros(J, J, 2*L+1, M+1);
for j= 1:J
   for jp= 1:J
      Aml= [[zeros(L, 1); 1; zeros(L, 1)] zeros(2*L+1, M)];
      for m= 1:M
         for l= -L:L
            Aml(L+1+l, m+1)= theta((m-1)*(2*L+1)*J+(L+l)*J+jp, j);
         end;
      end;
      AML(j, jp, :, :)= Aml;
   end;
end;
for l= 1:L
   AML(:, :, l, 1)= zeros(J);
   AML(:, :, 2*L+2-l, 1)= zeros(J);
end;
AML(:, :, L+1, 1)= eye(J);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sum(sum(sum(sum(abs(AML-AML0)))))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute b_{0, l}
BBH= zeros(J, J, N);
for l= 0:2*L
   for mp= 0:M
      for lp= -L:L
         BBH(:, :, l+1)= BBH(:, :, l+1) + AML(:, :, L+lp+1, mp+1)*ambi(:, :, 3*LMAX+1+l-lp, MMAX+1-mp);
      end;
   end;
end;
for l= -2*L:-1
   for mp= 0:M
      for lp= -L:L
         BBH(:, :, N+1+l)= BBH(:, :, N+1+l) + AML(:, :, L+lp+1, mp+1)*ambi(:, :, 3*LMAX+1+l-lp, MMAX+1-mp);
      end;
   end;
end;
bbH= ifft(BBH, N, 3);
for n= 0:N-1
   [V, D]= eig(bbH(:, :, n+1));
   bbH(:, :, n+1)= V*pospart(sqrt(D))*V';
end;
BBH= fft(bbH, N, 3);

B0L= zeros(J, J, 2*L+1);
for l= 0:L
   B0L(:, :, L+l+1)= BBH(:, :, l+1);
end;
for l= -L:-1
   B0L(:, :, L+l+1)= BBH(:, :, N+1+l);
end;
B0L= B0L/N;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sum(sum(sum(sum(abs(B0L-B0L0)))))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N= 256;
M= 3;
L= 2;
J= 2;
B= J;

MMAX= 2*M;
LMAX= 2*L;

filename= sprintf('vtfar%02d%02d%02d%04d', M, L, J, N)
load(filename);
AML0= AML;
B0L0= B0L;

AAMMLL= zeros(size(AML0));
BB00LL= zeros(size(B0L0));
AAMMBBII= zeros(J, J, N, N);
MM= 100;
for mm= 1:MM
   mm
   E= randn(N, J);
   X= vtfarma_gen(E, AML0, B0L0, B, 1/2);
   [AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);
   [AML, B0L]= vtfar_parest(ambi, M, L, N);
   AAMMBBII= AAMMBBII + AMBI;
   AAMMLL= AAMMLL + AML;
   BB00LL= BB00LL + B0L;
end;
AAMMLL= AAMMLL/MM;
BB00LL= BB00LL/MM;
AAMMBBII= AAMMBBII/MM;
[norm(AML0(:)-AAMMLL(:))/norm(AML0(:)) norm(B0L0(:)-BB00LL(:))/norm(B0L0(:))]

AMBI= vtfarma_ambi(AML0, B0L0, N, alpha);

for j= 1:J
   for jp= 1:J
      [j jp]
      norm(reshape(AAMMBBII(j, jp, :, :), N, N)-reshape(AMBI(j, jp, :, :), N, N))/norm(reshape(AMBI(j, jp, :, :), N, N))
      figure(1);tf_show(reshape(AAMMBBII(j, jp, :, :), N, N))
      figure(2);tf_show(reshape(AMBI(j, jp, :, :), N, N))
      figure(3);mesh(abs(reshape(AAMMBBII(j, jp, :, :), N, N)))
      figure(4);mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))
      pause
   end;
end;

ambi= AMBI(:, :, N/2+1-3*LMAX:N/2+1+3*LMAX, N/2+1-MMAX:N/2+1+MMAX);

[AML, B0L]= vtfar_parest(ambi, M, L, N);
[norm(AML0(:)-AML(:)) norm(B0L0(:)-B0L(:))]

[real(AML) AML0]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%