function [Bnm, B, REG, PMIN]= tvarma_est_cepsb(Aalpha, MMA, LMA, Aml)% function [Bnm, B, REG, PMIN]= tvarma_est_cepsb(Aalpha, MMA, LMA, Aml)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Estimates TVARMA models from Ralpha. size(Ralpha)= [N, N]. % Uses the complex TV cepstral recursion in IEEESPL vol9(4), apr02,% pp130. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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= 'n';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;%Psi= ones(N);Aalpha= tfarma_ambi(Aml, Bml, N, alpha);Bnmref=  param_expand(Bml, N);n=0m= 1mprime= 0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ralpha= ml_to_nm(Aalpha);% DimensionsN= size(Ralpha, 1);Anm=  param_expand(Aml, N);Bnm= [];for n= 0:N-1   [r_hat, REG, PMIN]= ti_ceps2(Ralpha(n+1, :).', 1e-4);%% Recursion in m% Initialization (m=0)%   BB= [exp(.5*r_hat(N/2+1)) zeros(1, MMA)];  BB= 1;  AA= Anm(n+1, :);% m<0   for m= 1:MMA      bb= m*(r_hat(N/2+m+1) + Anm(n+1, m+1));      if(m>1)         for mprime= 1:m-1            bb= bb + mprime*( AA(mprime+1)*BB(m-mprime+1) - BB(mprime+1)*AA(m-mprime+1) );         end;         bbb= 0;         for mprime= 1:m-1            for mpp= 0:m-mprime               bbb= bbb + AA(mpp+1)*BB(m-mprime-mpp+1);	             end;            bb= bb + mprime*r_hat(N/2+mprime+1)*bbb;         end;      end;%if(m>1)      BB(m+1)= bb/m;   end;   Bnm= [Bnm; BB];end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%max(max(abs(Bnm-Bnmref)))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%HMAmonic= [zeros(N, N/2) Bnm zeros(N, N/2-MMA-1)];HAR     = [zeros(N, N/2) Anm zeros(N, N/2-MMA-1)];Wmonic= abs(nm_to_nk(HMAmonic)./nm_to_nk(HAR)).^2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(1);mesh(abs(Wmonic))figure(2);mesh(abs(tfarma_wvsp(Aml, Bml, N, 1/2)))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%B0= sqrt(real(mean( (nm_to_nk(Ralpha) ./ Wmonic).' ).'));%figure(1);clf;plot(real(B0));%figure(2);clf;plot(real(param_expand(Bml, N)))Bnm= Bnm.*(B0*ones(1, MMA+1));B= param_compress(Bnm, LMA);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%max(max(abs(Bnm-Bnmref)))max(max(abs(Bml-B)))param_norm(param_expand(Bml, N), Bnm, 1)param_norm(Bml, B, 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%