octave_packages/tsa-4.2.4/sinvest1.m

   1 %SINVEST1 shows the parameters of a time series calculated by INVEST1
   2 % only called by INVEST1
   3
   4 %       $Id: sinvest1.m 5090 2008-06-05 08:12:04Z schloegl $
   5 %       Copyright (C) 1998-2002,2008 by Alois Schloegl <a.schloegl@ieee.org>
   6 %
   7 %    This program is free software: you can redistribute it and/or modify
   8 %    it under the terms of the GNU General Public License as published by
   9 %    the Free Software Foundation, either version 3 of the License, or
  10 %    (at your option) any later version.
  11 %
  12 %    This program is distributed in the hope that it will be useful,
  13 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 %    GNU General Public License for more details.
  16 %
  17 %    You should have received a copy of the GNU General Public License
  18 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  19
  20 Fs=flag_implicit_samplerate;
  21 M=size(AutoCorr,2);
  22 oo=M-1;
  23 while 1,
  24         K=menu('Select','Autocovariance ACOVF','Autocorrelation ACF', ...
  25                 'Partial ACF PACF','Coeff. of Determination R²', ...
  26                 'Error curve',...
  27                 'Autoregressive Parameters',...
  28                 'Information Criteria', ...
  29                 'Matched Filter', ...
  30                 'Log PSD and Phase', ...
  31                 'Poles', ...
  32                 'Inverse Filtering', ...
  33                 'Spectra H(f,p)', ...
  34                 'Entropy H=ln(det(R))', ...
  35                 'Histogram of MOPS', ...
  36                 'end');
  37         subplot(111);
  38         if K==1
  39                 plot(0:M,AutoCov);
  40                 title('Autocovariance function ACOVF(k)');
  41                 xlabel('Lag k');
  42         elseif K==2
  43                 if exist('OCTAVE_VERSION')==5             %%%%% Fuer OCTAVE
  44                         plot(1:M,AutoCorr);
  45                 elseif strcmp(version,'MIDEVA')           %%%%% fuer MatCom
  46                         plot(1:M,AutoCorr);
  47                 else                                      %%%%% fuer Matlab
  48                         %               size(AutoCorr),size(ACFsd)]
  49                         %                errorbar(ones(nr,1)*(1:M),AutoCorr,ACFsd,ACFsd);
  50                         %               errorbar(1:M,AutoCorr,AutoCorr-ACFsd,AutoCorr+ACFsd);
  51                         plot(1:M,AutoCorr,'b',[1,M],[-1,1;-1,1]/sqrt(min(NC)),'b:');
  52                         if exist('OCTAVE_VERSION')<5
  53                                 legend({'ACF','1/sqrt(N)'})
  54                         end;
  55                 end;
  56                 title('Autocorrelation function ACF(k)');
  57                 xlabel('Lag k');
  58
  59         elseif K==3
  60                 %rc=ARPMX(:,(1:M).*(2:M+1)/2);
  61                 %plot(1:M,PartACF);
  62                 if size(ARPMX,2)==2*Pmax,
  63                         RC=ARPMX(:,Pmax+1:2*Pmax);
  64                 else
  65                         RC=ARPMX(:,(1:M).*(2:M+1)/2);
  66                 end;
  67                 % according to http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc4463.htm
  68                 % is the 95% confidence interval 2/sqrt(N)
  69                 %plot(1:M,RC,'b',[1,M],[3,3;2,2;1,1;-1,-1;-2,-2;-3,-3]*(1/sqrt(min(NC))),'b:');
  70                 %legend({'Part. ACF','1/sqrt(N)','2/sqrt(N) = 95% confidence interval','3/sqrt(N)'})
  71                 plot(1:M,RC,'b',[1,M],[2,2;-2,-2]'*(1/sqrt(min(NC))),'b:');
  72                 if exist('OCTAVE_VERSION')<5
  73                         legend({'Part. ACF','2/sqrt(N) = 95% confidence interval'})
  74                 end;
  75                 title('Partial Autocorrelation function PACF(k)');
  76                 xlabel('Lag k');
  77         elseif K==4
  78                 plot(0:M,(E(1)-E)/E(1));
  79                 title('Determination of Regression R²=1-var{E}/var{Y}');
  80                 xlabel('Model order p');
  81         elseif K==5
  82                 plot(0:M,E,'r');
  83                 if exist('OCTAVE_VERSION')<5
  84                         v=axis; v(3)=min([v(3); 0]); axis(v);
  85                 end
  86                 title('Mean Square (prediction) Error decrease with increasing model order');
  87                 xlabel('Model order p');
  88         elseif K==6
  89                 plot(1:oo,ARPMX(:,oo/2*(oo-1)+(1:oo))','o-');
  90                 title( ['AutoRegressive Parameters with model order' int2str(oo)])
  91         elseif K==7
  92                 while 1
  93                         K=menu('Select Criterion',...
  94                                 'Final Predection Error FPE', ...
  95                                 'Akaike Information Criterion AIC', ...
  96                                 'Bayesian Akaike Information Criterion BIC', ...
  97                                 'Schwartz''s Bayesian Crit. SBC (=MDL)', ...
  98                                 'Parzen''s CAT Criterion ', ...
  99                                 'PHI Criterion', ...
 100                                 'end');
 101
 102                         if K==1
 103                                 tmp=min(2*optFPE+2,M);
 104                                 oo=optFPE;
 105                                 plot(0:tmp-1,FPE(1:tmp),'r',optFPE,FPE(optFPE+1),'ro');
 106                                 if exist('OCTAVE_VERSION')<5
 107                                         text(optFPE,FPE(optFPE+1),sprintf('%i',optFPE));
 108                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 109                                 end;
 110                                 title('Final Prediction Error FPE criterion');
 111                         elseif K==2
 112                                 tmp=min(2*optAIC+2,M);
 113                                 oo=optAIC;
 114                                 plot(0:tmp-1,AIC(1:tmp),'r',optAIC,AIC(optAIC+1),'ro');
 115                                 if exist('OCTAVE_VERSION')<5
 116                                         text(optAIC,AIC(optAIC+1),sprintf('%i',optAIC));
 117                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 118                                 end;
 119                                 title('Akaike''s Information Criterion AIC');
 120                         elseif K==3
 121                                 tmp=min(2*optBIC+2,M);
 122                                 oo=optBIC;
 123                                 plot(0:tmp-1,BIC(1:tmp),'r',optBIC,BIC(optBIC+1),'ro');
 124                                 if exist('OCTAVE_VERSION')<5
 125                                         text(optBIC,BIC(optBIC+1),sprintf('%i',optBIC));
 126                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 127                                 end;
 128                                 title('Bayesian Akaike Information Criterion BIC');
 129                         elseif K==4
 130                                 tmp=min(2*optSBC+2,M);
 131                                 oo=optSBC;
 132                                 plot(0:tmp-1,SBC(1:tmp),'r',optSBC,SBC(optSBC+1),'ro');
 133                                 if exist('OCTAVE_VERSION')<5
 134                                         text(optSBC,SBC(optSBC+1),sprintf('%i',optSBC));
 135                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 136                                 end;
 137                                 title('Schwartz''s Bayesian Criterion SBC');
 138                                 %elseif K==5
 139                                 %       tmp=min(2*optMDL,M);
 140                                 %       oo=optMDL;
 141                                 %       plot(0:tmp-1,MDL(1:tmp),'r',optMDL,MDL(optMDL+1),'ro');
 142                                 %       v=axis; v(3)=min([v(3); 0]); axis(v);
 143                                 %       text(optMDL,MDL(optMDL+1),sprintf('%i',optMDL))
 144                                 %       title('Minimal Description length Criterion MDL');
 145                         elseif K==5
 146                                 tmp=min(2*optCAT+2,M);
 147                                 oo=optCAT;
 148                                 plot(0:tmp-1,CATcrit(1:tmp),'r',optCAT,CATcrit(optCAT+1),'ro');
 149                                 if exist('OCTAVE_VERSION')<5
 150                                         text(optCAT,CATcrit(optCAT+1),sprintf('%i',optCAT));
 151                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 152                                 end;
 153                                 title('Parzen''s CAT Criterion ');
 154                         elseif K==6
 155                                 tmp=min(2*optPHI+2,M);
 156                                 oo=optPHI;
 157                                 plot(0:tmp-1,PHI(1:tmp),'r',optPHI,PHI(optPHI+1),'ro');
 158                                 if exist('OCTAVE_VERSION')<5
 159                                         text(optPHI,PHI(optPHI+1),sprintf('%i',optPHI));
 160                                         v=axis; v(3)=min([v(3); 0]); axis(v);
 161                                 end;
 162                                 title('Phi criterion ');
 163                         elseif K==7
 164                                 %[ARP,rc,res] =durlev(sum(AutoCov(:,1:(oo+1)),1));
 165                                 ARP=ARPMX(:,oo/2*(oo-1)+(1:oo));
 166                                 break;
 167                         end;    %IF
 168                 end;    %WHILE
 169         elseif K==8
 170                 % if model order p is given then the filter parameters A are
 171                 % A=[1; -earpyw(signal,p)]
 172                 % inverse filtering is    invfiltsignal=filter(A,1,signal);
 173
 174                 h=zeros(nr,512)';
 175                 w=zeros(nr,512)';
 176                 for k=1:nr,
 177                         tmp=freqz(sqrt(E(k,oo+1)),[1 -ARPMX(k,oo/2*(oo-1)+(1:oo))],512);
 178                         h(:,k)=tmp(:);
 179                 end
 180
 181                 plot((0:511)/512/2*Fs,abs(h));
 182                 %plot((0:511)/512/2,abs(freqz(1,[1 -ARP]')));
 183                 title('Matched Filter');
 184                 xlabel('Frequency f');
 185                 ylabel('|H(f)|');
 186         elseif K==9
 187                 h=zeros(nr,512)';
 188                 w=zeros(nr,512)';
 189                 for k=1:nr,
 190                         tmp=freqz(sqrt(E(k,oo+1)),[1 -ARPMX(k,oo/2*(oo-1)+(1:oo))],512);
 191                         h(:,k)=tmp(:);
 192                 end;
 193                 subplot(211);
 194                 semilogy((0:511)/512/2*Fs,abs(h'))
 195                 %semilogy((0:511)/512/2,abs(freqz(1,[1 -ARP]')))
 196                 title('Logarithmic Spectral Density Fct.');
 197                 subplot(212);
 198                 plot((0:511)/512/2*Fs,angle(h)');
 199                 %plot((0:511)/512/2,angle(freqz(1,[1 -ARP]')));
 200                 ylabel('rad');
 201         elseif K==10
 202                 %       clf;
 203                 %r = roots([1 -ARP]);
 204                 r = roots([1 -ARPMX(:,oo/2*(oo-1)+(1:oo))]);
 205                 t = 0:1/70:2*pi;
 206                 plot(cos(t), sin(t), 'b:',real(r), imag(r), 'rx');
 207
 208                 %       zplane([],[1 -ARPmx(oo+1,1:oo)]);
 209                 title('Pole Diagram');
 210                 xlabel('real(z)');
 211                 ylabel('imag(z)');
 212
 213                 MATLAB_VERSION = version;
 214                 if MATLAB_VERSION(1)=='4'
 215                         ax = gca;
 216                         tmp = get(ax,'Aspect');
 217                         set(ax,'Aspect',[tmp(1),1]);
 218                 elseif MATLAB_VERSION(1)=='5'
 219                         ax = gca;
 220                         tmp = get(ax,'DataAspectRatio');
 221                         set(ax,'PlotBoxAspectRatio',tmp);
 222                 end;
 223         elseif K==11
 224                 plot([Y(:) filter([1 -ARPMX(:,oo/2*(oo-1)+(1:oo))],1,Y(:))-max(Y)+min(Y)]);
 225         elseif K==12
 226                 %[tmp,ARPmx,PE]= acf2pacf(AutoCov(2:length(AutoCov))/AutoCov(1),AutoCov(1));
 227                 %[arp,rc,PE,ARPMX] = durlev(AutoCov);
 228                 %[tmp,ARPmx]=arp2pacf(AR);
 229                 N = Pmax-1; %2*oo-1; %2^(ceil(log(oo)/log(2)));
 230                 sdf = zeros(512,N); %length(AR));
 231                 for k = 1:N; %[k,size(sdf),N],%length(AR);
 232                         % [sdf(:,k),F] = freqz(1,[1 -ARPmx(k,1:k)],N);
 233                         [tmp,F] = freqz(1,[1 -ARPMX(1,k/2*(k+1)+(1:k))]',512);
 234                         sdf(:,k)=tmp(:);
 235                         % sdf(:,k)=sqrt(E(k+1))*sdf(:,k);
 236                 end;
 237                 mesh(F/2/pi*Fs,1:N,log10(abs(sdf)'));
 238                 zlabel('log10 |H(f,p)|');ylabel('model order p'); xlabel('frequency f [2*pi rad/s]');
 239                 if exist('OCTAVE_VERSION')<5
 240                         view(30,45);
 241                 end;
 242         elseif K==13
 243                 for k=1:M,
 244                         xxx=eig(toeplitz(AutoCov(1:k)/E(k)));
 245                         H(k) = .5*sum(log(xxx));
 246                         H1(k)= H(k)/(k);
 247                         %xxx1=eig(toeplitz(AutoCov(1:k)/E(k)*E(1)));
 248                         %xxx1=eig(toeplitz(AutoCorr(1:k)));
 249                         %H2(k) =.5*sum(log(xxx1));
 250                         %H3(k)=.5*sum(log(xxx1))/(k);
 251                         %if any(xxx<=0) fprintf(1,'COV positive definite for p up to %i\n',k-1); break; end;
 252                 end;
 253                 if 0
 254                         subplot(211);
 255                         plot(0:k-1,H2');
 256                         title('Entropy H(k) depending on number of coefficients')
 257                         subplot(212);
 258                         plot(0:k-1,H3');
 259                         title('Entropy rate H(k) depending on number of coefficients')
 260                 else
 261                         subplot(311);
 262                         plot(0:k-1,H');
 263                         title('Entropy H(k) depending on number of coefficients')
 264                         subplot(312);
 265                         xxx=(-diff(log(E(:))))/2; %[size(xxx)size(H) M k]
 266                         plot(0:k-1,H1','b',1:k,xxx(1:k),'g');
 267                         %plot(0:k-1,H1','b',0:k-1,log(E(1:k))-log(E(1)),'b',1:k,xxx(1:k),'g');
 268                         %plot(0:k-1,H1','b',1:k,xxx,'g',1:k-1,cumprod(1+xxx(1:M-1))./H(2:M)','r');
 269                         title('Entropy rate H(k) and Entropy difference H(k)-H(k+1)')
 270                         subplot(313);
 271                         plot(1:k,log(xxx(1:k)),'g');
 272                         title('LOG Entropy difference for p->p+1: H(k)-H(k+1)')
 273                 end;
 274         elseif K==14,
 275                 tmp = histo3(MOPS(1:max(1,size(MOPS,1)-1),:));
 276                 if exist('OCTAVE_VERSION')<5,
 277                         bar(tmp.X,tmp.H,'stacked');
 278                 else
 279                         bar(tmp.X,sum(tmp.H,2));
 280                 end;
 281                 xlabel('model order p')
 282                 if exist('OCTAVE_VERSION')<5
 283                         %legend({'FPE','AIC','BIC','SBC','MDL','CAT','PHI','JEW','HAR'});
 284                         legend('FPE','AIC','BIC','SBC','MDL','CAT','PHI','JEW','HAR');
 285                 end;
 286         elseif K==15
 287                 break;
 288         end;
 289 end;