Mercurial > forge

comparison main/signal/levinson.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children e26254f9e5dc
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 ## Copyright (C) 1999 Paul Kienzle
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 ##
17 ## Based on:
18 ## yulewalker.m
19 ## Copyright (C) 1995 (GPL)
20 ## Friedrich Leisch <Friedrich.Leisch@ci.tuwien.ac.at>
21
22 ## usage: [a, v, ref] = levinson (acf [, p])
23 ##
24 ## Use the Durbin-Levinson algorithm to solve:
25 ## toeplitz(acf(1:p)) * x = -acf(2:p+1).
26 ## The solution [1, x'] is the denominator of an all pole filter
27 ## approximation to the signal x which generated the autocorrelation
28 ## function acf.
29 ##
30 ## acf is the autocorrelation function for lags 0 to p.
31 ## p defaults to length(acf)-1.
32 ## Returns
33 ## a=[1, x'] the denominator filter coefficients.
34 ## v= variance of the white noise = square of the numerator constant
35 ## ref = reflection coefficients = coefficients of the lattice
36 ## implementation of the filter
37 ## Use freqz(sqrt(v),a) to plot the power spectrum.
38
39 ## Author: PAK <pkienzle@kienzle.powernet.co.uk>
40
41 ## TODO: Matlab doesn't return reflection coefficients and
42 ## TODO: errors in addition to the polynomial a.
43 ## TODO: What is the difference between aryule, levinson,
44 ## TODO: ac2poly, ac2ar, lpc, etc.?
45
46 function [a, v, ref] = levinson (acf, p)
47
48 if( columns (acf) > 1 ) acf=acf'; endif
49 if (nargin == 1) p = length(acf) - 1; endif
50
51 if nargout < 3 && p < 100
52 ## direct solution [O(p^3), but no loops so slightly faster for small p]
53 R = toeplitz(acf(1:p), conj(acf(1:p)));
54 a = R \ -acf(2:p+1);
55 a = [ 1, a' ];
56 v = sum(a'.*acf(1:p+1));
57 else
58 ## durbin-levinson [O(p^2), so significantly faster for large p]
59 ref = zeros (1, p);
60 g = acf(2) / acf(1);
61 a = [ g ];
62 v = ( 1 - g^2 ) * acf(1);
63 ref(1) = g;
64 for t = 2 : p
65 g = (acf(t+1) - a * acf(2:t)) / v;
66 a = [ g, a-g*a(t-1:-1:1) ];
67 v = v * ( 1 - g^2 ) ;
68 ref(t) = g;
69 endfor
70 a = [1, -a(p:-1:1)];
71 endif
72
73 endfunction