vtfarma_gen.m

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

function Y= vtfarma_gen(E, AML, BML, B, beta)
% function Y= vtfarma_gen(E, 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 a 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= 256;
M=  2;
L=  1;
F=  2;
B=  2;
beta= 1/2;

filename= sprintf('vtfar%02d%02d%02d%04d', M, L, F, N)
load(filename);
BML= B0L;
E0= randn(N, F);
E= E0;
BML= eye(F);

X= vtfarma_gen(randn(N, F), AML, BML)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% GENERATE FILES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N=  256;
MAR=  2;
LAR=  1;
MMA=  0;
LMA=LAR;
F=  2;
B=  2;
beta= 1/2;
AML= zeros(F, F, 2*LAR+1, MAR+1);
AML(:, :, LAR+1, 1)= eye(F);
B0L= zeros(F, F, 2*LAR+1, 1);
B0L(:, :, LAR+1)= randn(F) + j*randn(F);
B0L(:, :, LAR+1)= B0L(:, :, LAR+1)*B0L(:, :, LAR+1)';
%B0L(:, :, LAR+1)= diag(1:F);
for jj= 1:F
   for jp= 1:F
      [jj jp]
      aml= [[zeros(LAR, 1); 1; zeros(LAR, 1)] randn(2*LAR+1, MAR)+i*randn(2*LAR+1, MAR)];
%      aml= param_stabilize(param_hermite(aml), N, .8, 1)
      aml= param_stabilize(aml, N, .8, 1)
      AML(jj, jp, :, :)= reshape(aml, 1, 1, 2*LAR+1, MAR+1);
%      tf_show(tfarma_weyl(aml, 1, N, 1/2));drawnow
      mesh(abs(tfarma_weyl(aml, 1, N, 1/2)));drawnow
   end;
end;
for l= 1:LAR
   AML(:, :, l, 1)= zeros(F);
   AML(:, :, 2*LAR+2-l, 1)= zeros(F);
end;
AML(:, :, LAR+1, 1)= eye(F);

E0= randn(N, F);
X= vtfarma_gen(E0, AML, B0L, F, alpha);
plot(abs(X))
filename= sprintf('vtfar%02d%02d%02d%04d', MAR, LAR, F, N)
save(filename, 'AML', 'B0L');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% GENERATE FILES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if(nargin<=3)
   B= F;
end;

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

Y= [];
for n= 0:N-1
   yy= zeros(1, F);
   for m= 0:MMA
      yy= yy + (Bnm(:, :, n+1, m+1)*(E(N+n+1-m, :).')).';
   end;
   for m= 1:MAR
      if(n-m>=0)
         yy= yy - (Anm(:, :, n+1, m+1)*(Y(n+1-m, :).')).';
      end;
   end;
   Y= [Y; yy];
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y1= zeros(N, F);
for j= 1:F
   y= zeros(N, 1);
   for jp= 1:F
      if(abs(j-jp)<B)
         Aml= reshape(AML(j, jp, :, :), 2*LAR+1, MAR+1);
         Bml= reshape(BML(j, jp, :, :), 2*LMA+1, MMA+1);
         y= y + tfarma_gen(E0(:, jp), Aml, Bml, 1/2);
      end;
   end;
   Y1(:, j)= Y1(:, j) + y;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
F    =   2;
B    =   F;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta= 1/2;
filename= sprintf('vtfar%02d%02d%02d%04d', MAR, LAR, F, N)
load(filename);
BML= B0L;

E= randn(N, F);
HH= zeros(N*F);
for j= 1:F
   for jp= 1:F
      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, F);
Y= vtfarma_gen(E, AML, BML, B, beta);
norm(Y-Y0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AML(:,:,1,1)= [
     0     0
     0     0
];

AML(:,:,2,1)= [
     1     0
     0     1
];

AML(:,:,3,1)= [
     0     0
     0     0
];

AML(:,:,1,2)= [
  -0.0301 + 0.0801i   0.0422 + 0.1571i
  -0.0488 - 0.0455i   0.1999 + 0.0740i
];

AML(:,:,2,2)= [
   0.2794 - 0.0079i   0.1020 - 0.0698i
  -0.1316 - 0.1935i   0.3989 + 0.2268i
];

AML(:,:,3,2)= [
  -0.0735 - 0.2695i  -0.2218 + 0.1556i
   0.2284 - 0.0626i  -0.0291 + 0.1997i
];

AML(:,:,1,3)= [
  -0.1848 + 0.1263i   0.0352 + 0.2071i
  -0.0592 + 0.0581i   0.0468 + 0.1225i
];

AML(:,:,2,3)= [
   0.1046 + 0.0451i  -0.0107 + 0.0301i
   0.0745 + 0.0095i   0.0314 - 0.0086i
];

AML(:,:,3,3)= [
  -0.0042 + 0.0378i  -0.0433 + 0.1037i
  -0.0274 - 0.2137i  -0.0920 + 0.1374i
];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%