vtfarma_weyl.m

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

function Lalpha= vtfarma_weyl(AML, BML, N, B, alpha, tol)
% function Lalpha= vtfarma_weyl(AML, BML, N, B, alpha, tol)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the Weyl matrix of the sytem. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
N= 256;
J= 0;
B= J+1;
offset= 200;
M= 4;
L= 3;
filename= sprintf('chpar%02d%04d%02d%02d-%04d', J, N, M, L, offset)
load(filename)
BML= B0L;
tol= 0.01;


clear;tfpm
alpha= 1/2;
N= 128;
M0= 3;
L0= 2;
J0= 1;
B= J0;

filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)
load(filename);
BML= B0L;
tol= 0.01;

%AML= eye(J0);
%BML= eye(J0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

Anm= ifft(fftshift(aml, 3), [], 3);
clear aml
Ank= N*fftshift(fft(fftshift(Anm, 4), [], 4), 4);
clear Anm 
Bnm= ifft(fftshift(bml, 3), [], 3);
clear bml
Bnk= N*fftshift(fft(fftshift(Bnm, 4), [], 4), 4);
clear Bnm


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      figure(1);mesh(abs(reshape(Ank(j, jp, :, :), N, N)))
      figure(2);mesh(abs(reshape(Bnk(j, jp, :, :), N, N)))
%      figure(1);tf_show(reshape(Ank(j, jp, :, :), N, N))
%      figure(2);tf_show(reshape(Bnk(j, jp, :, :), N, N))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Regularized inversion
for n= 1:N
   for k= 1:N
      [Unk, Snk, Vnk]= svd(Ank(:, :, n, k));
%      Snk
      Snk= Snk+tol*Snk(1, 1)*eye(J);
%      Lalpha(:, :, n, k)= pinv(Ank(:, :, n, k), 1e-2)*Bnk(:, :, n, k);
%      Lalpha(:, :, n, k)= inv(Ank(:, :, n, k))*Bnk(:, :, n, k);
%      Lalpha(:, :, n, k)= Unk*inv(Snk)*Vnk'*Bnk(:, :, n, k);
      Lalpha(:, :, n, k)= inv(Vnk)'*inv(Snk)*inv(Unk)*Bnk(:, :, n, k);
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 1:J
   for jp= 1:J
      figure(2);mesh(abs(reshape(Lalpha(j, jp, :, :), N, N)))
      pause
   end;
end;
Lalpha0= tfarma_weyl(reshape(AML, 2*LAR+1, MAR+1), reshape(BML, 2*LMA+1, MMA+1), N, alpha);
figure(1);cesh(Lalpha0)
figure(2);cesh(reshape(Lalpha, N, N))
figure(3);cesh(Lalpha0-reshape(Lalpha, N, N))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%