tfma_papera.m

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

function tfma_papera(re_im, MM)
% function tfma_papera(re_im, MM)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates re_im-valued TFMA(MMA, LMA; N; 1/2) models. MM
% realizations are done to estimate the MSE of the parameter
% functions. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
re_im= 'i';
MM   = 50;

mm = 1;
n  = 6;
MMA= 4;
LMA= 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MAR  =   0;
LAR  =   0;
mo_no= 'n';
alpha= 1/2;
beta = 1/2;

xlabelstring= 'wing';
ylabelstring= 'winh';

for n= 8:8
   N= 2^n;
   for MMA= 4:4
      for LMA= 0:3

         tfpm_file_gen;
%---------------------
         for WIN_M= 1:1
            for WIN_L= 2:2
               resultname= sprintf('data/%04d/tfma_ceps_papera%d%d%d%d%s.mat', ...
	                           N, MMA, LMA, WIN_M, WIN_L, re_im)
               if(exist(resultname)~=2)
                  estimators= [
                   'WCm';%CTFC estimator, broad ambiguity  , tfma_est_cepsm.m
                   'WCl';%CTFC estimator, broad ambiguity  , tfma_est_cepsl.m
                   'KAD';%CTVC estimator of Kaderli et. al., tvma_est_cepsb.m
                  ];

                  MSEml= zeros(size(estimators, 1), 1);VARml= MSEml;BIQml= MSEml;
                  MSEmn= zeros(size(estimators, 1), 1);VARmn= MSEmn;BIQmn= MSEmn;

                  Psi= tf_multiwin(N, WIN_M*MMA, WIN_L*LMA, 6, min(3, N/32), 2);
%                  Psi= tf_window(N, 0, 0, WIN_L*LMA, WIN_M*MMA);
	 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
