function [A, B, INSTAB]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)% function [A, B, INSTAB]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)%   This file is part of the TFPM toolbox v0.3 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Estimates TFARMA(MAR, LAR; MMA, LMA; N; alpha) models WITH MONIC% TFMA PART 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_tfywu.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. %% IS FUCKING UNSTABLE!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N    =  64;MAR  =   3;LAR  =   3;MMA  =   2;LMA  = LAR;re_im= 'r';mo_no= 'n';tfpm_file_gen;%Bml(LMA+1, 1)= 6;%-------------beta = 1/2;alpha= 1/2;Aalpha= tfarma_ambi(Aml, Bml, N, alpha);KC=  2;LC= KC*(LAR+LMA)+1;MC= KC*(MAR+MMA)+1;Ayy= Aalpha(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1);%max(max(abs(Ayy)))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[CA, CB]= tfar_est_tfywu(Ayy, N);Cml= reshape(CA(:, :, MCmax), LL, MCmax+1);% A pure TFAR model?if(~MMA)   A= Cml(LC+1-LAR:LC+1+LAR, 1:MAR+1);   B= CB(LC+1-LMA:LC+1+LMA, 1, end);   return;end;% A pure TFMA model?%if(~MAR)%   A= nk_to_nm(nm_to_ml(1./param_weyl(Cml, N, 1/2) ))%end;% The TFMA partCCC= Cml(LC-3*LMA+1:LC+3*LMA+1, MAR-MMA+1:MAR+MMA+1);[BML, B0L]= tfar_est_tfywu(CCC, N);B= reshape(BML(:, :, MMA), 2*LMA+1, MMA+1)*CB(LC+1, 1, end);% 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)*...                       B1(LC+1+lprime, 1+mprime);         end;      end;      A(LAR+1+l, 1+m)= aml;   end;end;%Stable??TDIRA= param_tdir(param_expand(A, N));TDIRB= param_tdir(param_expand(B, N));if(max(max(abs([TDIRA TDIRB])))>.95)   INSTAB= 1;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N    = 256;MAR  =   2;LAR  =   2;MMA  =   1;LMA  = LAR;re_im= 'r';mo_no= 'n';tfpm_file_gen;%Bml(LMA+1, 1)= 2;%-------------beta = 1/2;alpha= 1/2;KC=  2;LC= KC*(LAR+LMA);MC= KC*(MAR+MMA);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);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]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)param_norm(A, Aml, 1)param_norm(B, Bml, 1)param_plot('exp', N, 1, 'r', 'g', A, 1, 0, 0, 1, 0, 0, B, 1, 0, 0, 2, 0, 0);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%