tfpm_motor.m

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

EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
PsiMA= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);
for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfma_fitB(sig, Mmax, Lmax, PsiMA, lambda);
   E= log(E);
   MDL= E + (log(N+1)+rho)  *(NA+NB+1/2)/N;
   AIC= E + 2               *(NA+NB    )/N;
   [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
   [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
   if(~LMDL)
      LMDL= 1;
   end;
   if(~LAIC)
      LAIC= 1;
   end;
   [MMDL LMDL MAIC LAIC]
   OMDL= [OMDL  [MMDL; LMDL]];
   OAIC= [OAIC  [MAIC; LAIC]];
   A= fft(corr_est(sig, sig, -1, alpha));
   A= [ A(N/2+1:N, :); A(1:N/2, :)];
   [BMDL, REG, PMIN]= tfma_est_cepsb(A.*conj(PsiMA), MMDL, LMDL);
   [BAIC, REG, PMIN]= tfma_est_cepsb(A.*conj(PsiMA), MAIC, LAIC);
   [BMDL, lambdamax, mmMDL, nrFDIR]= param_stabilize(BMDL, N, lambda);
   [BAIC, lambdamax, mmAIC, nrFDIR]= param_stabilize(BAIC, N, lambda);

   RMDL= tfarma_wvsp(1, BMDL, N, 1/2);
   RAIC= tfarma_wvsp(1, BAIC, N, 1/2);
   RD  = real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));