norm(Psi-sym_ambi(Psi))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                  BML_WCm= zeros(2*LMA+1, MMA+1, MM);
                  BML_WCl= BML_WCm;
                  BML_KAD= BML_WCm;

                  BMN_WCm= zeros(N      , MMA+1, MM);
                  BMN_WCl= BMN_WCm;
                  BMN_KAD= BMN_WCm;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                     sprintf('tfma_ceps(%2d, %2d; %4d) %s mm= %3d', MMA, LMA, N, re_im, mm)
                     e= randn(N, 1);
                     y= tfarma_gen(e, Aml, Bml, beta);
                     aalpha= ambi_est_cyc(y);
                     Aalpha= aalpha.*conj(Psi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
norm(aalpha-sym_ambi(aalpha))
norm(Aalpha-sym_ambi(Aalpha))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                     Ralpha= ml_to_nm(Aalpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(98);cesh(Aalpha)
figure(99);cesh(Ralpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                     BmlWCm= tfma_est_cepsm(Aalpha, MMA, LMA);
                     BmlWCl= tfma_est_cepsl(Aalpha, MMA, LMA);
                     BmnKAD= tvma_est_cepsb(Ralpha, MMA     );

                     BmnWCm= param_expand(BmlWCm, N, 'exp');
                     BmnWCl= param_expand(BmlWCl, N, 'exp');
                     BmlKAD= param_compress(BmnKAD, LMA);

                     BML_WCm(:, :, mm)= BmlWCm;
                     BML_WCl(:, :, mm)= BmlWCl;
                     BML_KAD(:, :, mm)= BmlKAD;
                     BMN_WCm(:, :, mm)= BmnWCm;
                     BMN_WCl(:, :, mm)= BmnWCl;
                     BMN_KAD(:, :, mm)= BmnKAD;
	    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  end;%for mm= 1:MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                  [M_BWCm, V_BWCm, B_BWCm]= param_mse(BML_WCm, Bml);
                  [M_BWCl, V_BWCl, B_BWCl]= param_mse(BML_WCl, Bml);
                  [M_BKAD, V_BKAD, B_BKAD]= param_mse(BML_KAD, Bml);

                  MSEml(1)= sum(sum(M_BWCm))/energy(Bml);
                  MSEml(2)= sum(sum(M_BWCl))/energy(Bml);
                  MSEml(3)= sum(sum(M_BKAD))/energy(Bml);

                  VARml(1)= sum(sum(V_BWCm))/energy(Bml);
                  VARml(2)= sum(sum(V_BWCl))/energy(Bml);
                  VARml(3)= sum(sum(V_BKAD))/energy(Bml);

                  BIQml(1)= sum(sum(B_BWCm))/energy(Bml);
                  BIQml(2)= sum(sum(B_BWCl))/energy(Bml);
                  BIQml(3)= sum(sum(B_BKAD))/energy(Bml);
                  MSEml
	 
                  [M_BWCm, V_BWCm, B_BWCm]= param_mse(BMN_WCm, Bmn);
                  [M_BWCl, V_BWCl, B_BWCl]= param_mse(BMN_WCl, Bmn);
                  [M_BKAD, V_BKAD, B_BKAD]= param_mse(BMN_KAD, Bmn);

                  MSEmn(1)= sum(sum(M_BWCm))/energy(Bmn);
                  MSEmn(2)= sum(sum(M_BWCl))/energy(Bmn);
                  MSEmn(3)= sum(sum(M_BKAD))/energy(Bmn);
		  
                  VARmn(1)= sum(sum(V_BWCm))/energy(Bmn);
                  VARmn(2)= sum(sum(V_BWCl))/energy(Bmn);
                  VARmn(3)= sum(sum(V_BKAD))/energy(Bmn);

                  BIQmn(1)= sum(sum(B_BWCm))/energy(Bmn);
                  BIQmn(2)= sum(sum(B_BWCl))/energy(Bmn);
                  BIQmn(3)= sum(sum(B_BKAD))/energy(Bmn);
                  MSEmn
 	 
                  save(resultname, 'estimators', 'xlabelstring', 'ylabelstring', ...
                                   'MSEml', 'VARml', 'BIQml', ...
                                   'MSEmn', 'VARmn', 'BIQmn', ...
                                   'MM');
               end 
            end
         end
      end
   end
end












%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   0;
LAR  =   0;
mo_no= 'n';
alpha= 1/2;
beta = 1/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nmin= 5;
nmax= 9;

MMA= 4;
LMA= 2;

s= 'm';




while s~='x'
for WIN_M= 1:1
   for WIN_L= 2:2

%                   EST  REIM N
      MSENml= zeros(3,   2,  nmax);BIQNml= MSENml;VARNml= MSENml;
      MSENmn= zeros(3,   2,  nmax);BIQNmn= MSENmn;VARNmn= MSENmn;

      for n= nmin:nmax
         N= 2^n;

%         resultnamer= sprintf('data/%04d/tfma_ceps_papera%d%d%d%d%s.mat', ...
%                              N, MMA, LMA, WIN_M, WIN_L, 'r')
%         load(resultnamer)
%         MSENml(:, 1, n)= MSEml;
%         BIQNml(:, 1, n)= BIQml;
%         VARNml(:, 1, n)= VARml;
%         MSENmn(:, 1, n)= MSEmn;
%         BIQNmn(:, 1, n)= BIQmn;
%         VARNmn(:, 1, n)= VARmn;

         resultnamei= sprintf('data/%04d/tfma_ceps_papera%d%d%d%d%s.mat', ...
                              N, MMA, LMA, WIN_M, WIN_L, 'i')
         load(resultnamei)
         MSENml(:, 2, n)= MSEml;
         BIQNml(:, 2, n)= BIQml;
         VARNml(:, 2, n)= VARml;
         MSENmn(:, 2, n)= MSEmn;
         BIQNmn(:, 2, n)= BIQmn;
         VARNmn(:, 2, n)= VARmn;
      end;

%% Plot MSE, BIQ, VAR over N