X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Farburg.m;fp=octave_packages%2Fsignal-1.1.3%2Farburg.m;h=a13ff20821de32e30ef7716bb647d1cf7cc662db;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/signal-1.1.3/arburg.m b/octave_packages/signal-1.1.3/arburg.m new file mode 100644 index 0000000..a13ff20 --- /dev/null +++ b/octave_packages/signal-1.1.3/arburg.m @@ -0,0 +1,228 @@ +%% Copyright (C) 2006 Peter V. Lanspeary +%% +%% 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 . + +%% [a,v,k] = arburg(x,poles,criterion) +%% +%% Calculate coefficients of an autoregressive (AR) model of complex data +%% "x" using the whitening lattice-filter method of Burg (1968). The inverse +%% of the model is a moving-average filter which reduces "x" to white noise. +%% The power spectrum of the AR model is an estimate of the maximum +%% entropy power spectrum of the data. The function "ar_psd" calculates the +%% power spectrum of the AR model. +%% +%% ARGUMENTS: +%% x %% [vector] sampled data +%% +%% poles %% [integer scalar] number of poles in the AR model or +%% %% limit to the number of poles if a +%% %% valid "stop_crit" is provided. +%% +%% criterion %% [optional string arg] model-selection criterion. Limits +%% %% the number of poles so that spurious poles are not +%% %% added when the whitened data has no more information +%% %% in it (see Kay & Marple, 1981). Recognised values are +%% %% 'AKICc' -- approximate corrected Kullback information +%% %% criterion (recommended), +%% %% 'KIC' -- Kullback information criterion +%% %% 'AICc' -- corrected Akaike information criterion +%% %% 'AIC' -- Akaike information criterion +%% %% 'FPE' -- final prediction error" criterion +%% %% The default is to NOT use a model-selection criterion +%% +%% RETURNED VALUES: +%% a %% [polynomial/vector] list of (P+1) autoregression coeffic- +%% %% ients; for data input x(n) and white noise e(n), +%% %% the model is +%% %% P+1 +%% %% x(n) = sqrt(v).e(n) + SUM a(k).x(n-k) +%% %% k=1 +%% +%% v %% [real scalar] mean square of residual noise from the +%% %% whitening operation of the Burg lattice filter. +%% +%% k %% [column vector] reflection coefficients defining the +%% %% lattice-filter embodiment of the model +%% +%% HINTS: +%% (1) arburg does not remove the mean from the data. You should remove +%% the mean from the data if you want a power spectrum. A non-zero mean +%% can produce large errors in a power-spectrum estimate. See +%% "help detrend". +%% (2) If you don't know what the value of "poles" should be, choose the +%% largest (reasonable) value you could want and use the recommended +%% value, criterion='AKICc', so that arburg can find it. +%% E.g. arburg(x,64,'AKICc') +%% The AKICc has the least bias and best resolution of the available +%% model-selection criteria. +%% (3) arburg runs in octave and matlab, does not depend on octave forge +%% or signal-processing-toolbox functions. +%% (4) Autoregressive and moving-average filters are stored as polynomials +%% which, in matlab, are row vectors. +%% +%% NOTE ON SELECTION CRITERION +%% AIC, AICc, KIC and AKICc are based on information theory. They attempt +%% to balance the complexity (or length) of the model against how well the +%% model fits the data. AIC and KIC are biassed estimates of the asymmetric +%% and the symmetric Kullback-Leibler divergence respectively. AICc and +%% AKICc attempt to correct the bias. See reference [4]. +%% +%% +%% REFERENCES +%% [1] John Parker Burg (1968) +%% "A new analysis technique for time series data", +%% NATO advanced study Institute on Signal Processing with Emphasis on +%% Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968. +%% +%% [2] Steven M. Kay and Stanley Lawrence Marple Jr.: +%% "Spectrum analysis -- a modern perspective", +%% Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981 +%% +%% [3] William H. Press and Saul A. Teukolsky and William T. Vetterling and +%% Brian P. Flannery +%% "Numerical recipes in C, The art of scientific computing", 2nd edition, +%% Cambridge University Press, 2002 --- Section 13.7. +%% +%% [4] Abd-Krim Seghouane and Maiza Bekara +%% "A small sample model selection criterion based on Kullback's symmetric +%% divergence", IEEE Transactions on Signal Processing, +%% Vol. 52(12), pp 3314-3323, Dec. 2004 + + +function [varargout] = arburg( x, poles, criterion ) +%% +%% Check arguments +if ( nargin < 2 ) + error( 'arburg(x,poles): Need at least 2 args.' ); +elseif ( ~isvector(x) || length(x) < 3 ) + error( 'arburg: arg 1 (x) must be vector of length >2.' ); +elseif ( ~isscalar(poles) || ~isreal(poles) || fix(poles)~=poles || poles<=0.5) + error( 'arburg: arg 2 (poles) must be positive integer.' ); +elseif ( poles >= length(x)-2 ) + %% lattice-filter algorithm requires "poles0 + error( 'arburg: arg 2 (poles) must be less than length(x)-2.' ); +elseif ( nargin>2 && ~isempty(criterion) && ... + (~ischar(criterion) || size(criterion,1)~=1 ) ) + error( 'arburg: arg 3 (criterion) must be string.' ); +else + %% + %% Set the model-selection-criterion flags. + %% is_AKICc, isa_KIC and is_corrected are short-circuit flags + if ( nargin > 2 && ~isempty(criterion) ) + is_AKICc = strcmp(criterion,'AKICc'); %% AKICc + isa_KIC = is_AKICc || strcmp(criterion,'KIC'); %% KIC or AKICc + is_corrected = is_AKICc || strcmp(criterion,'AICc'); %% AKICc or AICc + use_inf_crit = is_corrected || isa_KIC || strcmp(criterion,'AIC'); + use_FPE = strcmp(criterion,'FPE'); + if ( ~use_inf_crit && ~use_FPE ) + error( 'arburg: value of arg 3 (criterion) not recognised' ); + end + else + use_inf_crit = 0; + use_FPE = 0; + end + %% + %% f(n) = forward prediction error + %% b(n) = backward prediction error + %% Storage of f(n) and b(n) is a little tricky. Because f(n) is always + %% combined with b(n-1), f(1) and b(N) are never used, and therefore are + %% not stored. Not storing unused data makes the calculation of the + %% reflection coefficient look much cleaner :) + %% N.B. {initial v} = {error for zero-order model} = + %% {zero-lag autocorrelation} = E(x*conj(x)) = x*x'/N + %% E = expectation operator + N = length(x); + k = []; + if ( size(x,1) > 1 ) % if x is column vector + f = x(2:N); + b = x(1:N-1); + v = real(x'*x) / N; + else % if x is row vector + f = x(2:N).'; + b = x(1:N-1).'; + v = real(x*x') / N; + end + %% new_crit/old_crit is the mode-selection criterion + new_crit = abs(v); + old_crit = 2 * new_crit; + for p = 1:poles + %% + %% new reflection coeff = -2* E(f.conj(b)) / ( E(f^2)+E(b(^2) ) + last_k= -2 * (b' * f) / ( f' * f + b' * b); + %% Levinson-Durbin recursion for residual + new_v = v * ( 1.0 - real(last_k * conj(last_k)) ); + if ( p > 1 ) + %% + %% Apply the model-selection criterion and break out of loop if it + %% increases (rather than decreases). + %% Do it before we update the old model "a" and "v". + %% + %% * Information Criterion (AKICc, KIC, AICc, AIC) + if ( use_inf_crit ) + old_crit = new_crit; + %% AKICc = log(new_v)+p/N/(N-p)+(3-(p+2)/N)*(p+1)/(N-p-2); + %% KIC = log(new_v)+ 3 *(p+1)/N; + %% AICc = log(new_v)+ 2 *(p+1)/(N-p-2); + %% AIC = log(new_v)+ 2 *(p+1)/N; + %% -- Calculate KIC, AICc & AIC by using is_AKICc, is_KIC and + %% is_corrected to "short circuit" the AKICc calculation. + %% The extra 4--12 scalar arithmetic ops should be quicker than + %% doing if...elseif...elseif...elseif...elseif. + new_crit = log(new_v) + is_AKICc*p/N/(N-p) + ... + (2+isa_KIC-is_AKICc*(p+2)/N) * (p+1) / (N-is_corrected*(p+2)); + if ( new_crit > old_crit ) + break; + end + %% + %% (FPE) Final prediction error + elseif ( use_FPE ) + old_crit = new_crit; + new_crit = new_v * (N+p+1)/(N-p-1); + if ( new_crit > old_crit ) + break; + end + end + %% Update model "a" and "v". + %% Use Levinson-Durbin recursion formula (for complex data). + a = [ prev_a + last_k .* conj(prev_a(p-1:-1:1)) last_k ]; + else %% if( p==1 ) + a = last_k; + end + k = [ k; last_k ]; + v = new_v; + if ( p < poles ) + prev_a = a; + %% calculate new prediction errors (by recursion): + %% f(p,n) = f(p-1,n) + k * b(p-1,n-1) n=2,3,...n + %% b(p,n) = b(p-1,n-1) + conj(k) * f(p-1,n) n=2,3,...n + %% remember f(p,1) is not stored, so don't calculate it; make f(p,2) + %% the first element in f. b(p,n) isn't calculated either. + nn = N-p; + new_f = f(2:nn) + last_k * b(2:nn); + b = b(1:nn-1) + conj(last_k) * f(1:nn-1); + f = new_f; + end + end + %% end of for loop + %% + varargout{1} = [1 a]; + varargout{2} = v; + if ( nargout>=3 ) + varargout{3} = k; + end +end +end + +%%