tfpm_motor.m

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

clear;tfpm
N= 256;
mm0= 9;
MM= 100;
alpha= 1/2;
lambda= .98;
%PsiRS= tf_multiwin(N, 50, 10, 0, 2, 1);
%PsiRD= tf_multiwin(N, 30, 10, 0, 2, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% FOUR-COMPONENT SIGNAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('~/matlab/data/Motordaten/bmw1000.mat');
N0= 1468;
offset= 64;
C= c1(offset:2:offset+2*N-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];

[Bh,Ah] = butter(5, 12/N, 'high');
%[Bh,Ah] = butter(5, 22/N, 'high');
[Bl,Al] = butter(5, 120/N);
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(Bh, Ah, x);
%      x= filter(Bl, Al, x);
      x= x/sqrt(energy(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(PsiRD))))
      drawnow
      X= [X x];
%      pause
   end;
end;
%X= 10*X./max(max(X));
MM= size(X, 2)
EX= energy(X)/MM/N

%%Verrauschen der Daten
SNR= 3;%dB
n= sqrt(EX/10^(SNR/10))*randn(size(X));
energy(X)/energy(n)
Y= X+n;

save('motor.mat', 'X', 'Y', 'PsiRD', 'PsiRS');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TWO-COMPONENT SIGNAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('~/matlab/data/Motordaten/bmw1000.mat');
N0= 1468;
offset= 64;
C= c1(offset:2:offset+2*N-1, :);

WELL= [1 2 9 10 11 12 14 15 16 17 19 21 23 24 27 28 29 30 31 33 ...
      34 35 36 37 41 45 46 47 48 49 50 53 56 58 59 60 ...
       61 63 65 66 70 71 74 77 80 81 82 83 84 ...
       85 87 88 89 91 93 94 97 98 100 101 102 103 104 106 ...
       112 114 118 119 120 121 124 130 131 133 134 135 136 ...
       138 139 141 142 143 149 153 154 155 156 162 163 168 169 170 ...
       171 172 173 174 175 179 180 183];

[Bh,Ah] = butter(5, 22/N, 'high');
[Bl,Al] = butter(5, 120/N);
X= [];
for mm= 1:183
   if(sum(mm==WELL))
      mm
      x= C(:, mm);
      x= x-mean(x);
%      figure(1);plot(abs(fft(x).^2))
      x= filter(Bh, Ah, x);
      x= filter(Bl, Al, x);
      x= x/sqrt(energy(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(PsiRD))))
      drawnow
      X= [X x];
%      pause
   end;
end;
%X= 10*X./max(max(X));
MM= size(X, 2)
EX= energy(X)/MM/N

%%Verrauschen der Daten
SNR= 3;%dB
n= sqrt(EX/10^(SNR/10))*randn(size(X));
energy(X)/energy(n)
Y= X+n;

save('motorlow.mat', 'X', 'Y', 'PsiRD', 'PsiRS');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load motor
%load motorlow;

