function Lalpha= vtfarma_weyl(AML, BML, N, B, alpha, tol)% function Lalpha= vtfarma_weyl(AML, BML, N, B, alpha, tol)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Computes the Weyl matrix of the sytem. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmalpha= 1/2;N= 256;J= 0;B= J+1;offset= 200;M= 4;L= 3;filename= sprintf('chpar%02d%04d%02d%02d-%04d', J, N, M, L, offset)load(filename)BML= B0L;tol= 0.01;clear;tfpmalpha= 1/2;N= 128;M0= 3;L0= 2;J0= 1;B= J0;filename= sprintf('vtfar%02d%02d%02d%04d', M0, L0, J0, N)load(filename);BML= B0L;tol= 0.01;%AML= eye(J0);%BML= eye(J0);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Dimensions[J, J, LAR, MAR]= size(AML);LAR= (LAR-1)/2;MAR= MAR-1;[J, J, LMA, MMA]= size(BML);LMA= (LMA-1)/2;MMA= MMA-1;alpha= 1/2;aml= zeros(J, J, N, N);bml= zeros(J, J, N, N);for j= 1:J   for jp= 1:J      if(abs(j-jp)<B)         aml(j, jp, N/2+1-LAR:N/2+1+LAR, N/2+1:N/2+1+MAR)= AML(j, jp, :, :);         bml(j, jp, N/2+1-LMA:N/2+1+LMA, N/2+1:N/2+1+MMA)= BML(j, jp, :, :);      else         aml(j, jp, N/2+1-LAR:N/2+1+LAR, N/2+1:N/2+1+MAR)= zeros(size(AML(j, jp, :, :)));         bml(j, jp, N/2+1-LMA:N/2+1+LMA, N/2+1:N/2+1+MMA)= zeros(size(BML(j, jp, :, :)));      end;   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for j= 1:J   for jp= 1:J      figure(1);tf_show(reshape(aml(j, jp, :, :), N, N))      figure(2);tf_show(reshape(bml(j, jp, :, :), N, N))      pause   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Anm= ifft(fftshift(aml, 3), [], 3);clear amlAnk= N*fftshift(fft(fftshift(Anm, 4), [], 4), 4);clear Anm Bnm= ifft(fftshift(bml, 3), [], 3);clear bmlBnk= N*fftshift(fft(fftshift(Bnm, 4), [], 4), 4);clear Bnm%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for j= 1:J   for jp= 1:J      figure(1);mesh(abs(reshape(Ank(j, jp, :, :), N, N)))      figure(2);mesh(abs(reshape(Bnk(j, jp, :, :), N, N)))%      figure(1);tf_show(reshape(Ank(j, jp, :, :), N, N))%      figure(2);tf_show(reshape(Bnk(j, jp, :, :), N, N))      pause   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Regularized inversionfor n= 1:N   for k= 1:N      [Unk, Snk, Vnk]= svd(Ank(:, :, n, k));%      Snk      Snk= Snk+tol*Snk(1, 1)*eye(J);%      Lalpha(:, :, n, k)= pinv(Ank(:, :, n, k), 1e-2)*Bnk(:, :, n, k);%      Lalpha(:, :, n, k)= inv(Ank(:, :, n, k))*Bnk(:, :, n, k);%      Lalpha(:, :, n, k)= Unk*inv(Snk)*Vnk'*Bnk(:, :, n, k);      Lalpha(:, :, n, k)= inv(Vnk)'*inv(Snk)*inv(Unk)*Bnk(:, :, n, k);   end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for j= 1:J   for jp= 1:J      figure(2);mesh(abs(reshape(Lalpha(j, jp, :, :), N, N)))      pause   end;end;Lalpha0= tfarma_weyl(reshape(AML, 2*LAR+1, MAR+1), reshape(BML, 2*LMA+1, MMA+1), N, alpha);figure(1);cesh(Lalpha0)figure(2);cesh(reshape(Lalpha, N, N))figure(3);cesh(Lalpha0-reshape(Lalpha, N, N))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%