tfpm_channel.m

用于模拟时变非平稳的ARMA过程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
mm= 1;
%MEASURED= 1;MM= 64;load updown_h_rusk.mat;texxt= 'Meas';[B,A] = butter(5, .5);N= 1024;
%MEASURED= 1;MM= 64;load updown_h_rusk.mat;texxt= 'MEAS';N= 2048;
MEASURED= 0;MM= 10;texxt= 'Simu';N= 1024;


NFFT= 256;
alpha= 1/2;

PP= zeros(N);
for mm= 1:MM
   mm
   if(MEASURED)
      if(N==1024)
         x= filter(B, A, h_I2(113:end, mm));
         x= x(1:2:end);
      else
         x= h_I2(113:end, mm);
      end
      NORMA= energy(x);
      x= sqrt(N)*x/sqrt(energy(x));  
      x0= x;
   else
      Nh= 960;
      [h, W]= sim_nwssus_3;
      x= [randn((N-Nh)/2, 1)+j*randn((N-Nh)/2, 1); h; randn((N-Nh)/2, 1)+j*randn((N-Nh)/2, 1)];
      NORMA= energy(x);
      x0= sqrt(N)*x/sqrt(energy(x));
      v= (randn(N, 1)+j*randn(N, 1))/4;
      10*log10(energy(x)/energy(v))
      x= x0+v;
   end;

   [AmlEst, BmlEst, EstInfo]= tfarma20(x, 6, 6, 'AIC');
   filename= sprintf('%s%02d', texxt, mm)
   save(filename, 'AmlEst', 'BmlEst', 'EstInfo', 'NORMA', 'x', 'x0', 'N');
   EstInfo{2}
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT SIMU
%psi= real(weylback(fftshift(fft2(fftshift(Mask)))));
%Psiml=(fftshift(ifft2(fftshift(weylsymb(psi2)))))
%W= weylsymb(hb_comp(RR));                                                       
%weylsymb(hb_comp(Rx0))
%A= weylback(M);
%weylsymb(hbcomp(Ruc))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
load Simu03.mat;N= 1024;
NFFT= 512;

[MAREst, LAREst]= param_dim(AmlEst);
[MMAEst, LMAEst]= param_dim(BmlEst);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L1=96; % window length #1
L2=16; % window length #2
WSimu= wd(x, symmpad(hamming(L1), NFFT), hamming(L2));
WSimu= WSimu(:, L2/2:end-L2/2);
figure(2);clf;tf_show(WSimu)
colormap(flipud(gray))
axis([0 N-1, 0 NFFT-1])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTickLabel', ['A'; 'B'; 'C'; 'D'; 'E'])
text(40, 80, sprintf('SimuChWD'))
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/fig_SimuChWD.eps

