simu_tfarmawin.m

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

function simu_tfarmawin(MAR, LAR, MMA, LMA, N, re_im, MM)
% function simu_tfarmawin(MAR, LAR, MMA, LMA, N, re_im, MM)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MAR  =   3;
LAR  =   2;
MMA  = MAR-1;
LMA  = LAR;
re_im= 'r';
MM   = 3;
xx= 1;
yy= 1;
mm= 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mo_no= 'n';
N0= N;
N= 64;
tfpm_file_gen;
%-------------
N= N0;
%param_plot('exp', N, 0, 'r', 'g', Aml, 1, 1, 0, 3, 0, 0, Bml, 0, 0, 0, 0, 0, 0)
beta = 1/2;


exp= 'tfarmawin';
resultname= sprintf('data/%04d/%s%d%d%s.mat', N, exp, MAR, LAR, re_im)

% Window lengthes (wing)
XX= 0:2:6*MAR

% Window lengthes (winh)
YY= 0:2:6*LAR

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for xx= 1:xxmax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   for yy= 1:yymax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      AML_ywu= zeros(2*LAR+1, MAR+1, MM);
      BML_cep= zeros(2*LMA+1, MMA+1, MM);
      [Psi, Mask, v2]= tf_multiwin(N, XX(xx), YY(yy), 0, 2, 1);
%      figure(1);mesh(abs(Psi));drawnow
%      figure(2);mesh(abs(tfarma_ambi(Aml, Bml, N, 1/2)));drawnow
      
%      psi= Psi(N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);
      tf_show(Psi)
      psi= psi(N/2-3*LAR+1:N/2+3*LAR+1, N/2+MMA-MAR+1:N/2+MMA+MAR+1);
      text(10, 10, sprintf('%d', prod(size(v2))))
      drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         [N MAR LAR MMA LMA xx length(XX) yy length(YY) mm]
         e= randn(N, 1);
         y= tfarma_gen(e, Aml, Bml, beta);
         qyy= corr_est(y, y, -1, 1/2);
         ayy= fft(qyy);
         ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
         ayy= ayy(N/2-3*LAR+1:N/2+3*LAR+1, N/2+MMA-MAR+1:N/2+MMA+MAR+1);
	 
         AMLu= tfarma_est_tfywu(ayy.*psi, MMA, N);
         AML_ywu(:, :, mm)= [AMLu(:, 2:end, MAR) B0Lu(:, :, MAR)];
	 [TDIRnm, TDIRn0]= param_tdir(param_expand(AMLu(:, :, MAR), N));
         HYP(xx, yy, mm)= max(max(abs(TDIRnm)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end;%for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [Mywu, Vywu, Bywu]= param_mse(AML_ywu, [Aml(:, 2:end) B0l]);
      MSE(xx, yy)= sum(sum(Mywu))/sum(sum(abs([Aml(:, 2:end) B0l]).^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', 'lambda0', 'HYP')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
for n= 6:9
   simu_tfarmawin(3, 2, 2, 2, 2^n, 'r', 10);
end;

clear;tfpm;
for mar= 1:2
   simu_tfarmawin(mar, 2, mar-1, 2, 256, 'r', 10);
end;

clear;tfpm;
for mar= 4:6
   simu_tfarmawin(mar, 2, mar-1, 2, 256, 'r', 10);
end;

clear;tfpm;
for lar= 0:1
   simu_tfarmawin(3, lar, 2, 2, 256, 'r', 10);
end;

clear;tfpm;
for lar= 3:4
   simu_tfarmawin(3, lar, 2, 2, 256, 'r', 10);
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= 'tfarwin';
re_im= 'r';

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

%NULL SEARCH!!
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>lambda0, 3);
   HYP1= sum(HYP>      1, 3);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE , 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);
   
   [wing, winh]= find(mse==min(min(mse)));
   [2^n mar XX(wing) lar YY(winh) lambda0 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(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= 1:6
   resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im);
   load(resultname)

   HYP0= sum(HYP>lambda0, 3);
   HYP1= sum(HYP>      1, 3);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE , 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);
   
   [wing, winh]= find(mse==min(min(mse)));
   [2^n mar XX(wing) lar YY(winh) lambda0 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:3
   resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im);
   load(resultname)

   HYP0= sum(HYP>lambda0, 3);
   HYP1= sum(HYP>      1, 3);

   hyp0= conv2(HYP0, smooth);
   hyp1= conv2(HYP1, smooth);
   mse = conv2(MSE , 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);
   
   [wing, winh]= find(mse==min(min(mse)));
   [2^n mar XX(wing) lar YY(winh) lambda0 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



%%FULL SEARCH (NOT DONE!!)
for n= 6:8
   MLopt= zeros(MARmax, LARmax, 2);
   for mar= 1:MARmax
      for lar= 1:LARmax
         resultname= sprintf('data/%04d/%s%d%d%s.mat', 2^n, exp, mar, lar, re_im)
	 load(resultname)
	 [wing, winh]= find(MSE==min(min(MSE)));
	 [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(1);clf;contour(MSE);
%	 colormap(gray);
%	 colorbar
%         set(gca, 'XTick', 1:length(YY))
%         set(gca, 'XTicklabel', YY)
%         set(gca, 'YTick', 1:length(XX))
%         set(gca, 'YTicklabel', XX)
%	 text(wing, winh, '*') 
%	 text(mar, lar, '+';)

         figure(1);clf;mesh(MSE);
	 colormap(flipud(gray));
	 colorbar
         figure(2);clf;imagesc(MSE);
	 colormap(flipud(gray));
	 colorbar
         set(gca, 'XTick', 1:length(YY))
         set(gca, 'XTicklabel', YY)
         set(gca, 'YTick', 1:length(XX))
         set(gca, 'YTicklabel', XX)
	 text(winh, wing, '*') 
	 text(lar, mar, '+') 
	 
%         figure(3);mesh(MSE);colorbar
         drawnow
	 title(sprintf('%d, %d, %d', mar, lar, 2^n))
	 pause
      end;
   end;
   figure(n);clf
   subplot(2, 1, 1);
   plot(MLopt(:, :, 1))
   legend('1', '2', '3')
   subplot(2, 1, 2);
   plot(MLopt(:, :, 2)')
   legend('1', '2', '3', '4', '5')
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%