function [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N)% function [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % MONIC TFARMA SYSTEM!!!!!% MC, LC >= ??%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;MAR  =   2;LAR  =   2;MMA  =   1;LMA  = LAR;N    = 256;re_im= 'r';mo_no= 'm';tfpm_file_gen;%-------------beta = 1/2;alpha= 1/2;Bml= .2*Bml;KC=  4;MC= 2*MAR-1;LC= 2*LAR;Aalpha= tfarma_ambi(Aml, Bml, N, alpha);Ayy= Aalpha(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1);%y= tfarma_gen(randn(N, 1), Aml, Bml, beta);%qyy= corr_est(y, y, MC, 1/2);%ayy= fft(qyy);%ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];%Ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :);% Estimate the long TFAR model[CA, CB]= tfar_est_tfywu(Ayy, N);Cml= CA(:, :, end);D0l= CB(:, :, end);%Cml= ones(2*LC+1, 1)*(0:MC) + j*(-LC:LC)'*ones(1, MC+1);D0l= -1;P0= tfarma_wvsp(Aml, Bml, N, 1/2);figure(1);tf_show(P0)P1= tfarma_wvsp(Cml, D0l, N, 1/2);figure(2);tf_show(P1)figure(3);tf_show(abs(P0-P1))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%alpha= 1/2;% Dimensions[MC, LC]= param_dim(Cml);LL   = 2*LMA+1;% Fill the Toeplitz/block-Toeplitz generator rPMrPM= [];for mu= MMA-1:-1:-MMA+1   rPM= [rPM; fliplr(Cml(LC-2*LMA+1:LC+2*LMA+1, MAR+mu+1).')];end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%toeptoep(rPM)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RHS= [];for m= MAR+1:MAR+MMA   RHS= [RHS; -Cml(LC-LMA+1:LC+LMA+1, m+1)*param_get(D0l, 0, 0)];end;% Solve the systemBBstack= invert_waxkailath(rPM, RHS);B= [[zeros(LMA, 1); param_get(D0l, 0, 0); zeros(LMA, 1)] reshape(BBstack(:, :, end), LL, MMA)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Reshape to BBBB= zeros(2*LMA+1, 1+MMA, MMA);for m= 1:MMA   BB(:, :, m)= [[zeros(LMA, 1); param_get(D0l, 0, 0); zeros(LMA, 1)] reshape(BBstack(1:LL*m, 1, m), LL, m) zeros(LL, MMA-m)];end;param_norm(BB(:, :, end), Bml, 1)param_norm(B, Bml, 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The TFAR partB1= [ [zeros(LC-LMA, 1+MMA); B/param_get(D0l, 0, 0); zeros(LC-LMA, 1+MMA)] zeros(2*LC+1, MC-MMA) ];A= zeros(2*LAR+1, 1+MAR);A(LAR+1, 1)= 1;for m= 1:MAR   for l= -LAR:LAR      aml= Cml(LC+1+l, 1+m);      for mprime= 1:m         for lprime= -LAR:LAR            aml= aml + Cml(LC+1+l-lprime, 1+m-mprime)*...                       B1(LC+1+lprime, 1+mprime);         end;      end;      A(LAR+1+l, 1+m)= aml;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%param_norm(A, Aml, 0)param_norm(B, Bml, 0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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;MAR  =   3;LAR  =   2;MMA  =   2;LMA  =   2;N    = 256;re_im= 'i';mo_no= 'm';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;mm   =   1;KC=  2;LC= KC*(LAR);MC= KC*(MAR);MM= 1000;%Psi= tf_multiwin(N, MC, LC, 5, 1, 1);Psi= ones(N);psi= Psi(N/2+1-3*LC:N/2+1+3*LC, N/2+1-MC:N/2+1+MC);AA= zeros(size(Aml));BB= zeros(size(Bml));for mm= 1:MM   y= tfarma_gen(randn(N, 1), Aml, Bml, beta);   qyy= corr_est(y, y, MC, 1/2);   ayy= fft(qyy);   ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];   Ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :).*conj(psi);% Estimate the long TFAR model   [CA, CB]= tfar_est_tfywu(Ayy, N);   Cml= CA(:, :, end);   D0l= CB(:, :, end);   [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N);   [mm param_norm(A, Aml, 0) param_norm(B, Bml, 0)]   AA= AA+A;   BB= BB+B;end;AA/MMAmlBB/MMBml[param_norm(AA/MM, Aml, 0) param_norm(BB/MM, Bml, 0)]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%