ASimu= nm_to_ml(corr_est(x, x, -1, alpha));
[Psi, Mask, v2]= tf_multiwin(N, 50, 25, 0, 2, 1);%GOOD!!!
RSimu= nm_to_nk(ml_to_nm(ASimu.*conj(Psi)));
figure(7);clf;tf_show(rot90(flipud(RSimu(:, N/4+1:3*N/4)), 3))
colormap(flipud(gray))
axis([0 N-1, 0 NFFT-1])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTickLabel', ['A'; 'B'; 'C'; 'D'; 'E'])
text(40, 80, sprintf('SimuChRD'))
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/fig_SimuChRD.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PSimu= tfarma_wvsp(AmlEst, BmlEst, N, alpha);
figure(4);clf;tf_show(rot90(flipud(PSimu(:, N/4+1:3*N/4)), 3))
colormap(flipud(gray))
axis([0 N-1, 0 N/2-1])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTickLabel', ['A'; 'B'; 'C'; 'D'; 'E'])
text(40, 80, sprintf('SimuCh (%d, %d; %d, %d)', MAREst, LAREst, MMAEst, LMAEst))
line([256 256], [0 511], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
line([448 448], [0 511], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
line([768 768], [0 511], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/fig_SimuCh.eps

figure(5);clf;
subplot(1, 3, 1);
plot(log10(PSimu(256, :)), 'k')
axis([0 N-1 -1.6 1.6])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', ['-512'; '-256'; '0000'; ' 255'; ' 511'])
set(gca, 'YTick', -4:2)
set(gca, 'YTickLabel', ['F'; 'G'; 'H'; 'I'; 'J'; 'K'; 'L'])
grid
text(25, 1.5, 'schnitt1')

subplot(1, 3, 2);
plot(log10(PSimu(448, :)), 'k')
axis([0 N-1 -1.6 1.6])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', ['-512'; '-256'; '0000'; ' 255'; ' 511'])
set(gca, 'YTick', -4:2)
set(gca, 'YTickLabel', ['F'; 'G'; 'H'; 'I'; 'J'; 'K'; 'L'])
grid
text(25, 1.5, 'schnitt2')

subplot(1, 3, 3);
plot(log10(PSimu(768, :)), 'k')
axis([0 N-1 -1.6 1.6])
set(gca, 'XTick', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'XTickLabel', ['-512'; '-256'; '0000'; ' 255'; ' 511'])
set(gca, 'YTick', -4:2)
set(gca, 'YTickLabel', ['F'; 'G'; 'H'; 'I'; 'J'; 'K'; 'L'])
grid
text(25, 1.5, 'schnitt3')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/fig_SimuChSchnitte.eps





















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% script:
% Fit MEAS AND COMPUTE RS FROM THEM ALL!!!
% 
% TODO:
% lin/log plot
% surf
% offset= 600:1100
% N= 512
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;

load ../../data/METAMORP/updown_h_rusk.mat
figure(99);clf;mesh(abs(h_I2))

% CASE 1
tfpm
MMAX= 4;
LMAX= 3;
N= 256;
J= 5;
B= 2;
tapoffset= 30;
delta= 200;
startpoint= 0;
endpoint= 1800;

% CASE 2
tfpm
MMAX= 2;
LMAX= 2;
N= 256;
J= 7;
B= 3;
tapoffset= 28;
delta= 200;
startpoint= 0;
endpoint= 1800;


% CASE 3
tfpm
MMAX= 2;
LMAX= 2;
N= 512;
J= 4;
B= 2;
tapoffset= 30;
delta= 200;
startpoint= 600;
endpoint= 600;

% CASE 4
tfpm
MMAX= 3;
LMAX= 3;
N= 512;
J= 4;
B= 2;
tapoffset= 30;
delta= 200;
startpoint= 800;
endpoint= 800;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit MEAS AND COMPUTE RS THEM ALL!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for offset= startpoint:delta:endpoint
   figure(offset/100+1);
   X= 1000*h_I2(offset+1:offset+N, tapoffset:tapoffset+J);
%   X= sqrt(N*(J+1))*X/sqrt(energy(X));
   [offset energy(X)]
   mesh(abs(X))
%   for i= 1:J+1
%      subplot(J+1, 1, i)
%      specgram(X(:, i), 64, 1, [], 63)
%   end;
   drawnow
   
   [AMBI, ambi]= vtfar_ambiest(X, MMAX, LMAX);
   clear AMBI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for LAR= 1:LMAX
   for MAR= 1:MMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      [AML, B0L]= vtfar_parest2(ambi, MAR, LAR, B, N);
      filename= sprintf('chpar%02d%02d%04d%02d%02d-%04d', J, B, N, MAR, LAR, offset)
      save(filename, 'AML', 'B0L')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for MAR= 1:MMAX;
end;%for LAR= 0:LMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
clear h_I2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save Weyl
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for offset= startpoint:delta:endpoint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for LAR= 1:LMAX
   for MAR= 1:MMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      filename= sprintf('chpar%02d%02d%04d%02d%02d-%04d', J, B, N, MAR, LAR, offset)
      load(filename)
      BML= B0L;
      Lalpha= vtfarma_weyl(AML, BML, N, B, 1/2, .01);%!!!!!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 0:J
   for jp= 0:J
      figure(1);mesh(abs(reshape(Lalpha(j+1, jp+1, :, :), N, N)))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      filename= sprintf('chweyl%02d%02d%04d%02d%02d-%04d', J, B, N, MAR, LAR, offset)
      save(filename, 'Lalpha');
      clear Lalpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for MAR= 1:MMAX;
end;%for LAR= 0:LMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% CASE 3
tfpm
MMAX= 2;
LMAX= 2;
N= 512;
J= 4;
B= 2;
tapoffset= 30;
delta= 200;
startpoint= 600;
endpoint= 600;

% CASE 4
tfpm
MMAX= 2;
LMAX= 2;
N= 512;
J= 4;
B= 2;
tapoffset= 30;
delta= 200;
startpoint= 800;
endpoint= 800;

k= 0;

n= 0
n= N/8
n= 2*N/8
n= 3*N/8
n= 4*N/8
n= 5*N/8
n= 6*N/8
n= 7*N/8



for offset= startpoint:delta:endpoint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for LAR= 2:2
   for MAR= 2:2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%offset=  200;MAR= 4;LAR= 3;n=  15;
%offset=  800;MAR= 4;LAR= 3;n=  58;
%offset= 1200;MAR= 4;LAR= 3;n=   2;
%offset= 1400;MAR= 4;LAR= 3;n= 117;%!!!!!!!!!!!!

filename= sprintf('chpar%02d%02d%04d%02d%02d-%04d', J, B, N, MAR, LAR, offset)
load(filename);BML= B0L;
filename= sprintf('chweyl%02d%02d%04d%02d%02d-%04d', J, B, N, MAR, LAR, offset)
load(filename);

RS= nonWSSUS_RS(Lalpha);
RRSS= fftshift(ifft(fftshift(RS, 4), [], 4), 4);



rrss= fftshift(fft(RRSS(:, :, :, n+1), N, 3), 3);
clear Lalpha RS RRSS

figure(1);clf
subplot(2, 1, 1)
mesh(log10(.0001+abs(   rrss(:, :, k+N/2+1))))
colormap(flipud(gray.^4))
axis([1 N 1 J+1 -4 -1])

subplot(2, 1, 2)
imagesc(log10(.0001+abs(   rrss(:, :, k+N/2+1))))
axis([1 N 1 J+1])

%figure(1);surf(log10(.0001+abs(   rrss(:, :, k+N/2+1))))

%figure(1);mesh((.0001+abs(   rrss(:, :, k+N/2+1))))
%figure(1);surf((.0001+abs(   rrss(:, :, k+N/2+1))))


%colormap((gray))
%axis([1 N 1 J+1 0 2e-2])


filename= sprintf('%02d%02d-%03d', MAR, LAR, n)
text(180, 1, -.1, filename)
%beep
drawnow
%for r= 1:100000;end;
%beep


filename= sprintf('/users/mjachan/tex/prj/TF/figures/fig_RS%02d%02d%02d%02d%04d%04d-%03d.eps', J, B, MAR, LAR, N, offset, n)
%filename= sprintf('/users/mjachan/tex/prj/TF/figures/fig_RS%02d%02d%02d%02d%04d%04d-%03dlin.eps', J, B, MAR, LAR, N, offset, n)


print('-depsc2',  filename)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for MAR= 1:MMAX;
end;%for LAR= 0:LMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
beep
exit





























%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save inverse Weyl
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;

MMAX= 5;
LMAX= 3;
N= 256;
J= 5;
B= 2;
for offset= 0:200:1800
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for LAR= 1:LMAX
   for MAR= 1:MMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      filename= sprintf('chpar%02d%04d%02d%02d-%04d', J, N, MAR, LAR, offset)
      load(filename)
      BML= B0L;
      Linv= vtfarma_Linv(AML, BML, N, J+1, alpha);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j= 0:J
   for jp= 0:J
      figure(1);mesh(abs(reshape(Linv(j+1, jp+1, :, :), N, N)))
      pause
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      filename= sprintf('chLinv%02d%04d%02d%02d-%04d', J, N, MAR, LAR, offset)
      save(filename, 'Linv');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;%for MAR= 1:MMAX;
end;%for LAR= 0:LMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for l= 1:2*L+1
   for m= 1:M+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AML(:, :, l, m)= eye(J+1).*AML(:, :, l, m);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for l= 1:2*L+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BML(:, :, l, 1)= eye(J+1).*BML(:, :, l, 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%