function [B, sigma, REG, PMIN]= tfma_est_mepsb(Aalpha, MMA, LMA)% function [B, sigma, REG, PMIN]= tfma_est_mepsb(Aalpha, MMA, LMA)%   This file is part of the TFPM toolbox v0.5 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Estimates monic TFMA(MMA; LMA; N; 1/2) models and B0 from Aalpha.% size(Aalpha)= [N, N], Aalpha(N/2+1, N/2+1)= A_y[0, 0]. % Uses a cepstral recursion based on IEEESPL vol9(4), apr02, pp130. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;MAR  =   0;LAR  =   0;MMA  =   3;LMA  =   2;N    =  64;re_im= 'r';mo_no= 'm';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;%Bml(1:LMA, 1)= zeros(LMA, 1);%Bml(LMA+2:end, 1)= zeros(LMA, 1);sigma0= 1.4;Aalpha= tfarma_ambi(Aml, sigma0*Bml, N, alpha);Psi= tf_multiwin(N, MMA, 2*LMA, 5, min(3, N/32), 2);y= tfarma_gen(randn(N, 1), Aml, Bml, beta);aalpha= ambi_est_cyc(y);Aalpha= aalpha.*conj(Psi);figure(97);cesh(Aalpha)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DimensionsN= size(Aalpha, 1);[Aalpha_hat, REG, PMIN]= tf_ceps2(Aalpha, 1e-4);Aalpha_hat= Aalpha_hat/N;%<===========================%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(98);cesh(Aalpha_hat)REGPMIN%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recursion in m% Initialization (B0)A0= [Aalpha_hat(N/2+1:N, N/2+1); Aalpha_hat(1:N/2, N/2+1)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(97);clot(A0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%rx_hat= real(ifft(A0))*N;b0n= real(exp(rx_hat/2));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(97);clf;hold on   plot(b0n, 'r')   plot(real(Bmn(:, 1)), 'g')hold off;drawnow;RX_HAT= 2*log(real(Bmn(:, 1)));figure(98);clf;hold on   plot(rx_hat, 'r')   plot(RX_HAT, 'g')hold off;drawnow;figure(97);clf;   subplot(2, 1, 1);hold on;      plot(real(A0), 'r');      plot(real(fft(RX_HAT/N)), 'g');   subplot(2, 1, 2);hold on;      plot(imag(A0), 'r');      plot(imag(fft(RX_HAT/N)), 'g');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%b0l= fft(b0n)/N;sigma= b0l(1);B= [[zeros(LMA, 1); 1; zeros(LMA, 1)] zeros(2*LMA+1, MMA)];% m<0for m= 1:MMA   for l= -LMA:LMA      bb= 0;      for mprime= 0:m-1         for lprime= -LMA:LMA            bb= bb + (m-mprime)*param_get(B, mprime, lprime)*Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;	 end;      end;      B(LMA+1+l, m+1)= bb/m;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%param_norm(B, Bml, 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recursion in l% Initialization (l=0)A0= [Aalpha_hat(N/2+1, N/2+1:N) Aalpha_hat(N/2+1, 1:N/2)];rx_hat= real(ifft(A0))*N;bm0= real(exp(rx_hat/2));bk0= fft(bm0)/N;B1= [zeros(LMA, MMA+1); bk0(1:MMA+1); zeros(LMA, MMA+1)];for m= 0:MMA   for l= -LMA:-1      bb= 0;      for mprime= 0:m-1         for lprime= -LMA:LMA            bb= bb + (l-lprime)*param_get(B1, mprime, lprime)*...		Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;	 end;      end;      B1(LMA+1+l, m+1)= bb/l;   end;   for l= 1:LMA      bb= 0;      for mprime= 0:m-1         for lprime= -LMA:LMA            bb= bb + (l-lprime)*param_get(B1, mprime, lprime)*...		Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;	 end;      end;      B1(LMA+1+l, m+1)= bb/l;   end;end;param_norm(B, Bml, 1)param_norm(B1, Bml, 1)param_norm(B1, B, 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%