clear;tfpmN= 256;mm0= 9;MM= 100;alpha= 1/2;lambda= .98;%PsiRS= tf_multiwin(N, 50, 10, 0, 2, 1);load motor%load motorlow;R= X*X';RR= max(max(abs(R)))R= R+RR/100*eye(N);%figure(1);clf;mesh(R);A= nm_to_ml(ker_to_lag(R, -1, 1/2));%figure(2);clf;mesh(abs(A));P= ml_to_nk(A);%figure(1);clf;imagesc(rot90(abs(P)));PsiRD= tf_window(N, 0, 0, 20, 80);%PsiRD= tf_multiwin(N, 4, 4, 0, 2, 1);%figure(2);clf;mesh(abs(PsiRD));%figure(3);tf_show(rot90(abs(ml_to_nk(PsiRD))))AW= A.*conj(PsiRD);%figure(4);clf;mesh(abs(AW));PW= ml_to_nk(AW);figure(1);clf;imagesc(abs(rot90(PW(:, N/2+1:N))));colormap(flipud(gray))set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])text(220, 10, 'Porig')%print -deps2 /users/mjachan/tex/prj/TF/figures/fig_Porig.epsRW= real(lag_to_ker(ml_to_nm(AW), 1/2));norm(RW-RW')RW= (RW+RW')/2;%figure(3);clf;mesh(RW);HW= chol(RW)';%figure(4);clf;mesh(HW);HWalpha= ker_to_lag(HW, -1, 1/2);%figure(6);clf;mesh(rot90(HWalpha(:, N/2+1:N/2+64)));colormap(flipud(gray));axis([1 N 1 N/4 -.5 .9])figure(2);clf;subplot(2, 1, 1)mesh(log(.01+rot90(abs(HWalpha(:, N/2+1:N/2+64)))));colormap(flipud(gray));axis([1 N 1 N/4 -5 0])set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])set(gca, 'YTick', [1 N/8 N/4])set(gca, 'YTickLabel', ['tap63'; 'tap31'; 'tap00'])text(180, 10, 0, 'Horig')SW= nm_to_ml(HWalpha);subplot(2, 1, 2)mesh(log(.01+rot90(abs(SW(7*N/16+1:9*N/16, N/2+1:N/2+64)))));colormap(flipud(gray));axis([1 N/8 1 N/4 -5 5])%figure(7);clf;mesh(rot90(abs(SW(7*N/16+1:9*N/16, N/2+1:N/2+64))));colormap(flipud(gray).^8)%figure(7);clf;imagesc(rot90(abs(SW(7*N/16+1:9*N/16, N/2+1:N/2+64))));colormap(flipud(gray).^8);axis([1 N/8 1 N/4 0 110])set(gca, 'XTick', [1 N/16 N/8])set(gca, 'XTickLabel', ['-16'; '000'; ' 15'])set(gca, 'YTick', [1 N/8 N/4])set(gca, 'YTickLabel', ['tap63'; 'tap31'; 'tap00'])text(25, 10, 4, 'Sorig')%print -deps2 /users/mjachan/tex/prj/TF/figures/fig_SHorig.epsLW= nm_to_nk(HWalpha);%subplot(3, 1, 3)%imagesc(rot90(abs(LW(:, N/2+1:end))))%mesh(rot90(abs(LW(:, N/2+1:end))))MMIN= 1;MMAX= 30;LMIN= 0;LMAX= 20;Mopt= 13;Lopt= 6;Mopt= 14;Lopt= 4;err1= zeros(MMAX, LMAX+1);err2= zeros(MMAX, LMAX+1);for L= LMIN:LMAX   for M= MMIN:MMAX      [A, B]= tfarma_app_gkmu(SW, M, L, M-1, L, N);      [A1, lambdamax, mm, nrFDIR]= param_stabilize(A, N, .98, 1);      [B1, lambdamax, mm, nrFDIR]= param_stabilize(B, N, .98, 1);                  Halpha= tfarma_impr(A1, B1, N, 1/2);figure(2);clf;subplot(2, 1, 1)mesh(log(.01+rot90(abs(Halpha(:, N/2+1:N/2+64)))));colormap(flipud(gray));axis([1 N 1 N/4 -5 0])set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])set(gca, 'YTick', [1 N/8 N/4])set(gca, 'YTickLabel', ['tap63'; 'tap31'; 'tap00'])text(180, 10, 1, sprintf('H%d, %d-%03d', M, L, mm))            Salpha= tfarma_sprf(A1, B1, N, 1/2);subplot(2, 1, 2)mesh(log(.01+rot90(abs(Salpha(7*N/16+1:9*N/16, N/2+1:N/2+64)))));colormap(flipud(gray));axis([1 N/8 1 N/4 -5 5])%figure(7);clf;mesh(rot90(abs(SW(7*N/16+1:9*N/16, N/2+1:N/2+64))));colormap(flipud(gray).^8)%figure(7);clf;imagesc(rot90(abs(SW(7*N/16+1:9*N/16, N/2+1:N/2+64))));colormap(flipud(gray).^8);axis([1 N/8 1 N/4 0 110])set(gca, 'XTick', [1 N/16 N/8])set(gca, 'XTickLabel', ['-16'; '000'; ' 15'])set(gca, 'YTick', [1 N/8 N/4])set(gca, 'YTickLabel', ['tap63'; 'tap31'; 'tap00'])text(25, 10, 4, 'SML')%print -deps2 /users/mjachan/tex/prj/TF/figures/fig_SHapprox.eps      Lalpha= tfarma_weyl(A1, B1, N, 1/2);      Palpha= abs(Lalpha).^2;      figure(4);clf;imagesc(abs(rot90(Palpha(:, N/2+1:N))));colormap(flipud(gray))set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])text(220, 10, 'Porig')%print -deps2 /users/mjachan/tex/prj/TF/figures/fig_Papprox.eps      beep;      err1(M, L+1)= norm(PW-abs(Lalpha).^2)/norm(PW);      err2(M, L+1)= norm(PW-abs(Lalpha).^2, 'fro')/norm(PW, 'fro');%      pause   end;end;