function [AML, B0L]= vtfar_parest2(ambi, M, L, B, N)% function [AML, B0L]= vtfar_parest2(ambi, M, L, B, N)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % SPECIAL VERSION for SPAWC05 paper!!!!!!!!!!!!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmalpha= 1/2;N=  256;M=  2;L=  1;J=  4;B=  2;MMAX= 2*M;LMAX= 2*L;filename= sprintf('vtfar%02d%02d%02d%04d', M, L, J, N)load(filename);AML0= AML;B0L0= B0L;% Full width%AMBI= vtfarma_ambi(AML0, B0L0, N, J, alpha);% or bandedAMBI= vtfarma_ambi(AML0, B0L0, N, B, alpha);ambi= AMBI(:, :, N/2+1-3*LMAX:N/2+1+3*LMAX, N/2+1-MMAX:N/2+1+MMAX);m= 1;mp= 1;l= -L;lp= -L;i= 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for j= 1:J   for jp= 1:J      [j jp]      figure(2);clf;tf_show(reshape(AMBI(j, jp, :, :), N, N))      figure(4);clf;mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))      pause   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[J, J, LMAX, MMAX]= size(ambi);J= J-1;LMAX= (LMAX-1)/6;MMAX= (MMAX-1)/2;alpha= 1/2;NB= number_band(J+1, B);% Fill the TBTB MatrixF= zeros(M*(2*L+1)*NB);for m= 1:M for mp= 1:M  for l= -L:L   for lp= -L:L    Fml= zeros(NB);    j= 1;    for i= 0:J     taumin= max(i-B+1, 0);     taumax= min(i+B-1, J);     [taumin taumax];     Fmli= ambi(taumin+1:taumax+1, taumin+1:taumax+1, 3*LMAX+1+l-lp, MMAX+1+m-mp).';%     [size(Fmli) j j+taumax-taumin taumax-taumin+1]     Fml(j:j + (taumax-taumin), j:j + (taumax-taumin))= Fmli;     j= j + (taumax-taumin)+1;    end;    iidx= (m -1)*(2*L+1)*NB + (L-l )*NB+1:(m -1)*(2*L+1)*NB + (L-l +1)*NB;    jidx= (mp-1)*(2*L+1)*NB + (L-lp)*NB+1:(mp-1)*(2*L+1)*NB + (L-lp+1)*NB;    F(iidx, jidx)= Fml;   end;  end; end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[size(F) rank(F) cond(F)]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fill the RHS vectorf= zeros(M*(2*L+1)*NB, 1);for m= 1:M  for l= -L:L    fml= zeros(NB, 1);    j= 1;    for i= 0:J     taumin= max(i-B+1, 0);     taumax= min(i+B-1, J);     [taumin taumax];     fmli= ambi(i+1, taumin+1:taumax+1, 3*LMAX+1+l, MMAX+1+m).';%     [size(Fmli) j j+taumax-taumin taumax-taumin+1]     fml(j:j + (taumax-taumin))= fmli;     j= j + (taumax-taumin)+1;    end;    iidx= (m -1)*(2*L+1)*NB + (L-l )*NB+1:(m -1)*(2*L+1)*NB + (L-l +1)*NB;    f(iidx)= fml;  end;end;% The fast inversion;-)theta= -inv(F)*f;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%theta AML0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ReshapeAML= zeros(J+1, J+1, 2*L+1, M+1);j= 1;for m= 1:M for l= -L:L  for i= 0:J%     [m l tau taup]     taumin= max(i-B+1, 0);     taumax= min(i+B-1, J);     [taumin taumax];     AML(i+1, taumin+1:taumax+1, L+1-l, m+1)= theta(j:j + (taumax-taumin));     j= j + (taumax-taumin)+1;  end; end;end;for l= 1:L   AML(:, :, l, 1)= zeros(J+1);   AML(:, :, 2*L+2-l, 1)= zeros(J+1);end;AML(:, :, L+1, 1)= eye(J+1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[AML AML0]AML-AML0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%AMBI0= vtfarma_ambi(AML0, eye(J0), N, J0, alpha);AMBIE= vtfarma_ambi(AML , eye(J0), N, J0, alpha);for j= 1:J0   for jp= 1:J0      [j jp]      figure(1);clf;mesh(abs(reshape(AMBI0(j, jp, :, :), N, N)))      figure(2);clf;mesh(abs(reshape(AMBIE(j, jp, :, :), N, N)))      figure(3);clf;mesh(abs(reshape(AMBI0(j, jp, :, :)-AMBIE(j, jp, :, :), N, N)))      pause   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute b_{0, l}BBH= zeros(J+1, J+1, N);for l= 0:2*L   for mp= 0:M      for lp= -L:L         BBH(:, :, l+1)= BBH(:, :, l+1) + AML(:, :, L+lp+1, mp+1)*ambi(:, :, 3*LMAX+1+l-lp, MMAX+1-mp);      end;   end;end;for l= -2*L:-1   for mp= 0:M      for lp= -L:L         BBH(:, :, N+1+l)= BBH(:, :, N+1+l) + AML(:, :, L+lp+1, mp+1)*ambi(:, :, 3*LMAX+1+l-lp, MMAX+1-mp);      end;   end;end;bbH= ifft(BBH, N, 3);for n= 0:N-1   [V, D]= eig(bbH(:, :, n+1));   bbH(:, :, n+1)= V*pospart(sqrt(D))*V';end;BBH= fft(bbH, N, 3);B0L= zeros(J+1, J+1, 2*L+1);for l= 0:L   B0L(:, :, L+l+1)= BBH(:, :, l+1);end;for l= -L:-1   B0L(:, :, L+l+1)= BBH(:, :, N+1+l);end;B0L= B0L/N;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sum(sum(sum(sum(abs(B0L-B0L0)))))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmalpha= 1/2;N= 256;M= 2;L= 1;J= 2;B= J;MMAX= 2*M;LMAX= 2*L;filename= sprintf('vtfar%02d%02d%02d%04d', M, L, J, N)load(filename);AML0= AML;B0L0= B0L;AAMMLL= zeros(size(AML0));BB00LL= zeros(size(B0L0));AAMMBBII= zeros(J, J, N, N);MM= 100;for mm= 1:MM   mm   E= randn(N, J);   X= vtfarma_gen(E, AML0, B0L0, B, 1/2);   [AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);   [AML, B0L]= vtfar_parest2(ambi, M, L, B, N);   AAMMBBII= AAMMBBII + AMBI;   AAMMLL= AAMMLL + AML;   BB00LL= BB00LL + B0L;end;AAMMLL= AAMMLL/MM;BB00LL= BB00LL/MM;AAMMBBII= AAMMBBII/MM;[norm(AML0(:)-AAMMLL(:))/norm(AML0(:))][norm(AML0(:)-AAMMLL(:))/norm(AML0(:)) norm(B0L0(:)-BB00LL(:))/norm(B0L0(:))]AMBI= vtfarma_ambi(AML0, B0L0, N, alpha);for j= 1:J   for jp= 1:J      [j jp]      norm(reshape(AAMMBBII(j, jp, :, :), N, N)-reshape(AMBI(j, jp, :, :), N, N))/norm(reshape(AMBI(j, jp, :, :), N, N))      figure(1);tf_show(reshape(AAMMBBII(j, jp, :, :), N, N))      figure(2);tf_show(reshape(AMBI(j, jp, :, :), N, N))      figure(3);mesh(abs(reshape(AAMMBBII(j, jp, :, :), N, N)))      figure(4);mesh(abs(reshape(AMBI(j, jp, :, :), N, N)))      pause   end;end;ambi= AMBI(:, :, N/2+1-3*LMAX:N/2+1+3*LMAX, N/2+1-MMAX:N/2+1+MMAX);[AML, B0L]= vtfar_parest(ambi, M, L, N);[norm(AML0(:)-AML(:)) norm(B0L0(:)-B0L(:))][real(AML) AML0]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%