simu_win.m

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

function simu_win(M, L, N, re_im, MM)
% function simu_win(M, L, N, re_im, MM)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates <re_im>-valued TF(M, L, N) models to find optimum
% window lengthes. <MM> realizations are done to estimate the MSE. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
M    =   2;
L    =   2;
re_im= 'i';
MM   = 2;
xx= 3;
yy= 3;
mm= 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mo_no= 'n';

% Generate TFMA/TFAR model:
MAR= M;
LAR= L;
MMA= M;
LMA= L;
N0= N;
N= 128;
tfpm_file_gen;
%-------------
N= N0;
BmlB= Bml./param_get(Bml, 0, 0);
AmlA= Aml;
BmlA= BmlB(:, 1);
parA= [AmlA(:, 2:end) BmlA];
parB= BmlB;

% Generate TFARMA model:
MAR= M;
LAR= L;
MMA= M-1;
LMA= L;
N0= N;
N= 64;
tfpm_file_gen;
%-------------
N= N0;
AmlC= Aml;
BmlC= Bml./param_get(Bml, 0, 0);
parC= [AmlC(:, 2:end) BmlC];

TDIRB= param_tdir(param_expand(BmlB, N));
TDIRA= param_tdir(param_expand(AmlA, N));
TDIRC= [ param_tdir(param_expand(AmlC, N)) param_tdir(param_expand(BmlC, N)) ];
lambda0B= max(max(abs(TDIRB)))
lambda0A= max(max(abs(TDIRA)))
lambda0C= max(max(abs(TDIRC)))

exp= 'tfwin';
resultname= sprintf('data/%04d/%s%d%d%s.mat', N, exp, M, L, re_im)

% Window lengthes (wing)
XX= 0:4*M;

% Window lengthes (winh)
YY= 0:4*L;

xxmax= length(XX);
yymax= length(YY);

MSE= zeros(xxmax, yymax, 8);
HYP= zeros(xxmax, yymax, 8, MM);
COM= zeros(xxmax, yymax);


