tfpm_cornel.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
load cornel.mat;
N= 512;
alpha= 1/2;

x1= roi1(20:end).';
x1= sqrt(N)*x1/sqrt(energy(x1));
v1= (randn(N, 1)+j*randn(N, 1))/4;
10*log10(energy(x1)/energy(v1))
x1= x1+v1;
figure(1);clot(x1)

x2= roi2(20:end).';
x2= sqrt(N)*x2/sqrt(energy(x2));
v2= (randn(N, 1)+j*randn(N, 1))/4;
10*log10(energy(x2)/energy(v2))
x2= x2+v2;
figure(2);clot(x2)

[Aml1, Bml1, EstInfo1]= tfarma20(x1, 10, 5, 'MDL', 2);
[Aml2, Bml2, EstInfo2]= tfarma20(x2, 10, 5, 'MDL', 2);
save cornel_MDL_2.mat Aml1 Bml1 EstInfo1 Aml2 Bml2 EstInfo2

[Aml1, Bml1, EstInfo1]= tfarma20(x1, 10, 5, 'MDL', 3);
[Aml2, Bml2, EstInfo2]= tfarma20(x2, 10, 5, 'MDL', 3);
save cornel_MDL_3.mat Aml1 Bml1 EstInfo1 Aml2 Bml2 EstInfo2

[Aml1, Bml1, EstInfo1]= tfarma20(x1, 10, 5, 'AIC', 2);
[Aml2, Bml2, EstInfo2]= tfarma20(x2, 10, 5, 'AIC', 2);
save cornel_AIC_2.mat Aml1 Bml1 EstInfo1 Aml2 Bml2 EstInfo2

[Aml1, Bml1, EstInfo1]= tfarma20(x1, 10, 5, 'AIC', 3);
[Aml2, Bml2, EstInfo2]= tfarma20(x2, 10, 5, 'AIC', 3);
save cornel_AIC_3.mat Aml1 Bml1 EstInfo1 Aml2 Bml2 EstInfo2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
N= 512;
NFFT= 64;
alpha= 1/2;
load cornel.mat;
x1= roi1(20:end).';
x1= sqrt(N)*x1/sqrt(energy(x1));
x2= roi2(20:end).';
x2= sqrt(N)*x2/sqrt(energy(x2));
load cornel1.mat
load cornel2.mat

% PLOT tfarma20
P1= tfarma_wvsp(Aml1, Bml1, N, alpha);
figure(1);clf;
subplot(2, 1, 1);tf_show(rot90(flipud(P1), 3))
subplot(2, 1, 2);specgram(x1, NFFT, [], symmpad(hamming(64), NFFT), NFFT-1);

P2= tfarma_wvsp(Aml2, Bml2, N, alpha);
figure(2);clf;
subplot(2, 1, 1);tf_show(rot90(flipud(P2), 3))
subplot(2, 1, 2);specgram(x2, NFFT, [], symmpad(hamming(64), NFFT), NFFT-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% PLOT A SIMULATED SINGLE REALIZATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
NFFT= 256;
alpha= 1/2;
N= 1024;
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)];
x= sqrt(N)*x/sqrt(energy(x));
x= x+sqrt(energy(x)/N)*(randn(N, 1)+j*randn(N, 1))/8;
[AmlEst, BmlEst, EstInfo]= tfarma20(x, 6, 10, 'AIC');
[MA, LA]= param_dim(AmlEst);
[MB, LB]= param_dim(BmlEst);
[MA LA MB LB]
save simulated.mat AmlEst BmlEst% NEVER THROW AWAY!!!!!!!!!!!!!!
save meas01.mat AmlEst BmlEst% NEVER THROW AWAY!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT SIMU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
N= 1024;
alpha= 1/2;
load simulated.mat;



PSimu= tfarma_wvsp(AmlEst, BmlEst, N, alpha);

figure(4);clf;tf_show(rot90(flipud(PSimu(:, N/4+1:3*N/4)), 3))
%axis([1 N, 1 N])
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, 60, 'SimuCh')
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 -3.2 2.2])
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 -3.2 2.2])
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 -3.2 2.2])
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

TDIRnm= param_tdir(AmlEst, N);
figure(6);clf;hold on
plot(TDIRnm(:, 1), 'Linewidth', 1, 'Linestyle', '.', 'Color', 'k')
plot(TDIRnm(:, 2), 'Linewidth', 1, 'Linestyle', '.', 'Color', 'r')
plot(TDIRnm(:, 3), 'Linewidth', 1, 'Linestyle', '.', 'Color', 'b')
plot(TDIRnm(:, 4), 'Linewidth', 1, 'Linestyle', '.', 'Color', 'c')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT MEAS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
N= 1024;
alpha= 1/2;
load meas01.mat;



PMeas= tfarma_wvsp(AmlEst, BmlEst, N, alpha);
figure(4);clf;tf_show(rot90(flipud(PMeas), 3))
colormap(flipud(gray))
axis([0 N-1, 0 N-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, 160, 'MeasCh')
line([256 256], [0 1023], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
line([512 512], [0 1023], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
line([768 768], [0 1023], 'Linewidth', 1, 'Linestyle', '--', 'Color', 'k')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/fig_MeasCh.eps

figure(5);clf;
subplot(1, 3, 1);
plot(log10(PMeas(256, :)), 'k')
axis([0 N-1 -3.2 2.2])
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(PMeas(448, :)), 'k')
axis([0 N-1 -3.2 2.2])
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(PMeas(768, :)), 'k')
axis([0 N-1 -3.2 2.2])
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_MeasChSchnitte.eps