function [E, Eprime]= tfarma_fit2(x, Mmax, Lmax, Psi, lambda)% function [E, Eprime]= tfarma_fit2(x, Mmax, Lmax, Psi, lambda)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Searches TFARMA(M, L, M-1, L) models for M= 1:Mmax, L= 0:Lmax. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N    = 256;MAR  =   2;LAR  =   2;MMA  =   1;LMA  =   2;re_im= 'r';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;lambda= .9800;Mmax= max([2*MAR 2]);Lmax= max([2*LAR 2]);%Psi= ones(N);Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);x= tfarma_gen(randn(N, 1), Aml, Bml, beta);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%alpha= 1/2;beta = 1/2;N= length(x);Ax= fft(corr_est(x, x, -1, alpha));Ax= [ Ax(N/2+1:N, :); Ax(1:N/2, :)].*conj(Psi);E= 99999*ones(Mmax, Lmax+1, Lmax+1);Eprime= 99999*ones(Mmax, Lmax+1, Lmax+1);NA= (1:Mmax)'*(2*(0:Lmax)+1);NB= NA - 1;for L= 0:Lmax   M= 1;      AA= tfarma_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2+(M-1)-M+1:N/2+(M-1)+M+1), M-1, N);      [Aml1, lambdamax, mm, nrFDIR]= param_stabilize(AA(:, 1:M+1, M), N, lambda);      Aml1= Aml1./param_get(Aml1, 0, 0);      Aml1(1:L, 1)= zeros(L, 1);      Aml1(L+2:end, 1)= zeros(L, 1);      uM= tfarma_inv(x, Aml1, 1);      Au= fft(corr_est(uM, uM, -1, alpha));      Au= [ Au(N/2+1:N, :); Au(1:N/2, :) ].*conj(Psi);      [B, REG, PMIN]= tfma_est_cepsb(Au, M-1, L);      [Bml1, lambdamax, mm, nrFDIR]= param_stabilize(B, N, lambda);% Inverse TFMA Filtering      eprime= tfarma_inv(uM, 1, Bml1/param_get(Bml1, 0, 0));      e     = tfarma_inv(uM, 1, Bml1);      E(M, L+1, 1)     = e'*e/N;      Eprime(M, L+1, 1)= eprime'*eprime/N;      for M= 2:Mmax      AA= tfarma_est_tfywu(Ax(N/2-3*L+1:N/2+3*L+1, N/2+(M-1)-M+1:N/2+(M-1)+M+1), M-1, N);      for LL= 0:L% Inverse TFAR Filtering         Mask= ones(2*L+1, M+1);         Mask(1:L-LL, end)= zeros(L-LL, 1);         Mask(L+2+LL:end, end)= zeros(L-LL, 1);%         figure(1);plot(abs(param_tdir(AA(:, 1:M+1, M), N)))%         figure(2);plot(abs(param_tdir(AA(:, 1:M+1, M).*Mask, N)))         [Aml1, lambdamax, mm, nrFDIR]= param_stabilize(AA(:, 1:M+1, M).*Mask, N, lambda);         Aml1= Aml1./param_get(Aml1, 0, 0);         Aml1(1:L, 1)= zeros(L, 1);         Aml1(L+2:end, 1)= zeros(L, 1);%         figure(3);plot(abs(param_tdir(Aml1, N)))         uM= tfarma_inv(x, Aml1, 1);         Au= fft(corr_est(uM, uM, -1, alpha));         Au= [ Au(N/2+1:N, :); Au(1:N/2, :) ].*conj(Psi);         [B, REG, PMIN]= tfma_est_cepsb(Au, M-1, L);         [Bml1, lambdamax, mm, nrFDIR]= param_stabilize(B.*Mask(:, 2:end), N, lambda);% Inverse TFMA Filtering         eprime= tfarma_inv(uM, 1, Bml1/param_get(Bml1, 0, 0));         e     = tfarma_inv(uM, 1, Bml1);         E(M, L+1, L-LL+1)     = e'*e/N;         Eprime(M, L+1, L-LL+1)= eprime'*eprime/N;%	 Aml1%	 Bml1%	 drawnow      end;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;for n= 8:8   N= 2^n;   re_im= 'i';   mo_no= 'n';   alpha= 1/2;   beta = 1/2;   lambda= .9800;   MM= 100;   for M= 5:5      for L= 2:2         MAR  =   M;         LAR  =   L;         MMA= MAR-1;         LMA  = LAR;         tfpm_file_gen;%-------------------------         Bml= Bml/param_get(Bml, 0, 0);         Mmax= max([2*M 2]);         Lmax= max([2*L 2]);         [Psi, Mask, v2]= tf_multiwin(N, 2*Mmax+2, 2*Lmax+2, 5, 0, 0);%         Psi= ones(N);         EE= zeros(Mmax, Lmax+1, Lmax+1, MM);         EEprime= zeros(Mmax, Lmax+1, Lmax+1, MM);         for mm= 1:MM            x= tfarma_gen(randn(N, 1), Aml, Bml/param_get(Bml, 0, 0), 1/2);            [E, Eprime]= tfarma_fit2(x, Mmax, Lmax, Psi, lambda);	    EE(:, :, :, mm)= E;	    EEprime(:, :, :, mm)= Eprime;            [N MAR LAR mm]	 end;         exp= 'ord2C';         resultname= sprintf('data/%04d/%s%d%d%s%s.mat', N, exp, M, L, re_im, mo_no)         save(resultname, 'lambda', 'Mmax', 'Lmax', 'MM', 'EE', 'EEprime')      end;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%