simu_tfmawin.m

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

function simu_tfmawin(MMA, LMA, N, re_im, MM)
% function simu_tfmawin(MMA, LMA, 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 TFMA(<MMA>, <LMA>; <N>; 1/2) models to find
% optimum window lengthes. <MM> realizations are done to estimate the
% MSE. Uses the estimator tfma_est_cepsb.m. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MMA  =   3;
LMA  =   2;
re_im= 'i';
MM   =   3;
xx= 1;
yy= 1;
mm= 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MAR  =   0;
LAR  =   0;
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;

[TDIRnm, TDIRn0]= param_tdir(param_expand(Bml, N));
lambda0= max(max(abs(TDIRnm)));

exp= 'tfmawin';
resultname= sprintf('data/%04d/%s%d%d%s.mat', N, exp, MMA, LMA, re_im)

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

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   for yy= 1:yymax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      BML_cep= zeros(2*LMA+1, MMA+1, MM);
      Psi= 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         [N MMA LMA xx 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+1-3*LAR:N/2+1+3*LAR, :);
%         figure(3);mesh(abs(ayy));drawnow
%         figure(4);mesh(abs(ayy.*psi));drawnow
         B= tfma_est_cepsb(ayy.*conj(Psi), MMA, LMA);
         BML_cep(:, :, mm)= B;
	 [TDIRnm, TDIRn0]= param_tdir(param_expand(B, N));
         HYP(xx, yy, mm)= max(max(abs(TDIRnm)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end;%for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [Mcep, Vcep, Bcep]= param_mse(BML_cep, Bml);
      MSE(xx, yy)= sum(sum(Mcep))/sum(sum(abs(Bml).^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_tfmawin(3, 2, 2^n, 'i', 20);
end;

clear;tfpm;
for mar= 1:2
   simu_tfmawin(mar, 2, 256, 'i', 20);
end;

clear;tfpm;
for mar= 4:6
   simu_tfmawin(mar, 2, 256, 'i', 20);
end;

clear;tfpm;
for lar= 0:1
   simu_tfmawin(3, lar, 256, 'i', 20);
end;

clear;tfpm;
for lar= 3:4
   simu_tfmawin(3, lar, 256, 'i', 20);
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= 'tfmawin';
re_im= 'i';

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:4
   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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%