figure(1);clf;
%subplot(3, 1, 1);cla
subplot(2, 1, 1);cla
plot(X(:, mm0), 'Linewidth', 1, 'Linestyle', '-', 'Color', 'k')
axis([1 N, -.35 .35])
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', -.3:.1:.3)
set(gca, 'YTickLabel', ['A'; 'B'; 'C'; 'D'; 'E'; 'F'; 'G'])
text(220, .27, 'xi')
%print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorx.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3, 1, 2);cla;hold on
Si= (abs(fft(X(:, mm0))).^2);
ri= ifft(Si);
ri= real(ri.*[PsiRD(N/2+1, N/2+1:N).'; PsiRD(N/2+1, 1:N/2).']);
Wi= fft(ri);
Wi= real([Wi(N/2+1:N); Wi(1:N/2)])/N;
Si= [Si(N/2+1:N); Si(1:N/2)]/N;
plot(Si, 'Linewidth', 1, 'Linestyle', '-', 'Color', 'r')
plot(Wi, 'Linewidth', 2, 'Linestyle', '-', 'Color', 'k')
axis([1 N, 0 2])
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [-128 -64 0 63 127])
text(220, .55, 'PSDi')
box
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(2);clf;hold on
%subplot(3, 1, 3);cla;hold on
subplot(2, 1, 2);cla;hold on
S= (sum(abs(fft(X)).^2.')/MM).';
r= ifft(S);
r= real(r.*[PsiRS(N/2+1, N/2+1:N).'; PsiRS(N/2+1, 1:N/2).']);
W= fft(r);
W= real([W(N/2+1:N); W(1:N/2)])/N;
S= [S(N/2+1:N); S(1:N/2)]/N;
%plot(W, 'Linewidth', 1, 'Linestyle', '-', 'Color', 'r')
plot(S(N/2+1:N), 'Linewidth', 2, 'Linestyle', '-', 'Color', 'k')
axis([1 N/2, 0 .06])
set(gca, 'XTick', [1 N/4 N/2])
set(gca, 'XTickLabel', ['K'; 'L'; 'M'])
set(gca, 'YTick', 0:.02:.06)
set(gca, 'YTickLabel', ['U'; 'V'; 'W'; 'X'])
text(110, .05, 'PSD')
box
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorx.eps


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

RXengine= real(nm_to_nk(ml_to_nm(AX.*conj(PsiRS))));
RYengine= real(nm_to_nk(ml_to_nm(AY.*conj(PsiRS))));
%RXengine= real(nm_to_nk(ml_to_nm(AX)));
%RYengine= real(nm_to_nk(ml_to_nm(AY)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);subplot(3, 1, [1 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/4 N/2 ])
set(gca, 'YTickLabel', [127 63 0])
text(220, 10, 'RS')
subplot(3, 1, 3);plot(sum(Rengine)/N/MM, 'Linewidth', 2, 'Linestyle', '-', 'Color', 'k')
axis([1 N, 0 2])
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [-128 -64 0 63 127])
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRx.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);clf;
tf_show(rot90(RXengine(:, 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, 'RSX')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRx.eps

figure(3);clf;
tf_show(rot90(RYengine(:, 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, 'RSY')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRy.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);clf;%subplot(3, 1, [1 2])
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/4 N/2 ])
set(gca, 'YTickLabel', [127 63 0])
set(gca, 'YTickLabel', ['R'; 'S'; 'T'])
text(220, 10, 'RDi')
%subplot(3, 1, 3);plot(sum(RRD)/N, 'Linewidth', 2, 'Linestyle', '-', 'Color', 'k')
%axis([1 N, 0 2])
%set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
%set(gca, 'XTickLabel', [-128 -64 0 63 127])
print -deps2 -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)
[min(KUn) max(KUn)]










figure(3);clf;
tf_show(rot90(RD(:, 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, 'RD')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorRD.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFMA model across ensemble by MDL/AIC (A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 18;
Lmax= 3;
NB= ((1:Mmax)'+1)*(2*(0:Lmax)+1);
NA= zeros(Mmax, Lmax+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;
EE= zeros(Mmax, Lmax+1, MM);
OAIC= [];
OMDL= [];
PsiMA= tf_multiwin(N, Mmax, Lmax, 0, 2, 1);
mm= mm0;
for mm= 1:MM
   mm
   sig= X(:, mm);
   E= tfma_fitA(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-A-MDL')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorMA-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])
set(gca, 'YTickLabel', ['R'; 'S'; 'T'])
text(220, 10, 'MA-A-AIC')
print -deps2 -tiff -r300 /users/mjachan/tex/prj/TF/figures/motorMA-A-AIC.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify TFMA model across ensemble by MDL/AIC (B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mmax= 18;
Lmax= 3;
NB= ((1:Mmax)'+1)*(2*(0:Lmax)+1);
NA= zeros(Mmax, Lmax+1);
rho= 1-log(12);
lambda= .9800;
MM= 100;