motor.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate matrix X containing 100 selected realizations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
load('~/matlab/data/Motordaten/bmw1000.mat');
MM= 100;
N0= 1468;
N= 256;
offset= 64;
alpha= 1/2;
mm0= 59;

C= c1(offset:2:offset+2*N-1, :);
[PsiRD, Mask, v2]= tf_multiwin(N, 50, 10, 0, 2, 1);

NOTT= [13 38 51 61 73 82 83 87 92 102 107 110 114 116 118 122 126 128 138 140 144 146 148 150 151 158 160 161 162 163 166 172 176 177 178 179 181 183 4 5 12 18 21 23 24 25 29 35 37 39 40 42 43 44 52 53 55 56 63 64 71 75 85 94 96 97 99 105 108 109 117 124 133 137 143 147 157 159 164 165 173 174 180];

[B,A] = butter(5, 12/N, 'high');
X= [];
for mm= 1:183
   if(~sum(mm==NOTT))
      mm
      x= C(:, mm);
      x= x-mean(x);
%      figure(1);plot(abs(fft(x).^2))
      x= filter(B, A, x);
%      figure(2);plot(abs(fft(x).^2))
%      figure(3);plot(x)
%      figure(4);tf_show(nm_to_nk(ml_to_nm(nm_to_ml(ker_to_lag(x*x', -1, alpha)).*conj(Psi))))
      drawnow
      X= [X x];
%      pause
   end;
end;
X= X./max(max(X));
MM= size(X, 2)

figure(1);clf;
subplot(2, 1, 1);cla
plot(X(:, mm0), 'Linewidth', 1, 'Linestyle', '-', 'Color', 'k')
axis([1 N, -1.1 1.1])
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])

subplot(2, 1, 2);cla
S= mean(abs(fft(X)).^2.');
S= [S(N/2+1:N) S(1:N/2)];
plot(S, 'Linewidth', 1, 'Linestyle', '-', 'Color', 'k')
axis([1 N, 0 15])
set(gca, 'XTick', [1 N/2+1 N])
set(gca, 'XTickLabel', [-128 0 127])
print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorx.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The true (but a little bit smoothed) RS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A= nm_to_ml(ker_to_lag(X*X', -1, alpha));
%figure(99);tf_show(A)
%figure(98);tf_show(Psi)

Rengine= real(nm_to_nk(ml_to_nm(A.*conj(PsiRD))));

figure(2);tf_show(rot90(Rengine(:, 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/2 ])
set(gca, 'YTickLabel', [127 0])
print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRx.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The RD of x_mm0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ARD= nm_to_ml(ker_to_lag(X(:, mm0)*X(:, mm0)', -1, alpha)).*conj(PsiRD);
%figure(97);tf_show(ARD)
RRD= real(nm_to_nk(ml_to_nm(ARD)));

figure(3);tf_show(rot90(RRD(:, 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/2 ])
set(gca, 'YTickLabel', [127 0])
print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRxi.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2nd, 3rd, 4th moments over time n
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VAn= [];
SKn= [];
KUn= [];
for n= 0:N-1
%   hist(X(n+1, :), 10);
   VAn= [VAn var(X(n+1, :))];
   SKn= [SKn skewness(X(n+1, :))];
   KUn= [KUn kurtosis(X(n+1, :))];
%   pause
end;
figure(97);plot(VAn)
figure(98);plot(SKn)
figure(99);plot(KUn)

%----------------------------------------------------------------
if(0)
%----------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find optimum kernel by MSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 30;
Lmax= 30;
MSE= zeros(Mmax, Lmax);

for M= 1:5:Mmax
   for L= 1:5:Lmax
      Psi= tf_multiwin(N, M, L, 0, 2, 1);
      for mm= 1:30
	 mm
	 R= X(:, mm)*X(:, mm)';
	 A= nm_to_ml(ker_to_lag(R, -1, alpha)).*conj(Psi);
	 Rmm= real(nm_to_nk(ml_to_nm(A)));
         MSE(M, L)= MSE(M, L) + norm(Rmm-Rengine)/norm(Rengine)
	 imagesc(MSE(1:5:end, 1:5:end));
	 drawnow
      end;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFAR model across ensemble by AIC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 5;
Lmax= 8;
mm= 1;
Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);
ARAIChist= zeros(Mmax, Lmax);
for mm= 1:20
   mm
   [AmlAIC, BmlAIC]= tfar_fit_aic(X(:, mm), Mmax, Lmax, Psi);
   [MAIC, LAIC]= param_dim(AmlAIC);
   ARAIChist(MAIC, LAIC)= ARAIChist(MAIC, LAIC)+1;
   figure(99);clf;imagesc(ARAIChist);drawnow
end
ARAIChist
[MARAIC, LARAIC]=find(ARAIChist==max(max(ARAIChist)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFMA model across ensemble by AIC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 10;
Lmax= 10;
mm= 1;
AIChist= zeros(Mmax, Lmax);
Psi= tf_multiwin(N, 3*Mmax, 4*Lmax, 0, 2, 1);
for mm= 1:20
   mm
   BmlAIC= tfma_fit_aic(X(:, mm), Mmax, Lmax, Psi);
   [MAIC, LAIC]= param_dim(BmlAIC);
   AIChist(MAIC, LAIC)= AIChist(MAIC, LAIC)+1;
   figure(99);clf;imagesc(AIChist);drawnow
end
AIChist
[MMAAIC, LMAAIC]=find(AIChist==max(max(AIChist)))
%----------------------------------------------------------------
end;
%----------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The RD, TFMA, TFAR, TFARMA of x_mm0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MMA_MA= 12;
LMA_MA= 4;

MAR_AR= 8;
LAR_AR= 4;

MAR_ARMA= 4;
LAR_ARMA= 4;
MMA_ARMA= MAR_ARMA-1;
LMA_ARMA= LAR_ARMA;

PsiMA  = tf_multiwin(N, 3*MMA_MA  , 4*LMA_MA  , 0, 2, 1);
PsiAR  = tf_multiwin(N, 3*MAR_AR  , 4*LAR_AR  , 0, 2, 1);
PsiARMA= tf_multiwin(N, 3*MAR_ARMA, 4*LAR_ARMA, 0, 2, 1);

maxMA     = [];
maxAR     = [];
maxARMA_AR= [];
maxARMA_MA= [];
%for mm0= 1:100
   R= X(:, mm0)*X(:, mm0)';

   ARD  = nm_to_ml(ker_to_lag(R, -1, alpha)).*conj(PsiRD  );
   RRD  = real(nm_to_nk(ml_to_nm(ARD)));

   AMA  = nm_to_ml(ker_to_lag(R, -1, alpha)).*conj(PsiMA  );
   BmlMA= tfma_est_cepsb(AMA, MMA_MA, LMA_MA);
   TDIRMA= param_tdir(param_expand(BmlMA, N));
   RMA  = tfarma_wvsp(    1  , BmlMA  , N, alpha);
   maxMA= [maxMA; max(max(abs(TDIRMA)))]
   figure(4);tf_show(rot90(RMA(:, 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/2 ])
   set(gca, 'YTickLabel', [127 0])
   print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRMA.eps

   AAR  = nm_to_ml(ker_to_lag(R, -1, alpha)).*conj(PsiAR  );
   [AA, BB]= tfar_est_tfywu(AAR(N/2-3*LAR_AR+1:N/2+3*LAR_AR+1, N/2-MAR_AR+1:N/2+MAR_AR+1), N);
   AmlAR= AA(:, :, MAR_AR);
   BmlAR= BB(:, :, MAR_AR);
   TDIRAR= param_tdir(param_expand(AmlAR, N));
   RAR  = tfarma_wvsp(AmlAR  , BmlAR  , N, alpha);
   maxAR= [maxAR; max(max(abs(TDIRAR)))]
   figure(5);tf_show(rot90(RAR(:, 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/2 ])
   set(gca, 'YTickLabel', [127 0])
   print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRAR.eps

   AARMA= nm_to_ml(ker_to_lag(R, -1, alpha)).*conj(PsiARMA);
   AA= tfarma_est_tfywu(AARMA(N/2-3*LAR_ARMA+1:N/2+3*LAR_ARMA+1, N/2+MMA_ARMA-MAR_ARMA+1:N/2+MMA_ARMA+MAR_ARMA+1), MMA_ARMA, N);
   AmlARMA= AA(:, :, MAR_ARMA);
   BmlARMA= tfarma_est_cepsb(AARMA, AmlARMA, MMA_ARMA, LMA_ARMA);
   TDIRAR= param_tdir(param_expand(AmlARMA, N));
   TDIRMA= param_tdir(param_expand(BmlARMA, N));
   RARMA= tfarma_wvsp(AmlARMA, BmlARMA, N, alpha);
   maxARMA_AR= [maxARMA_AR; max(max(abs(TDIRAR)))]
   maxARMA_MA= [maxARMA_MA; max(max(abs(TDIRMA)))]
   figure(6);tf_show(rot90(RARMA(:, 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/2 ])
   set(gca, 'YTickLabel', [127 0])
   print -depsc2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRARMA.eps
   
   drawnow;
%end;













%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check all 3 models incl RD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Psi= tf_multiwin(N, MAR, LAR, 0, 2, 1);
SS= zeros(N);
PP= zeros(N);
for mm= 1:MM
   mm
   Qalpha= corr_est(X(:, mm), X(:, mm), -1, alpha);
   Aalpha= fft(Qalpha);
   Aalpha= [ Aalpha(N/2+1:N, :); Aalpha(1:N/2, :)].*conj(Psi);
   Ayy= Aalpha(N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);
   [AA, BB]= tfar_est_tfywu(Ayy, N);
   P= tfarma_wvsp(AA(:, :, end), BB(:, :, end), N, alpha);
   figure(1);mesh(P)
   S= real(ml_to_nm(nm_to_nk(Aalpha)));
   figure(2);mesh(S)
   drawnow
   SS= SS+S;
   PP= PP+P;
end
figure(1);mesh(abs(PP))
figure(2);mesh(abs(SS))