%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Generate matrix X containing 100 selected realizations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;tfpmload('~/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);claplot(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);claS= 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, :))];%   pauseend;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);drawnowendARAIChist[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);drawnowendAIChist[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;endfigure(1);mesh(abs(PP))figure(2);mesh(abs(SS))