vtfarma_inv.m

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

function E= vtfarma_inv(Y, AML, BML, B, beta)
% function E= vtfarma_inv(Y, AML, BML, B, beta)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Generates the output of the inverse TFARMA(MAR, LAR, MMA, LMA; N; beta)
% model. The TF model parameters a_ml resp. b_ml are given in the
% (MAR+1) x (2LAR+1) resp. (MMA+1) x (2LMA+1) matrices Aml
% resp. Bml, where N is the length of the input signal vector e. 

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

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

E0= randn(N, J);
%BML= eye(J);
Y= vtfarma_gen(E0, AML, BML, B, beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions
N= size(Y, 1);
[J, J, LAR, MAR]= size(BML);
LAR= (LAR-1)/2;
MAR= MAR-1;
[J, J, LMA, MMA]= size(AML);
LMA= (LMA-1)/2;
MMA= MMA-1;
Y= [Y; Y; Y];

% Generate c_m[n]
bml= zeros(J, J, N, MAR+1);
aml= zeros(J, J, N, MMA+1);
for j= 1:J
   for jp= 1:J
      if(abs(j-jp)<B)
         bml(j, jp, N/2+1-LAR:N/2+1+LAR, :)= BML(j, jp, :, :);
         aml(j, jp, N/2+1-LMA:N/2+1+LMA, :)= AML(j, jp, :, :);
      else
         bml(j, jp, N/2+1-LAR:N/2+1+LAR, :)= zeros(size(BML(j, jp, :, :)));
         aml(j, jp, N/2+1-LMA:N/2+1+LMA, :)= zeros(size(AML(j, jp, :, :)));
      end;
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      figure(1);tf_show(reshape(bml(j, jp, :, :), N, MAR+1))
      figure(2);tf_show(reshape(aml(j, jp, :, :), N, MMA+1))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Anm= N*ifft(fftshift(aml, 3), [], 3);
clear aml
Bnm= N*ifft(fftshift(bml, 3), [], 3);
clear bml
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      figure(1);tf_show(reshape(Bnm(j, jp, :, :), N, MAR+1))
      figure(2);tf_show(reshape(Anm(j, jp, :, :), N, MMA+1))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

E= [];
for n= 0:2*N-1
   if(n>=N)
      nn= n-N;
   else
      nn= n;
   end;
   bb= inv(Bnm(:, :, nn+1));
   ee= zeros(1, J);
   for m= 0:MMA
      ee= ee + (bb*Anm(:, :, nn+1, m+1)*(Y(N+nn+1-m, :).')).';
   end;
   for m= 1:MAR
%      if(n-m>=0)
         ee= ee - (bb*Bnm(:, :, nn+1, m+1)*(E(nn+1-m, :).')).';
%      end;
   end;
   E= [E; ee];
end;
E= E(N+1:end, :);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E-E0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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  =   2;
LAR  =   1;
MMA  =   0;
LMA  = LAR;
N    = 512;
J    =   2;
B    =   J;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta= 1/2;
filename= sprintf('vtfar%02d%02d%02d%04d', MAR, LAR, J, N)
load(filename);
BML= B0L;

E= randn(N, J);
HH= zeros(N*J);
for j= 1:J
   for jp= 1:J
      Aml= reshape(AML(j, jp, :, :), 2*LAR+1, MAR+1);
      B0l= reshape(B0L(j, jp, :)   , 2*LMA+1, 1);
      H= lag_to_ker(tfarma_impr(Aml, B0l, N, alpha), alpha);
      HH((j-1)*N+1:j*N, (jp-1)*N+1:jp*N)= H;
   end;
end;


Y0= reshape(HH*E(:), N, J);
Y= vtfarma_gen(E, AML, BML, B, beta);
norm(Y-Y0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%