X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Ftsa-4.2.4%2Fsinvest1.m;fp=octave_packages%2Ftsa-4.2.4%2Fsinvest1.m;h=9b5a0499193d53575c6615bf99d5e0f542b81a21;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/tsa-4.2.4/sinvest1.m b/octave_packages/tsa-4.2.4/sinvest1.m new file mode 100644 index 0000000..9b5a049 --- /dev/null +++ b/octave_packages/tsa-4.2.4/sinvest1.m @@ -0,0 +1,289 @@ +%SINVEST1 shows the parameters of a time series calculated by INVEST1 +% only called by INVEST1 + +% $Id: sinvest1.m 5090 2008-06-05 08:12:04Z schloegl $ +% Copyright (C) 1998-2002,2008 by Alois Schloegl +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + +Fs=flag_implicit_samplerate; +M=size(AutoCorr,2); +oo=M-1; +while 1, + K=menu('Select','Autocovariance ACOVF','Autocorrelation ACF', ... + 'Partial ACF PACF','Coeff. of Determination R²', ... + 'Error curve',... + 'Autoregressive Parameters',... + 'Information Criteria', ... + 'Matched Filter', ... + 'Log PSD and Phase', ... + 'Poles', ... + 'Inverse Filtering', ... + 'Spectra H(f,p)', ... + 'Entropy H=ln(det(R))', ... + 'Histogram of MOPS', ... + 'end'); + subplot(111); + if K==1 + plot(0:M,AutoCov); + title('Autocovariance function ACOVF(k)'); + xlabel('Lag k'); + elseif K==2 + if exist('OCTAVE_VERSION')==5 %%%%% Fuer OCTAVE + plot(1:M,AutoCorr); + elseif strcmp(version,'MIDEVA') %%%%% fuer MatCom + plot(1:M,AutoCorr); + else %%%%% fuer Matlab + % size(AutoCorr),size(ACFsd)] + % errorbar(ones(nr,1)*(1:M),AutoCorr,ACFsd,ACFsd); + % errorbar(1:M,AutoCorr,AutoCorr-ACFsd,AutoCorr+ACFsd); + plot(1:M,AutoCorr,'b',[1,M],[-1,1;-1,1]/sqrt(min(NC)),'b:'); + if exist('OCTAVE_VERSION')<5 + legend({'ACF','1/sqrt(N)'}) + end; + end; + title('Autocorrelation function ACF(k)'); + xlabel('Lag k'); + + elseif K==3 + %rc=ARPMX(:,(1:M).*(2:M+1)/2); + %plot(1:M,PartACF); + if size(ARPMX,2)==2*Pmax, + RC=ARPMX(:,Pmax+1:2*Pmax); + else + RC=ARPMX(:,(1:M).*(2:M+1)/2); + end; + % according to http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc4463.htm + % is the 95% confidence interval 2/sqrt(N) + %plot(1:M,RC,'b',[1,M],[3,3;2,2;1,1;-1,-1;-2,-2;-3,-3]*(1/sqrt(min(NC))),'b:'); + %legend({'Part. ACF','1/sqrt(N)','2/sqrt(N) = 95% confidence interval','3/sqrt(N)'}) + plot(1:M,RC,'b',[1,M],[2,2;-2,-2]'*(1/sqrt(min(NC))),'b:'); + if exist('OCTAVE_VERSION')<5 + legend({'Part. ACF','2/sqrt(N) = 95% confidence interval'}) + end; + title('Partial Autocorrelation function PACF(k)'); + xlabel('Lag k'); + elseif K==4 + plot(0:M,(E(1)-E)/E(1)); + title('Determination of Regression R²=1-var{E}/var{Y}'); + xlabel('Model order p'); + elseif K==5 + plot(0:M,E,'r'); + if exist('OCTAVE_VERSION')<5 + v=axis; v(3)=min([v(3); 0]); axis(v); + end + title('Mean Square (prediction) Error decrease with increasing model order'); + xlabel('Model order p'); + elseif K==6 + plot(1:oo,ARPMX(:,oo/2*(oo-1)+(1:oo))','o-'); + title( ['AutoRegressive Parameters with model order' int2str(oo)]) + elseif K==7 + while 1 + K=menu('Select Criterion',... + 'Final Predection Error FPE', ... + 'Akaike Information Criterion AIC', ... + 'Bayesian Akaike Information Criterion BIC', ... + 'Schwartz''s Bayesian Crit. SBC (=MDL)', ... + 'Parzen''s CAT Criterion ', ... + 'PHI Criterion', ... + 'end'); + + if K==1 + tmp=min(2*optFPE+2,M); + oo=optFPE; + plot(0:tmp-1,FPE(1:tmp),'r',optFPE,FPE(optFPE+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optFPE,FPE(optFPE+1),sprintf('%i',optFPE)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Final Prediction Error FPE criterion'); + elseif K==2 + tmp=min(2*optAIC+2,M); + oo=optAIC; + plot(0:tmp-1,AIC(1:tmp),'r',optAIC,AIC(optAIC+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optAIC,AIC(optAIC+1),sprintf('%i',optAIC)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Akaike''s Information Criterion AIC'); + elseif K==3 + tmp=min(2*optBIC+2,M); + oo=optBIC; + plot(0:tmp-1,BIC(1:tmp),'r',optBIC,BIC(optBIC+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optBIC,BIC(optBIC+1),sprintf('%i',optBIC)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Bayesian Akaike Information Criterion BIC'); + elseif K==4 + tmp=min(2*optSBC+2,M); + oo=optSBC; + plot(0:tmp-1,SBC(1:tmp),'r',optSBC,SBC(optSBC+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optSBC,SBC(optSBC+1),sprintf('%i',optSBC)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Schwartz''s Bayesian Criterion SBC'); + %elseif K==5 + % tmp=min(2*optMDL,M); + % oo=optMDL; + % plot(0:tmp-1,MDL(1:tmp),'r',optMDL,MDL(optMDL+1),'ro'); + % v=axis; v(3)=min([v(3); 0]); axis(v); + % text(optMDL,MDL(optMDL+1),sprintf('%i',optMDL)) + % title('Minimal Description length Criterion MDL'); + elseif K==5 + tmp=min(2*optCAT+2,M); + oo=optCAT; + plot(0:tmp-1,CATcrit(1:tmp),'r',optCAT,CATcrit(optCAT+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optCAT,CATcrit(optCAT+1),sprintf('%i',optCAT)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Parzen''s CAT Criterion '); + elseif K==6 + tmp=min(2*optPHI+2,M); + oo=optPHI; + plot(0:tmp-1,PHI(1:tmp),'r',optPHI,PHI(optPHI+1),'ro'); + if exist('OCTAVE_VERSION')<5 + text(optPHI,PHI(optPHI+1),sprintf('%i',optPHI)); + v=axis; v(3)=min([v(3); 0]); axis(v); + end; + title('Phi criterion '); + elseif K==7 + %[ARP,rc,res] =durlev(sum(AutoCov(:,1:(oo+1)),1)); + ARP=ARPMX(:,oo/2*(oo-1)+(1:oo)); + break; + end; %IF + end; %WHILE + elseif K==8 + % if model order p is given then the filter parameters A are + % A=[1; -earpyw(signal,p)] + % inverse filtering is invfiltsignal=filter(A,1,signal); + + h=zeros(nr,512)'; + w=zeros(nr,512)'; + for k=1:nr, + tmp=freqz(sqrt(E(k,oo+1)),[1 -ARPMX(k,oo/2*(oo-1)+(1:oo))],512); + h(:,k)=tmp(:); + end + + plot((0:511)/512/2*Fs,abs(h)); + %plot((0:511)/512/2,abs(freqz(1,[1 -ARP]'))); + title('Matched Filter'); + xlabel('Frequency f'); + ylabel('|H(f)|'); + elseif K==9 + h=zeros(nr,512)'; + w=zeros(nr,512)'; + for k=1:nr, + tmp=freqz(sqrt(E(k,oo+1)),[1 -ARPMX(k,oo/2*(oo-1)+(1:oo))],512); + h(:,k)=tmp(:); + end; + subplot(211); + semilogy((0:511)/512/2*Fs,abs(h')) + %semilogy((0:511)/512/2,abs(freqz(1,[1 -ARP]'))) + title('Logarithmic Spectral Density Fct.'); + subplot(212); + plot((0:511)/512/2*Fs,angle(h)'); + %plot((0:511)/512/2,angle(freqz(1,[1 -ARP]'))); + ylabel('rad'); + elseif K==10 + % clf; + %r = roots([1 -ARP]); + r = roots([1 -ARPMX(:,oo/2*(oo-1)+(1:oo))]); + t = 0:1/70:2*pi; + plot(cos(t), sin(t), 'b:',real(r), imag(r), 'rx'); + + % zplane([],[1 -ARPmx(oo+1,1:oo)]); + title('Pole Diagram'); + xlabel('real(z)'); + ylabel('imag(z)'); + + MATLAB_VERSION = version; + if MATLAB_VERSION(1)=='4' + ax = gca; + tmp = get(ax,'Aspect'); + set(ax,'Aspect',[tmp(1),1]); + elseif MATLAB_VERSION(1)=='5' + ax = gca; + tmp = get(ax,'DataAspectRatio'); + set(ax,'PlotBoxAspectRatio',tmp); + end; + elseif K==11 + plot([Y(:) filter([1 -ARPMX(:,oo/2*(oo-1)+(1:oo))],1,Y(:))-max(Y)+min(Y)]); + elseif K==12 + %[tmp,ARPmx,PE]= acf2pacf(AutoCov(2:length(AutoCov))/AutoCov(1),AutoCov(1)); + %[arp,rc,PE,ARPMX] = durlev(AutoCov); + %[tmp,ARPmx]=arp2pacf(AR); + N = Pmax-1; %2*oo-1; %2^(ceil(log(oo)/log(2))); + sdf = zeros(512,N); %length(AR)); + for k = 1:N; %[k,size(sdf),N],%length(AR); + % [sdf(:,k),F] = freqz(1,[1 -ARPmx(k,1:k)],N); + [tmp,F] = freqz(1,[1 -ARPMX(1,k/2*(k+1)+(1:k))]',512); + sdf(:,k)=tmp(:); + % sdf(:,k)=sqrt(E(k+1))*sdf(:,k); + end; + mesh(F/2/pi*Fs,1:N,log10(abs(sdf)')); + zlabel('log10 |H(f,p)|');ylabel('model order p'); xlabel('frequency f [2*pi rad/s]'); + if exist('OCTAVE_VERSION')<5 + view(30,45); + end; + elseif K==13 + for k=1:M, + xxx=eig(toeplitz(AutoCov(1:k)/E(k))); + H(k) = .5*sum(log(xxx)); + H1(k)= H(k)/(k); + %xxx1=eig(toeplitz(AutoCov(1:k)/E(k)*E(1))); + %xxx1=eig(toeplitz(AutoCorr(1:k))); + %H2(k) =.5*sum(log(xxx1)); + %H3(k)=.5*sum(log(xxx1))/(k); + %if any(xxx<=0) fprintf(1,'COV positive definite for p up to %i\n',k-1); break; end; + end; + if 0 + subplot(211); + plot(0:k-1,H2'); + title('Entropy H(k) depending on number of coefficients') + subplot(212); + plot(0:k-1,H3'); + title('Entropy rate H(k) depending on number of coefficients') + else + subplot(311); + plot(0:k-1,H'); + title('Entropy H(k) depending on number of coefficients') + subplot(312); + xxx=(-diff(log(E(:))))/2; %[size(xxx)size(H) M k] + plot(0:k-1,H1','b',1:k,xxx(1:k),'g'); + %plot(0:k-1,H1','b',0:k-1,log(E(1:k))-log(E(1)),'b',1:k,xxx(1:k),'g'); + %plot(0:k-1,H1','b',1:k,xxx,'g',1:k-1,cumprod(1+xxx(1:M-1))./H(2:M)','r'); + title('Entropy rate H(k) and Entropy difference H(k)-H(k+1)') + subplot(313); + plot(1:k,log(xxx(1:k)),'g'); + title('LOG Entropy difference for p->p+1: H(k)-H(k+1)') + end; + elseif K==14, + tmp = histo3(MOPS(1:max(1,size(MOPS,1)-1),:)); + if exist('OCTAVE_VERSION')<5, + bar(tmp.X,tmp.H,'stacked'); + else + bar(tmp.X,sum(tmp.H,2)); + end; + xlabel('model order p') + if exist('OCTAVE_VERSION')<5 + %legend({'FPE','AIC','BIC','SBC','MDL','CAT','PHI','JEW','HAR'}); + legend('FPE','AIC','BIC','SBC','MDL','CAT','PHI','JEW','HAR'); + end; + elseif K==15 + break; + end; +end;