Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

[90d51d]: inst / arburg.m Maximize Restore History

Download this file

arburg.m    245 lines (239 with data), 9.7 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
## Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
##
## 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 <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{a}, @var{v}, @var{k}] =} arburg (@var{x}, @var{poles})
## @deftypefnx {Function File} {[@var{a}, @var{v}, @var{k}] =} arburg (@var{x}, @var{poles}, @var{criterion})
##
## Calculate coefficients of an autoregressive (AR) model of complex data
## @var{x} using the whitening lattice-filter method of Burg (1968). The
## inverse of the model is a moving-average filter which reduces @var{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 @code{ar_psd}
## calculates the power spectrum of the AR model.
##
## ARGUMENTS:
## @itemize
## @item
## @var{x}
## sampled data
## @item
## @var{poles}
## number of poles in the AR model or limit to the number of poles if a
## valid @var{criterion} is provided.
## @item
## @var{criterion}
## 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
## @end itemize
##
## RETURNED VALUES:
## @itemize
## @item
## @var{a}
## list of (P+1) autoregression coefficients; for data input @math{x(n)} and
## white noise @math{e(n)}, the model is
##
## @example
## @group
## P+1
## x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
## k=1
## @end group
## @end example
##
## @var{v}
## mean square of residual noise from the whitening operation of the Burg
## lattice filter.
## @item
## @var{k}
## reflection coefficients defining the lattice-filter embodiment of the model
## @end itemize
##
## 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) 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
##
## @seealso{ar_psd}
## @end deftypefn
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 "poles<length(x)"
## AKICc and AICc require "length(x)-poles-2">0
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' );
endif
else
use_inf_crit = 0;
use_FPE = 0;
endif
##
## 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;
endif
## 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;
endif
##
## (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;
endif
endif
## 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;
endif
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;
endif
endfor
varargout{1} = [1 a];
varargout{2} = v;
if ( nargout>=3 )
varargout{3} = k;
endif
endif
endfunction