function [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA)% function [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Estimates TFARMA(MAR, LAR; MMA, LMA; N; alpha) models using the% method of Graupe/Krause/Moore (IEEEAC 20, feb75, pp104). In this% method: MAR>MMA, LAR>LMA! The ambiguity function Ayy needs to be% given up to MCmax, the maximum M-dimension of the long underlying% TFAR model. Its L-dimension is LC. The estimator tfar_est_tfywe.m% is used in here. INSTAB is zero, if the long TFAR model in stable% (no pole magnitude > 1) otherwise the maximum value of the impulse% response. % MC, LC >= ??%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;MAR  =   1;LAR  =   1;MMA  =   1;LMA  =   1;N    = 256;re_im= 'r';mo_no= 'n';tfpm_file_gen;%-------------beta = 1/2;alpha= 1/2;KC=  2;LC= KC*(LAR+LMA);MC= KC*(MAR+MMA);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, :);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%INSTAB= 0;alpha= 1/2;% DimensionsMCmax= (size(Ayy, 2)-1)/2;LC   = (size(Ayy, 1)-1)/6;LL   = 2*LC+1;% Estimate the long TFAR model[Cml, CB]= tfar_est_tfywe(Ayy, N);%param_plot('exp', N, 0, 'r', 'g', Cml, 1, 1, 2, 3, 4, 5, 1, 0, 0, 0, 0, 0, 0);[TDIRnm, TDIRn0]= param_tdir(param_expand(Cml, N));if(max(max(abs(TDIRnm))) > 1)   HAR= tfarma_impr(Cml, 1, N, alpha);   INSTAB= max(max(abs(HAR)));end;   % The TFMA partCCC= [      Cml(LC-3*LMA+1:LC+3*LMA+1, MAR-MMA+1:MAR+MMA+1)     ];[B, B0L]= tfar_est_tfywe(CCC, N);% The TFAR partB1= [ [zeros(LC-LMA, 1+MMA); B; zeros(LC-LMA, 1+MMA)] zeros(2*LC+1, MCmax-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)*...                       exp(2*pi*j/N*mprime*(l-lprime))*...                       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)INSTAB%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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  =   1;LAR  =   1;MMA  =   1;LMA  =   1;N    = 256;re_im= 'i';mo_no= 'm';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;KC=  2;LC= KC*(LAR+LMA);MC= KC*(MAR+MMA);MM= 100;Psi= tf_multiwin(N, MC, LC, 5, 1, 1);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);   [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA);   [mm param_norm(A, Aml, 0) param_norm(B, Bml, 0) INSTAB]   AA= AA+A;   BB= BB+B;end;AA/MMBB/MM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%