xx=1;yy=1;
      [Psi, Mask, v2]= tf_multiwin(N, XX(xx), YY(yy), 0, 2, 1);
      psi= Psi(N/2-3*L+1:N/2+3*L+1, N/2-M+1:N/2+M+1);
      COM(xx, yy)= length(v2);
      
      BML_B1= zeros(2*L+1, M+1, MM);
      BML_B2= zeros(2*L+1, M+1, MM);
      
      AML_A1= zeros(2*L+1, M+1, MM);
      AML_A2= zeros(2*L+1, M+1, MM);
      
      AML_C0 = zeros(2*L+1, M  , MM);
      BML_C1a= zeros(2*L+1, M  , MM);
      BML_C1b= zeros(2*L+1, M  , MM);
      BML_C2 = zeros(2*L+1, M  , MM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         [N M L xx length(XX) yy length(YY) mm]
         e= randn(N, 1);
         yB= tfarma_gen(e,    1, BmlB, 1/2);
         yA= tfarma_gen(e, AmlA, 1   , 1/2);
         yC= tfarma_gen(e, AmlC, BmlB, 1/2);

         aB= fft(corr_est(yB, yB, -1, 1/2));
         aA= fft(corr_est(yA, yA, -1, 1/2));
         aC= fft(corr_est(yC, yC, -1, 1/2));

         aB= [aB(N/2+1:N, :); aB(1:N/2, :)].*conj(Psi);
         aA= [aA(N/2+1:N, :); aA(1:N/2, :)].*conj(Psi);
         aC= [aC(N/2+1:N, :); aC(1:N/2, :)].*conj(Psi);

	 B1= tfma_est_cepsb(aB, M, L);
	 BB= tfma_est_lin(yB, aB, Psi, M, L, 2*M, L);B2= BB(:, :, end);

         BML_B1(:, :, mm)= B1;
         BML_B2(:, :, mm)= B2;
         TDIRB1= param_tdir(param_expand(B1             , N));
         TDIRB2= param_tdir(param_expand(B2             , N));
	 
	 [A, B, REG, PMIN]= tfar_est_cepsb(aA, M, L);A1= [A(:, 2:end) B];
	 [AA, BB]= tfar_est_tfywu(aA(N/2+1-3*L:N/2+1+3*L, N/2+1-M:N/2+1+M), N);A2= [AA(:, 2:end, end) BB(:, :, end)];

         AML_A1(:, :, mm)= A1;
         AML_A2(:, :, mm)= A2;
         TDIRA1= param_tdir(param_expand(A1(:, 1:end-1) , N));
         TDIRA2= param_tdir(param_expand(A2(:, 1:end-1) , N));
	 
	 A0= tfarma_est_tfywu(aC(N/2-3*L+1:N/2+3*L+1, N/2+M-1-M+1:N/2+M-1+M+1), M-1, N);
	 A0= A0(:, :, end);
	 B1a= tfarma_est_cepsb(aC, A0, M-1, L);
	 B1b= tfarma_est_ceps1(yC, Psi, A0, M-1, L);
	 BB2= tfarma_est_lin(yC, Psi, A0, M-1, L, 2*M, L);
	 B2= BB2(:, :, end);
	 
         AML_C0(:, :, mm) = A0(:, 2:end);
         BML_C1a(:, :, mm)= B1a;
         BML_C1b(:, :, mm)= B1b;
         BML_C2(:, :, mm) = B2;
         TDIRA0 = param_tdir(param_expand(A0 , N));
         TDIRC1a= param_tdir(param_expand(B1a, N));
         TDIRC1b= param_tdir(param_expand(B1b, N));
         TDIRC2 = param_tdir(param_expand(B2 , N));
	 
         HYP(xx, yy, 1, mm)= max(max(abs(TDIRB1)));
	 HYP(xx, yy, 2, mm)= max(max(abs(TDIRB2)));
	 HYP(xx, yy, 3, mm)= max(max(abs(TDIRA1)));
	 HYP(xx, yy, 4, mm)= max(max(abs(TDIRA2)));
         HYP(xx, yy, 5, mm)= max(max(abs(TDIRA0)));
	 HYP(xx, yy, 6, mm)= max(max(abs(TDIRC1a)));
	 HYP(xx, yy, 7, mm)= max(max(abs(TDIRC1b)));
	 HYP(xx, yy, 8, mm)= max(max(abs(TDIRC2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end;%for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      MSE(xx, yy, 1)= sum(sum(param_mse(BML_B1, parB)))/sum(sum(abs(parB).^2));
      MSE(xx, yy, 2)= sum(sum(param_mse(BML_B2, parB)))/sum(sum(abs(parB).^2));
      
      MSE(xx, yy, 3)= sum(sum(param_mse(AML_A1, parA)))/sum(sum(abs(parA).^2));
      MSE(xx, yy, 4)= sum(sum(param_mse(AML_A2, parA)))/sum(sum(abs(parA).^2));

      MSE(xx, yy, 5)= sum(sum(param_mse(AML_C0, parC(:, 2:M+1))))/sum(sum(abs(parC(:, 2:M+1)).^2));
      MSE(xx, yy, 6)= sum(sum(param_mse(BML_C1a, parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));
      MSE(xx, yy, 7)= sum(sum(param_mse(BML_C1b, parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));
      MSE(xx, yy, 8)= sum(sum(param_mse(BML_C2 , parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for xx= 2:xxmax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   for yy= 2:yymax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [Psi, Mask, v2]= tf_multiwin(N, XX(xx), YY(yy), 0, 2, 1);
      psi= Psi(N/2-3*L+1:N/2+3*L+1, N/2-M+1:N/2+M+1);
      COM(xx, yy)= length(v2);
      
      BML_B1= zeros(2*L+1, M+1, MM);
      BML_B2= zeros(2*L+1, M+1, MM);
      
      AML_A1= zeros(2*L+1, M+1, MM);
      AML_A2= zeros(2*L+1, M+1, MM);
      
      AML_C0 = zeros(2*L+1, M  , MM);
      BML_C1a= zeros(2*L+1, M  , MM);
      BML_C1b= zeros(2*L+1, M  , MM);
      BML_C2 = zeros(2*L+1, M  , MM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         [N M L xx length(XX) yy length(YY) mm]
         e= randn(N, 1);
         yB= tfarma_gen(e,    1, BmlB, 1/2);
         yA= tfarma_gen(e, AmlA, 1   , 1/2);
         yC= tfarma_gen(e, AmlC, BmlB, 1/2);

         aB= fft(corr_est(yB, yB, -1, 1/2));
         aA= fft(corr_est(yA, yA, -1, 1/2));
         aC= fft(corr_est(yC, yC, -1, 1/2));

         aB= [aB(N/2+1:N, :); aB(1:N/2, :)].*conj(Psi);
         aA= [aA(N/2+1:N, :); aA(1:N/2, :)].*conj(Psi);
         aC= [aC(N/2+1:N, :); aC(1:N/2, :)].*conj(Psi);

	 B1= tfma_est_cepsb(aB, M, L);
	 BB= tfma_est_lin(yB, aB, Psi, M, L, 2*M, L);B2= BB(:, :, end);

         BML_B1(:, :, mm)= B1;
         BML_B2(:, :, mm)= B2;
         TDIRB1= param_tdir(param_expand(B1             , N));
         TDIRB2= param_tdir(param_expand(B2             , N));
	 
	 [A, B, REG, PMIN]= tfar_est_cepsb(aA, M, L);A1= [A(:, 2:end) B];
	 [AA, BB]= tfar_est_tfywu(aA(N/2+1-3*L:N/2+1+3*L, N/2+1-M:N/2+1+M), N);A2= [AA(:, 2:end, end) BB(:, :, end)];

         AML_A1(:, :, mm)= A1;
         AML_A2(:, :, mm)= A2;
         TDIRA1= param_tdir(param_expand(A1(:, 1:end-1) , N));
         TDIRA2= param_tdir(param_expand(A2(:, 1:end-1) , N));
	 
	 A0= tfarma_est_tfywu(aC(N/2-3*L+1:N/2+3*L+1, N/2+M-1-M+1:N/2+M-1+M+1), M-1, N);
	 A0= A0(:, :, end);
	 B1a= tfarma_est_cepsb(aC, A0, M-1, L);
	 B1b= tfarma_est_ceps1(yC, Psi, A0, M-1, L);
	 BB2= tfarma_est_lin(yC, Psi, A0, M-1, L, 2*M, L);
	 B2= BB2(:, :, end);
	 
         AML_C0(:, :, mm) = A0(:, 2:end);
         BML_C1a(:, :, mm)= B1a;
         BML_C1b(:, :, mm)= B1b;
         BML_C2(:, :, mm) = B2;
         TDIRA0 = param_tdir(param_expand(A0 , N));
         TDIRC1a= param_tdir(param_expand(B1a, N));
         TDIRC1b= param_tdir(param_expand(B1b, N));
         TDIRC2 = param_tdir(param_expand(B2 , N));
	 
         HYP(xx, yy, 1, mm)= max(max(abs(TDIRB1)));
	 HYP(xx, yy, 2, mm)= max(max(abs(TDIRB2)));
	 HYP(xx, yy, 3, mm)= max(max(abs(TDIRA1)));
	 HYP(xx, yy, 4, mm)= max(max(abs(TDIRA2)));
         HYP(xx, yy, 5, mm)= max(max(abs(TDIRA0)));
	 HYP(xx, yy, 6, mm)= max(max(abs(TDIRC1a)));
	 HYP(xx, yy, 7, mm)= max(max(abs(TDIRC1b)));
	 HYP(xx, yy, 8, mm)= max(max(abs(TDIRC2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end;%for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      MSE(xx, yy, 1)= sum(sum(param_mse(BML_B1, parB)))/sum(sum(abs(parB).^2));
      MSE(xx, yy, 2)= sum(sum(param_mse(BML_B2, parB)))/sum(sum(abs(parB).^2));
      
      MSE(xx, yy, 3)= sum(sum(param_mse(AML_A1, parA)))/sum(sum(abs(parA).^2));
      MSE(xx, yy, 4)= sum(sum(param_mse(AML_A2, parA)))/sum(sum(abs(parA).^2));

      MSE(xx, yy, 5)= sum(sum(param_mse(AML_C0, parC(:, 2:M+1))))/sum(sum(abs(parC(:, 2:M+1)).^2));
      MSE(xx, yy, 6)= sum(sum(param_mse(BML_C1a, parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));
      MSE(xx, yy, 7)= sum(sum(param_mse(BML_C1b, parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));
      MSE(xx, yy, 8)= sum(sum(param_mse(BML_C2 , parC(:, M+1:end))))/sum(sum(abs(parC(:, M+1:end)).^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for yy= 1:winhmax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;%for xx= 1:wingmax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(5);imagesc(MSE(:, :));colorbar 
%figure(6);mesh(MSE(:, :));colorbar 

save(resultname, 'MSE', 'XX', 'YY', 'MM', 'lambda0B', 'lambda0B', 'lambda0C', 'HYP', 'COM')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
for n= 6:9
   simu_win(3, 2, 2^n, 'i', 100);
end;

clear;tfpm;
for m= 2:2
   simu_win(m, 2, 256, 'i', 100);
end;

clear;tfpm;
for m= 4:6
   simu_win(m, 2, 256, 'i', 100);
end;

clear;tfpm;
for l= 0:1
   simu_win(3, l, 256, 'i', 100);
end;

clear;tfpm;
for l= 3:4
   simu_win(3, l, 256, 'i', 100);
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;
exp= 'tfwin';
re_im= 'i';

smooth= [.025 .025 .025;
	 .025 .800 .025; 
	 .025 .025 .025];

smooth= [0 0 0;
	 0 1 0; 
	 0 0 0];


%NULL SEARCH!!

ESTI= 6;
mar= 3;
lar= 2;
for n= 6:9
   resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im)
   load(resultname)

   HYP0= sum(HYP(:, :, ESTI, :)>lambda0B, 4);
   HYP1= sum(HYP(:, :, ESTI, :)>       1, 4);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE(:, :, ESTI) , smooth);
   hyp0= hyp0(2:end-1, 2:end-1);
   hyp1= hyp1(2:end-1, 2:end-1);
   mse = mse (2:end-1, 2:end-1);
   
%   mse= MSE(:, :, ESTI);
%   hyp0= HYP0;
%   hyp1= HYP1;
   
   [wing, winh]= find(mse==min(min(mse(2:end, 2:end))));
   [2^n mar XX(wing) lar YY(winh) MSE(wing, winh) MSE(1, 1) HYP0(wing, winh) HYP0(1, 1) HYP1(wing, winh) HYP1(1, 1)]
%   [mar lar XX(wing) YY(winh) XX(end) YY(end) min(min(MSE))]
   MLopt(mar, lar, 1)= XX(wing);
   MLopt(mar, lar, 2)= YY(winh);
   figure(n);clf;
   subplot(3, 2, 1);mesh(mse);
   subplot(3, 2, 2);imagesc(mse);colorbar
   subplot(3, 2, 3);mesh(hyp0)
   subplot(3, 2, 4);imagesc(hyp0);colorbar
   subplot(3, 2, 5);mesh(hyp1)
   subplot(3, 2, 6);imagesc(hyp1);colorbar
end;

n= 8;
lar= 2;
for mar= 2:6
   resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im);
   load(resultname)

   HYP0= sum(HYP(:, :, ESTI, :)>lambda0B, 4);
   HYP1= sum(HYP(:, :, ESTI, :)>       1, 4);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE(:, :, ESTI) , smooth);
   hyp0= hyp0(2:end-1, 2:end-1);
   hyp1= hyp1(2:end-1, 2:end-1);
   mse = mse (2:end-1, 2:end-1);
   
%   mse= MSE(:, :, ESTI);
%   hyp0= HYP0;
%   hyp1= HYP1;
   
   [wing, winh]= find(mse==min(min(mse(2:end, 2:end))));
   [2^n mar XX(wing) lar YY(winh) lambda0B MSE(1, 1) MSE(wing, winh) HYP0(1, 1) HYP0(wing, winh) HYP1(1, 1) HYP1(wing, winh)]
%   [mar lar XX(wing) YY(winh) XX(end) YY(end) min(min(MSE))]
   MLopt(mar, lar, 1)= XX(wing);
   MLopt(mar, lar, 2)= YY(winh);
   figure(mar);clf;
   subplot(3, 2, 1);mesh(mse);
   subplot(3, 2, 2);imagesc(mse);colorbar
   subplot(3, 2, 3);mesh(hyp0)
   subplot(3, 2, 4);imagesc(hyp0);colorbar
   subplot(3, 2, 5);mesh(hyp1)
   subplot(3, 2, 6);imagesc(hyp1);colorbar
end

n= 8;
mar= 3;
for lar= 1:4
   resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im);
   load(resultname)

   HYP0= sum(HYP(:, :, ESTI, :)>lambda0B, 4);
   HYP1= sum(HYP(:, :, ESTI, :)>       1, 4);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE(:, :, ESTI) , smooth);
   hyp0= hyp0(2:end-1, 2:end-1);
   hyp1= hyp1(2:end-1, 2:end-1);
   mse = mse (2:end-1, 2:end-1);
   
%   mse= MSE(:, :, ESTI);
%   hyp0= HYP0;
%   hyp1= HYP1;
   
   [wing, winh]= find(mse==min(min(mse(2:end, 2:end))));
   [2^n mar XX(wing) lar YY(winh) lambda0B MSE(1, 1) MSE(wing, winh) HYP0(1, 1) HYP0(wing, winh) HYP1(1, 1) HYP1(wing, winh)]
%   [mar lar XX(wing) YY(winh) XX(end) YY(end) min(min(MSE))]
   MLopt(mar, lar, 1)= XX(wing);
   MLopt(mar, lar, 2)= YY(winh);
   figure(lar);clf;
   subplot(3, 2, 1);mesh(mse);
   subplot(3, 2, 2);imagesc(mse);colorbar
   subplot(3, 2, 3);mesh(hyp0)
   subplot(3, 2, 4);imagesc(hyp0);colorbar
   subplot(3, 2, 5);mesh(hyp1)
   subplot(3, 2, 6);imagesc(hyp1);colorbar
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%