function [B, REG, PMIN]= tfarma_est_cepsb(Aalpha, Aml, MMA, LMA)% function [B, REG, PMIN]= tfarma_est_cepsb(Aalpha, Aml, MMA, LMA)%   This file is part of the TFPM toolbox v0.5 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Estimates TFMA(MMA; LMA; N; 1/2) models 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;N    = 256;MAR  =   3;LAR  =   2;MMA  =   2;LMA  =   2;re_im= 'i';mo_no= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta= 1/2;Aalpha= tfarma_ambi(Aml, Bml, N, alpha);tf_show(Aalpha)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DimensionsN= size(Aalpha, 1);[MAR, LAR]= param_dim(Aml);[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 (m=0)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;B= [[b0l(N-LMA+1:N); b0l(1:LMA+1)] zeros(2*LMA+1, MMA)];for m= 1:MMA   for l= -LMA:LMA      bb= 0;      for mprime= 0:m-1         for lprime= -LAR:LAR            d= 0;	    for mpp= 0:mprime               for lpp= -LAR:LAR                  d= d+param_get(Aml, mpp, lpp)*param_get(B, mprime-mpp, lprime-lpp);	       end;	    end;            bb= bb + (m-mprime)/m*d*Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;	 end;      end;      for mprime= 1:m-1         for lprime= -LAR:LAR            bb= bb + (m-mprime)/m*(param_get(Aml, m-mprime, l-lprime)*param_get(B, mprime, lprime) - ...                param_get(Aml, mprime, lprime)*param_get(B, m-mprime, l-lprime));	 end;      end;      for lprime= -LAR:LAR         bb= bb+param_get(Aml, m, lprime)*param_get(B, 0, l-lprime);      end;      B(LMA+1+l, m+1)= bb;   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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%