function [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMIN, MMAX, LMIN, LMAX, BMIN, BMAX, crit)% function [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMIN, MMMAX, LMIN, LMAX, BMIN, BMAX, crit)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % An automatic procedure to identify VTFAR(M, L) models. The% maximum search space defaults to MMAX, LMAX= N/16, The order% estimation IC defaults to the AIC. The signal is NOT normalized% to energy N! % EstInfo= {SigHat, InfoCr, NrIter, NrPara, EstInn, PENALTY, AMBI};%% TODO: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmalpha= 1/2;N=  64;M0=  3;L0=  2;J0=  1;B0=  1;MMIN= 1;MMAX= M0+1;LMIN= 0;LMAX= L0+1;BMIN= 1;BMAX= J0;crit= 'AIC';filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)load(filename);AML0= AML;B0L0= B0L;%AML0= eye(J0);B0L0= eye(J0);E0= randn(N, J0);X= vtfarma_gen(E0, AML0, B0L0, B0, 1/2);L= L0;M= M0;[AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMAX, LMIN, LMAX, BMIN, BMAX, crit)y= tfarma_gen(E0, reshape(AML, 2*L+1, M+1), reshape(B0L, 2*L+1, 1), 1/2);norm(X-y)/norm(X)% THE CHANNELclear;tfpmalpha= 1/2;load updown_h_rusk.matfigure(1);clf;subplot(2, 1, 1);mesh(abs(h_I2))N= 512;J= 5;MMIN= 1;MMAX= 6;LMIN= 0;LMAX= 3;BMIN= 1;BMAX= 3;crit= 'AIC';for offset= 0:200:1800   X= 1000*h_I2(offset+1:offset+N, 29:29+J);   subplot(2, 1, 2);mesh(abs(X/1000));drawnow   [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMIN, MMAX, LMIN, LMAX, BMIN, BMAX, crit);   [Mest, Lest, Best]   filename= sprintf('ident%02d%02d%02d-%04d', J, BMIN, BMAX, offset)   save(filename, 'AML', 'B0L', 'InfoCr', 'Varian', 'Mest', 'Lest', 'Best');end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The penalty%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%switch(lower(crit)) case 'mdl'  PENALTY= log(N+1)+rho; case 'aic'  PENALTY= 2; otherwise  PENALTY= str2num(crit);end;[N, J]= size(X);alpha= 1/2;% Compute the EXAFsX= X/sqrt(energy(X));[AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for j= 1:J   for jp= 1:J      [j jp]      figure(1);clf;tf_show(reshape(AMBI(j, jp, :, :), N, N))      figure(2);clf;mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))      pause   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PrepareInfoCr= 9999999*ones(MMAX, LMAX+1, J);NrPara= InfoCr;Varian= InfoCr;for m= 1:MMAX   for l= 0:LMAX      for b= 1:J         NrPara(m, l+1, b)= ((m+1)*(.7*2*l+1)-1)*number_band(J, b);      end;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for L= LMIN:LMAX   for M= MMIN:MMAX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for L= 1:1%   for M= 2:2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      [M L]      [AML, B0L]= vtfar_parest(ambi, M, L, N);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[AML0 AML]sum(sum(sum(sum(abs(AML-AML0)))))/sum(sum(sum(sum(abs(AML0)))))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the MoNiC Model!!!!!!!!!!!!!      b0L= zeros(size(B0L));      bbb= inv(B0L(:, :, L+1, 1));      for l= -L:L         b0L(:, :, L+1+l, 1)= B0L(:, :, L+1+l, 1)*bbb;      end;      for B= BMIN:BMAX%         clear AMBI%	 pack%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         clear Linv         Linv= vtfarma_Linv(AML, b0L, N, B, alpha);         V= zeros(J);         for tau= 1:J            for taup= 1:J               for t= 1:J                  for tp= 1:J%                     [M L B tau taup t tp]                     V(tau, taup)= V(tau, taup) + sum(sum(...                     reshape(Linv(tau, t, :, :), N, N).*...                     reshape(conj(Linv(taup, tp, :, :)), N, N).*...                     nm_to_nk(ml_to_nm(reshape(AMBI(t, tp, :, :), N, N)))));                  end;               end;            end;         end;         V= V/N^2;         e= real(eig(V));         e= e.*(e>0);	 if(prod(e))	    Varian(M, L+1, B)= log(prod(e));	 else	    Varian(M, L+1, B)= 99999999;	 end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         Ehat= vtfarma_inv(X, AML, b0L, B, 1/2);         V= Ehat'*Ehat/N;	 Varian(M, L+1, B)= log(real(det(V)));      end;%for B= 0:J-1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   end;%for M= 1:MMAX;end;%for L= 0:LMAX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%InfoCr= Varian + PENALTY*NrPara/N;MinCr= max(max(max(InfoCr)));Best= 0;for b= 1:J   if(min(min(InfoCr(:, :, b)))<MinCr)      [Mest, Lest]= find(InfoCr(:, :, b)==min(min(InfoCr(:, :, b))));      MinCr= min(min(InfoCr(:, :, b)));      Best= b;   end;end;Lest= Lest-1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[M0 L0 B0 Mest Lest Best]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[AML, B0L]= vtfar_parest(ambi, Mest, Lest, N);for j= 1:J   for jp= 1:J      if(abs(j-jp)>=Best)         AML(j, jp, :, :)= zeros(2*Lest+1, Mest+1);         B0L(j, jp, :, :)= zeros(2*Lest+1, 1);      end;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmalpha= 1/2;N=  256;M0= 2;L0= 1;J0= 4;B0= 2;MMAX= M0+2;LMAX= L0+1;crit= 'AIC';filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)load(filename);AML0= AML;B0L0= B0L;MM= 100;MLBest= [];for mm= 1:MM   X= vtfarma_gen(randn(N, J0), AML0, B0L0, B0, 1/2);   [AML, B0L, Mest, Lest, Best, InfoCr, Varian]= vtfar20(X, MMAX, 0, LMAX, 1, J0, crit);   [mm Mest Lest Best M0 L0 B0]   MLBest= [MLBest; [Mest Lest Best]];end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%