vtfar_parest2.m

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

function [AML, B0L]= vtfar_parest2(ambi, M, L, B, N)
% function [AML, B0L]= vtfar_parest2(ambi, M, L, B, N)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% SPECIAL VERSION for SPAWC05 paper!!!!!!!!!!!!!

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

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

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

% Full width
%AMBI= vtfarma_ambi(AML0, B0L0, N, J, alpha);

% or banded
AMBI= vtfarma_ambi(AML0, B0L0, N, B, alpha);

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

m= 1;
mp= 1;
l= -L;
lp= -L;
i= 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      [j jp]
      figure(2);clf;tf_show(reshape(AMBI(j, jp, :, :), N, N))
      figure(4);clf;mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

% Fill the TBTB Matrix
F= zeros(M*(2*L+1)*NB);
for m= 1:M
 for mp= 1:M
  for l= -L:L
   for lp= -L:L
    Fml= zeros(NB);
    j= 1;
    for i= 0:J
     taumin= max(i-B+1, 0);
     taumax= min(i+B-1, J);
     [taumin taumax];
     Fmli= ambi(taumin+1:taumax+1, taumin+1:taumax+1, 3*LMAX+1+l-lp, MMAX+1+m-mp).';
%     [size(Fmli) j j+taumax-taumin taumax-taumin+1]
     Fml(j:j + (taumax-taumin), j:j + (taumax-taumin))= Fmli;
     j= j + (taumax-taumin)+1;
    end;
    iidx= (m -1)*(2*L+1)*NB + (L-l )*NB+1:(m -1)*(2*L+1)*NB + (L-l +1)*NB;
    jidx= (mp-1)*(2*L+1)*NB + (L-lp)*NB+1:(mp-1)*(2*L+1)*NB + (L-lp+1)*NB;
    F(iidx, jidx)= Fml;
   end;
  end;
 end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[size(F) rank(F) cond(F)]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fill the RHS vector
f= zeros(M*(2*L+1)*NB, 1);
for m= 1:M
  for l= -L:L
    fml= zeros(NB, 1);
    j= 1;
    for i= 0:J
     taumin= max(i-B+1, 0);
     taumax= min(i+B-1, J);
     [taumin taumax];
     fmli= ambi(i+1, taumin+1:taumax+1, 3*LMAX+1+l, MMAX+1+m).';
%     [size(Fmli) j j+taumax-taumin taumax-taumin+1]
     fml(j:j + (taumax-taumin))= fmli;
     j= j + (taumax-taumin)+1;
    end;
    iidx= (m -1)*(2*L+1)*NB + (L-l )*NB+1:(m -1)*(2*L+1)*NB + (L-l +1)*NB;
    f(iidx)= fml;
  end;
end;



% The fast inversion;-)
theta= -inv(F)*f;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta 
AML0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Reshape
AML= zeros(J+1, J+1, 2*L+1, M+1);
j= 1;
for m= 1:M
 for l= -L:L
  for i= 0:J
%     [m l tau taup]
     taumin= max(i-B+1, 0);
     taumax= min(i+B-1, J);
     [taumin taumax];
     AML(i+1, taumin+1:taumax+1, L+1-l, m+1)= theta(j:j + (taumax-taumin));
     j= j + (taumax-taumin)+1;
  end;
 end;
end;


for l= 1:L
   AML(:, :, l, 1)= zeros(J+1);
   AML(:, :, 2*L+2-l, 1)= zeros(J+1);
end;
AML(:, :, L+1, 1)= eye(J+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[AML AML0]
AML-AML0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AMBI0= vtfarma_ambi(AML0, eye(J0), N, J0, alpha);
AMBIE= vtfarma_ambi(AML , eye(J0), N, J0, alpha);
for j= 1:J0
   for jp= 1:J0
      [j jp]
      figure(1);clf;mesh(abs(reshape(AMBI0(j, jp, :, :), N, N)))
      figure(2);clf;mesh(abs(reshape(AMBIE(j, jp, :, :), N, N)))
      figure(3);clf;mesh(abs(reshape(AMBI0(j, jp, :, :)-AMBIE(j, jp, :, :), N, N)))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% Compute b_{0, l}
BBH= zeros(J+1, J+1, 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+1, J+1, 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= 2;
L= 1;
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_parest2(ambi, M, L, B, 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(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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%