function [PML, Pre, Post, NrIter]= param_2dfact(Cml, N, re_im, mo_no)% function [PML, Pre, Post, NrIter]= param_2dfact(Cml, N, re_im, mo_no)%   This file is part of the TFPM toolbox v1.0 (c)%   michael.jachan@tuwien.ac.at and underlies the GPL.% % Computes the approximate factorization of the system Cml of order% (M, L) into (1, 1) and LTI or LFI systems. Pre/Post: GM-phd,% p43. Without phase factor --> underspread. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpm;N    = 128;MAR  =   2;LAR  =   2;MMA  =   1;LMA  =   1;re_im= 'r';mo_no= 'm';tfpm_file_gen;%-------------alpha= 1/2;beta = 1/2;Cml= param_tconv(Aml, Bml, N, beta);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Dimensions[M, L]= param_dim(Cml);thetaC= param_stack_ml(Cml);PML= zeros(3, 2, min(M, L));M1= min(1, M);L1= min(1, L);M2= M-M1;L2= L-L1;tfeps= .005125;% Step0: InitNrIter= 0;Pre= [];Post= [];% Random initialization:Mask1= 1./(1+(abs(-L1:L1)+1)'*(1:M1+1)).^2;Mask2= 1./(1+(abs(-L2:L2)+1)'*(1:M2+1)).^2;Pml= param_rand(M1, L1, N, 'i', 'n');% Primitive polyEml= param_rand(M2, L2, N, 'i', 'n');% Residual polyPml= Pml.*Mask1;Emk= Eml.*Mask2;if(re_im=='r')   Pml= param_hermite(Pml);   Eml= param_hermite(Eml);end;thetaP= param_stack_ml(Pml);thetaE= param_stack_ml(Eml);Bml= conv2(Pml, Eml);B= param_2dconv_B(Pml, Eml);E= param_2dconv_E(Pml, Eml);% Step1: 1st Iteration   Eml= param_destack_ml(inv(B'*B)*B'*thetaC, M2, L2);   Eml(:, 1)= param_hermite(Eml(:, 1));   E= param_2dconv_E(Pml, Eml);   Pml= param_destack_ml(inv(E'*E)*E'*thetaC, M1, L1);   Pml(:, 1)= param_hermite(Pml(:, 1));   B= param_2dconv_B(Pml, Eml);   Bml= conv2(Pml, Eml);   Error= norm(Cml-Bml)/norm(Cml);%   Error= [Error norm(Cml-Bml)/norm(Cml) ];Error= Error(end-5:end);   NrIter= NrIter + 1;% Correction   corr= Pml(L1+1, 1);   Pml= Pml/corr;   Eml= Eml*corr;   [Error tfeps]% Step2: 2nd Iteration   Eml= param_destack_ml(inv(B'*B)*B'*thetaC, M2, L2);   Eml(:, 1)= param_hermite(Eml(:, 1));   E= param_2dconv_E(Pml, Eml);   Pml= param_destack_ml(inv(E'*E)*E'*thetaC, M1, L1);   Pml(:, 1)= param_hermite(Pml(:, 1));   B= param_2dconv_B(Pml, Eml);   Bml= conv2(Pml, Eml);   Error= norm(Cml-Bml)/norm(Cml);%   Error= [Error norm(Cml-Bml)/norm(Cml) ];Error= Error(end-5:end);   NrIter= NrIter + 1;% Correction   corr= Pml(L1+1, 1);   Pml= Pml/corr;   Eml= Eml*corr;   [Error tfeps]% Iteration:while(Error>tfeps)   Eml= param_destack_ml(inv(B'*B)*B'*thetaC, M2, L2);   Eml(:, 1)= param_hermite(Eml(:, 1));   E= param_2dconv_E(Pml, Eml);   Pml= param_destack_ml(inv(E'*E)*E'*thetaC, M1, L1);   Pml(:, 1)= param_hermite(Pml(:, 1));   B= param_2dconv_B(Pml, Eml);   Bml= conv2(Pml, Eml);   Error= norm(Cml-Bml)/norm(Cml);%   Error= [Error norm(Cml-Bml)/norm(Cml) ];Error= Error(end-5:end);   NrIter= NrIter + 1;% Correction   corr= Pml(L1+1, 1);   Pml= Pml/corr;   Eml= Eml*corr;   figure(98);clf;mesh(abs(Pml))   figure(99);clf;mesh(abs(Eml))   drawnow   [Error tfeps]end;PML(:, :, 1)= Pml;if( (M2==1) & (L2==1) )   PML(:, :, 2)= Eml;   returnelse   [PML(:, :, 2:end), xxx, yyy, nriter]= param_2dfact(Eml, N, re_im, mo_no);   NrIter= NrIter + nriterend;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%beta= 1/2;Cml= param_rand(5, 5, N, 'i', 'n');Cml= conv2(Aml, conv2(Aml, conv2(Aml, Aml)))figure(1);clf;mesh(abs(Cml))[PML, Pre, Post, NrIter]= param_2dfact(Cml, N, re_im, mo_no)Cml1= param_itconv(PML, N, beta);norm(Cml1-Cml)/norm(Cml)figure(1);clf;mesh(abs(Cml))figure(2);clf;mesh(abs(Cml1))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%