function [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, crit)% function [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, 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 TFARMA(M, L) models. The% maximum TFAR search space defaults to MMAX, LMAX= N/8, The order% estimation IC defaults to the AIC. % TODO: M/LMAX even%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N    = 256;MAR  =   0;LAR  =   2;MMA  =   3;LMA  =   2;re_im= 'r';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;%Psi= ones(N);%Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);crit= 'AIC';MMAX= 4*max(MAR, MMA);LMAX= 4*max(LAR, LMA);x= tfarma_gen(randn(N, 1), Aml, Bml, beta);[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, '1.5');[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX);[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX);[AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x);%nargin= 4;MMAX= 5;LMAX= 5;MM= 100;NRITER= zeros(MMAX, LMAX+1);for mm= 1:MM   x= tfarma_gen(randn(N, 1), Aml, Bml, beta);   [AmlEst, BmlEst, SigmaHat, NrIter]= tfarma(x, MMAX, LMAX, 'MDL');   [MA, LA]= param_dim(AmlEst);   [MB, LB]= param_dim(BmlEst);   [mm MA LA MB LB]   NRITER= NRITER + NrIter(:, :, 1);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The nargin sory%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(nargin<=3)   crit= 'AIC';end;if(nargin<=2)   LMAX= floor(N/8);end;if(nargin<=1)   MMAX= floor(N/8);   LMAX= floor(N/8);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Some constants, etc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%alpha = 1/2;beta  = 1/2;N     = length(x);lambda= .9800;rho   = 1-log(12);switch(lower(crit)) case 'mdl'  PENALTY= log(N+1)+rho; case 'aic'  PENALTY= 2; otherwise  PENALTY= str2num(crit);end;MMAX= min(MMAX, floor(N/8));LMAX= min(LMAX, floor(N/8));SigmaHat= 9999999999*ones(MMAX, LMAX+1, 3);NrIter  = zeros(MMAX, LMAX+1, 3);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonparametric signal statistics%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ax= fft(corr_est(x, x, -1, 1/2));Ax= [Ax(N/2+1:N, :); Ax(1:N/2, :)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TFAR Identification (1:Mmax)x(0:LMAX)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sprintf('Searching up to TFAR(%d, %d) using the "%s"', MMAX, LMAX, crit)TFAR= {};TFARunstab= {};for L= 0:LMAX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% 127/32/32%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   if(~L)      Mmax= N/2-1;   else      Mmax= min(floor(N/4/L), MMAX);   end;%if(~L)   [L Mmax]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   warning off MATLAB:divideByZero   Mmax= min(floor(N/4/L), MMAX);   warning on MATLAB:divideByZero%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TAPER DESIGN (CANCEL???)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   [Psi, Mask, v2]= tf_multiwin(N, 3*Mmax, 4*L, 0, 2, 1);   v2   subplot(1, N/4+1, L+1)   tf_show(Psi);drawnow%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   [ATFAR, BTFAR]= tfar_est_tfywu(Ax(N/2+1-3*L:N/2+1+3*L, N/2-Mmax+1:N/2+Mmax+1), N);   for M= 1:Mmax      [Aml1, lambdamax, mmTFAR, nrFDIR]= param_stabilize(ATFAR(:, 1:M+1, M), N, lambda);      [M L mmTFAR];      B= BTFAR(:, :, M);      TFAR{M, L+1}= {Aml1 B};      TFARunstab{M, L+1}= {Aml1 B};      eprime= tfarma_inv(x, Aml1, B/param_get(B, 0, 0));      SigmaHat(M, L+1, 1)= eprime'*eprime/N;      NrIter(M, L+1, 1)= mmTFAR;   end;end;%for L= 0:LMAXNATFAR= (1:MMAX)'*(2*(0:LMAX)+1);NBTFAR= ones(MMAX, 1)*(2*(0:LMAX)+1);ICTFAR= log(SigmaHat(:, :, 1)) + PENALTY*(NATFAR+NBTFAR)/N;[MTFAR, LTFAR] = find(ICTFAR==min(min(ICTFAR)));LTFAR= LTFAR-1;AmlTFAR= TFAR{MTFAR, LTFAR+1}{1};BmlTFAR= TFAR{MTFAR, LTFAR+1}{2};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TFMA Identification (1:Mmax/2)x(0:LMAX/2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sprintf('Searching up to TFMA(%d, %d) using the "%s"', MMAX/2, LMAX/2, crit)TFMA= {};TFMAunstab= {};for L= 0:LMAX/2   warning off MATLAB:divideByZero   Mmax= min(floor(N/8/L), MMAX/2);   warning on MATLAB:divideByZero   for M= 1:Mmax      BB= tfma_est_dmu(TFAR{2*M, 2*L+1}{1}, TFAR{2*M, 2*L+1}{2}, M, L, N);      [Bml1, lambdamax, mmTFMA, nrFDIR]= param_stabilize(BB(:, 1:M+1, M+1), N, lambda);      [M L mmTFMA];      TFMA{M, L+1}= Bml1;      TFMAunstab{M, L+1}= BB(:, 1:M+1, M+1);            eprime= tfarma_inv(x, 1, Bml1/param_get(Bml1, 0, 0));      SigmaHat(M, L+1, 2)= eprime'*eprime/N;      NrIter(M, L+1, 2)= mmTFMA;   end;end;NATFMA= zeros(MMAX/2, 1)*(2*(0:LMAX/2)+1);NBTFMA= (1:MMAX/2)'*(2*(0:LMAX/2)+1);ICTFMA= log(SigmaHat(1:Mmax, 1:LMAX/2+1, 2)) + PENALTY*(NATFMA+NBTFMA)/N;[MTFMA, LTFMA] = find(ICTFMA==min(min(ICTFMA)));LTFMA= LTFMA-1;AmlTFMA= [zeros(LTFMA, 1); 1; zeros(LTFMA, 1)];BmlTFMA= TFMA{MTFMA, LTFMA+1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TFARMA Identification (1:Mmax/2)x(0:LMAX/2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sprintf('Searching up to TFARMA(%d, %d) using the "%s"', MMAX/2, LMAX/2, crit)TFARMA= {};TFARMAunstab= {};for L= 0:LMAX/2   warning off MATLAB:divideByZero   Mmax= min(floor(N/8/L), MMAX/2);   warning on MATLAB:divideByZero   for M= 1:Mmax      [AA, BB]= tfarma_est_gkmu(TFAR{2*M, 2*L+1}{1}, TFAR{2*M, 2*L+1}{2}, M, L, M-1, L, N);      [Bml1, lambdamax, mmTFMA, nrFDIR]= param_stabilize(BB(:, 1:M+1, M+1), N, lambda);      [M L mmTFMA];      TFMA{M, L+1}= Bml1;      TFMAunstab{M, L+1}= BB(:, 1:M+1, M+1);            eprime= tfarma_inv(x, 1, Bml1/param_get(Bml1, 0, 0));      SigmaHat(M, L+1, 2)= eprime'*eprime/N;      NrIter(M, L+1, 2)= mmTFMA;   end;end;NATFMA= zeros(MMAX/2, 1)*(2*(0:LMAX/2)+1);NBTFMA= (1:MMAX/2)'*(2*(0:LMAX/2)+1);ICTFMA= log(SigmaHat(1:Mmax, 1:LMAX/2+1, 2)) + PENALTY*(NATFMA+NBTFMA)/N;[MTFMA, LTFMA] = find(ICTFMA==min(min(ICTFMA)));LTFMA= LTFMA-1;AmlTFMA= [zeros(LTFMA, 1); 1; zeros(LTFMA, 1)];BmlTFMA= TFMA{MTFMA, LTFMA+1};%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Model Type Selection%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%P0= tfarma_wvsp(Aml, Bml, N, alpha);P1= tfarma_wvsp(AmlTFAR, BmlTFAR, N, alpha);P2= tfarma_wvsp(AmlTFMA, BmlTFMA, N, alpha);%figure(1);tf_show(P0)figure(2);tf_show(P1)figure(3);tf_show(P2)if(ICTFAR(MTFAR, LTFAR+1)<ICTFMA(MTFMA, LTFMA+1))   sprintf('TFAR!')   AmlEst= AmlTFAR;   BmlEst= BmlTFAR;else   sprintf('TFMA!')   AmlEst= AmlTFMA;   BmlEst= BmlTFMA;end;