%   figure(1);clf;tf_show(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;tf_show(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;tf_show(rot90(RD(:, 1:N/2)))
%   figure(1);clf;mesh(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;mesh(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;mesh(rot90(RD(:, 1:N/2)))
%   drawnow
end;
figure(1);clf;
tf_show(rot90(RMDL(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
set(gca, 'YTickLabel', ['R'; 'S'; 'T'])
text(220, 10, 'MA-B-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorMA-B-MDL.eps

figure(2);clf;
tf_show(rot90(RAIC(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
set(gca, 'YTickLabel', ['R'; 'S'; 'T'])
text(220, 10, 'MA-B-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorMA-B-AIC.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFAR model across ensemble by MDL/AIC (A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 18;
Lmax= 3;
NA= (1:Mmax)'*(2*(0:Lmax)+1);
NB= ones(Mmax, 1)*(2*(0:Lmax)+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;
EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
PsiAR= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);

for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfar_fitA(sig, Mmax, Lmax, PsiAR, lambda);
   E= log(E);
   MDL= E + (log(N+1)+rho)  *(NA+NB+1/2)/N;
   AIC= E + 2               *(NA+NB    )/N;
   [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
   [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
   if(~LMDL)
      LMDL= 1;
   end;
   if(~LAIC)
      LAIC= 1;
   end;
   [MMDL LMDL MAIC LAIC]
   OMDL= [OMDL  [MMDL; LMDL]];
   OAIC= [OAIC  [MAIC; LAIC]];
   A= fft(corr_est(sig, sig, -1, alpha));
   A= [ A(N/2+1:N, :); A(1:N/2, :)];
   [AMDL, BMDL, REG, PMIN]= tfar_est_cepsb(A.*conj(PsiAR), MMDL, LMDL);
   [AAIC, BAIC, REG, PMIN]= tfar_est_cepsb(A.*conj(PsiAR), MAIC, LAIC);
   [AMDL, lambdamax, mmMDL, nrFDIR]= param_stabilize(AMDL, N, lambda);
   [AAIC, lambdamax, mmAIC, nrFDIR]= param_stabilize(AAIC, N, lambda);

   RMDL= tfarma_wvsp(AMDL, BMDL, N, 1/2);
   RAIC= tfarma_wvsp(AAIC, BAIC, N, 1/2);
   RD  = real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));
%   figure(1);clf;mesh(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;mesh(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;mesh(rot90(RD(:, 1:N/2)))
%   figure(1);clf;tf_show(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;tf_show(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;tf_show(rot90(RD(:, 1:N/2)))
   drawnow
end;
figure(1);clf;
tf_show(rot90(RMDL(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'AR-A-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorAR-A-MDL.eps

figure(2);clf;
tf_show(rot90(RAIC(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'AR-A-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorAR-A-AIC.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFAR model across ensemble by MDL/AIC (B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 18;
Lmax= 3;
NA= (1:Mmax)'*(2*(0:Lmax)+1);
NB= ones(Mmax, 1)*(2*(0:Lmax)+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;
EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
PsiAR= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);

for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfar_fitB(sig, Mmax, Lmax, PsiAR, lambda);
   E= log(E);
   MDL= E + (log(N+1)+rho)  *(NA+NB+1/2)/N;
   AIC= E + 2               *(NA+NB    )/N;
   [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
   [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
   if(~LMDL)
      LMDL= 1;
   end;
   if(~LAIC)
      LAIC= 1;
   end;
   [MMDL LMDL MAIC LAIC]
   OMDL= [OMDL  [MMDL; LMDL]];
   OAIC= [OAIC  [MAIC; LAIC]];
   A= fft(corr_est(sig, sig, -1, alpha));
   A= [ A(N/2+1:N, :); A(1:N/2, :)];
   [AMDL, BMDL, REG, PMIN]= tfar_est_cepsb(A.*conj(PsiAR), MMDL, LMDL);
   [AAIC, BAIC, REG, PMIN]= tfar_est_cepsb(A.*conj(PsiAR), MAIC, LAIC);
   [AMDL, lambdamax, mmMDL, nrFDIR]= param_stabilize(AMDL, N, lambda);
   [AAIC, lambdamax, mmAIC, nrFDIR]= param_stabilize(AAIC, N, lambda);

   RMDL= tfarma_wvsp(AMDL, BMDL, N, 1/2);
   RAIC= tfarma_wvsp(AAIC, BAIC, N, 1/2);
   RD  = real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));
%   figure(1);clf;mesh(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;mesh(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;mesh(rot90(RD(:, 1:N/2)))
%   figure(1);clf;tf_show(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;tf_show(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;tf_show(rot90(RD(:, 1:N/2)))
%   drawnow
end;
figure(1);clf;
tf_show(rot90(RMDL(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'AR-B-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorAR-B-MDL.eps

figure(2);clf;
tf_show(rot90(RAIC(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'AR-B-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorAR-B-AIC.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFARMA model across ensemble by MDL/AIC (A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 12;
Lmax= 3;
NA= (1:Mmax)'*(2*(0:Lmax)+1);
NB= (1:Mmax)'*(2*(0:Lmax)+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;
EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
[PsiAA, mask, v2]= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);
for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfarma_fitA(sig, Mmax, Lmax, PsiAA, lambda);
   E= log(E);
   MDL= E + (log(N+1)+rho)  *(NA+NB+1/2)/N;
   AIC= E + 2               *(NA+NB    )/N;
   [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
   [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
   if(~LMDL)
      LMDL= 1;
   end;
   if(~LAIC)
      LAIC= 1;
   end;
   [MMDL LMDL MAIC LAIC]
   OMDL= [OMDL  [MMDL; LMDL]];
   OAIC= [OAIC  [MAIC; LAIC]];
   A= fft(corr_est(sig, sig, -1, alpha));
   A= [ A(N/2+1:N, :); A(1:N/2, :)];

   AAMDL= tfarma_est_tfywu(A(N/2-3*LMDL+1:N/2+3*LMDL+1, N/2:N/2+2*MMDL), MMDL-1, N);
   [AMDL, lambdamax, mm1, nrFDIR]= param_stabilize(AAMDL(:, :, end), N, lambda);
   [BMDL, REG, PMIN]= tfarma_est_cepsb(A.*conj(PsiAA), AMDL, MMDL-1, LMDL);
   [BMDL, lambdamax, mm2, nrFDIR]= param_stabilize(BMDL, N, lambda);

   AAAIC= tfarma_est_tfywu(A(N/2-3*LAIC+1:N/2+3*LAIC+1, N/2:N/2+2*MAIC), MAIC-1, N);
   [AAIC, lambdamax, mm1, nrFDIR]= param_stabilize(AAAIC(:, :, end), N, lambda);
   [BAIC, REG, PMIN]= tfarma_est_cepsb(A.*conj(PsiAA), AAIC, MAIC-1, LAIC);
   [BAIC, lambdamax, mm2, nrFDIR]= param_stabilize(BAIC, N, lambda);

   RMDL= tfarma_wvsp(AMDL, BMDL, N, 1/2);
   RAIC= tfarma_wvsp(AAIC, BAIC, N, 1/2);
   RD  = real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));
%   figure(1);clf;tf_show(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;tf_show(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;tf_show(rot90(RD(:, 1:N/2)))
%   figure(1);clf;mesh(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;mesh(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;mesh(rot90(RD(:, 1:N/2)))
%   drawnow
end;
figure(1);clf;
tf_show(rot90(RMDL(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'ARMA-A-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorARMA-A-MDL.eps

figure(2);clf;
tf_show(rot90(RAIC(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'ARMA-A-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorARMA-A-AIC.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFARMA model across ensemble by MDL/AIC (B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 12;
Lmax= 3;
NA= (1:Mmax)'*(2*(0:Lmax)+1);
NB= (1:Mmax)'*(2*(0:Lmax)+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;
EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
[PsiAA, mask, v2]= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);
for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfarma_fitB(sig, Mmax, Lmax, PsiAA, lambda);
   E= log(E);
   MDL= E + (log(N+1)+rho)  *(NA+NB+1/2)/N;
   AIC= E + 2               *(NA+NB    )/N;
   [MMDL, LMDL] = find(MDL==min(min(MDL)));LMDL= LMDL-1;
   [MAIC, LAIC] = find(AIC==min(min(AIC)));LAIC= LAIC-1;
   if(~LMDL)
      LMDL= 1;
   end;
   if(~LAIC)
      LAIC= 1;
   end;
   [MMDL LMDL MAIC LAIC]
   OMDL= [OMDL  [MMDL; LMDL]];
   OAIC= [OAIC  [MAIC; LAIC]];
   A= fft(corr_est(sig, sig, -1, alpha));
   A= [ A(N/2+1:N, :); A(1:N/2, :)];

   AAMDL= tfarma_est_tfywu(A(N/2-3*LMDL+1:N/2+3*LMDL+1, N/2:N/2+2*MMDL), MMDL-1, N);
   [AMDL, lambdamax, mm1, nrFDIR]= param_stabilize(AAMDL(:, :, end), N, lambda);
   [BMDL, REG, PMIN]= tfarma_est_cepsb(A.*conj(PsiAA), AMDL, MMDL-1, LMDL);
   [BMDL, lambdamax, mm2, nrFDIR]= param_stabilize(BMDL, N, lambda);

   AAAIC= tfarma_est_tfywu(A(N/2-3*LAIC+1:N/2+3*LAIC+1, N/2:N/2+2*MAIC), MAIC-1, N);
   [AAIC, lambdamax, mm1, nrFDIR]= param_stabilize(AAAIC(:, :, end), N, lambda);
   [BAIC, REG, PMIN]= tfarma_est_cepsb(A.*conj(PsiAA), AAIC, MAIC-1, LAIC);
   [BAIC, lambdamax, mm2, nrFDIR]= param_stabilize(BAIC, N, lambda);

   RMDL= tfarma_wvsp(AMDL, BMDL, N, 1/2);
   RAIC= tfarma_wvsp(AAIC, BAIC, N, 1/2);
   RD  = real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));
%   figure(1);clf;tf_show(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;tf_show(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;tf_show(rot90(RD(:, 1:N/2)))
%   figure(1);clf;mesh(rot90(RMDL(:, 1:N/2)))
%   figure(2);clf;mesh(rot90(RAIC(:, 1:N/2)))
%   figure(3);clf;mesh(rot90(RD(:, 1:N/2)))
%   drawnow
end;
figure(1);clf;
tf_show(rot90(RMDL(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'ARMA-B-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorARMA-B-MDL.eps

figure(2);clf;
tf_show(rot90(RAIC(:, 1:N/2)))
colormap(flipud(gray))
axis([1 N, 1 N/2])
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 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'ARMA-B-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorARMA-B-AIC.eps