1 %% Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
3 %% This program is free software; you can redistribute it and/or modify it under
4 %% the terms of the GNU General Public License as published by the Free Software
5 %% Foundation; either version 3 of the License, or (at your option) any later
8 %% This program is distributed in the hope that it will be useful, but WITHOUT
9 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 %% You should have received a copy of the GNU General Public License along with
14 %% this program; if not, see <http://www.gnu.org/licenses/>.
17 %% [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)
19 %% Calculate Burg maximum-entropy power spectral density.
20 %% The functions "arburg" and "ar_psd" do all the work.
21 %% See "help arburg" and "help ar_psd" for further details.
24 %% All but the first two arguments are optional and may be empty.
25 %% x %% [vector] sampled data
27 %% poles %% [integer scalar] required number of poles of the AR model
29 %% freq %% [real vector] frequencies at which power spectral density
31 %% %% [integer scalar] number of uniformly distributed frequency
32 %% %% values at which spectral density is calculated.
35 %% Fs %% [real scalar] sampling frequency (Hertz) [default=1]
38 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
39 %% Control-string arguments can be in any order after the other arguments.
42 %% range %% 'half', 'onesided' : frequency range of the spectrum is
43 %% %% from zero up to but not including sample_f/2. Power
44 %% %% from negative frequencies is added to the positive
45 %% %% side of the spectrum.
46 %% %% 'whole', 'twosided' : frequency range of the spectrum is
47 %% %% -sample_f/2 to sample_f/2, with negative frequencies
48 %% %% stored in "wrap around" order after the positive
49 %% %% frequencies; e.g. frequencies for a 10-point 'twosided'
50 %% %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
51 %% %% 'shift', 'centerdc' : same as 'whole' but with the first half
52 %% %% of the spectrum swapped with second half to put the
53 %% %% zero-frequency value in the middle. (See "help
54 %% %% fftshift". If "freq" is vector, 'shift' is ignored.
55 %% %% If model coefficients "ar_coeffs" are real, the default
56 %% %% range is 'half', otherwise default range is 'whole'.
58 %% method %% 'fft': use FFT to calculate power spectral density.
59 %% %% 'poly': calculate spectral density as a polynomial of 1/z
60 %% %% N.B. this argument is ignored if the "freq" argument is a
61 %% %% vector. The default is 'poly' unless the "freq"
62 %% %% argument is an integer power of 2.
64 %% plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
65 %% %% specifies the type of plot. The default is 'plot', which
66 %% %% means linear-linear axes. 'squared' is the same as 'plot'.
67 %% %% 'dB' plots "10*log10(psd)". This argument is ignored and a
68 %% %% spectrum is not plotted if the caller requires a returned
71 %% criterion %% [optional string arg] model-selection criterion. Limits
72 %% %% the number of poles so that spurious poles are not
73 %% %% added when the whitened data has no more information
74 %% %% in it (see Kay & Marple, 1981). Recognised values are
75 %% %% 'AKICc' -- approximate corrected Kullback information
76 %% %% criterion (recommended),
77 %% %% 'KIC' -- Kullback information criterion
78 %% %% 'AICc' -- corrected Akaike information criterion
79 %% %% 'AIC' -- Akaike information criterion
80 %% %% 'FPE' -- final prediction error" criterion
81 %% %% The default is to NOT use a model-selection criterion
84 %% If return values are not required by the caller, the spectrum
85 %% is plotted and nothing is returned.
86 %% psd %% [real vector] power-spectral density estimate
87 %% f_out %% [real vector] frequency values
90 %% This function is a wrapper for arburg and ar_psd.
91 %% See "help arburg", "help ar_psd".
93 function [psd,f_out]=pburg(x,poles,varargin)
96 error( 'pburg: need at least 2 args. Use "help pburg"' );
98 nvarargin=length(varargin);
101 %% Search for a "criterion" arg. If found, remove it
102 %% from "varargin" list and feed it to arburg instead.
103 for iarg = 1: nvarargin
104 arrgh = varargin{iarg};
105 if ( ischar(arrgh) && ( strcmp(arrgh,'AKICc') ||...
106 strcmp(arrgh,'KIC') || strcmp(arrgh,'AICc') ||...
107 strcmp(arrgh,'AIC') || strcmp(arrgh,'FPE') ) )
117 [ar_coeffs,residual]=arburg(x,poles,criterion);
119 ar_psd(ar_coeffs,residual,varargin{:});
120 elseif ( nargout==1 )
121 psd = ar_psd(ar_coeffs,residual,varargin{:});
122 elseif ( nargout>=2 )
123 [psd,f_out] = ar_psd(ar_coeffs,residual,varargin{:});
129 %! rand('seed',2038014164);
130 %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
131 %! signal = detrend(filter(0.70181,a,rand(1,16384)));
132 %! % frequency shift by modulating with exp(j.omega.t)
133 %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
136 %! pburg(signal,3,[],Fs);
137 %! input('Onesided 3-pole spectrum. Press ENTER', 's' );
138 %! pburg(signal,4,[],Fs,'whole');
139 %! input('Twosided 4-pole spectrum of same data. Press ENTER', 's' );
140 %! pburg(signal,5,128,Fs,'shift', 'semilogy');
141 %! input('Twosided, centred zero-frequency, 5-pole. Press ENTER', 's' );
142 %! pburg(skewed,7,128,Fs,'AKICc','shift','semilogy');
143 %! input('Complex data, AKICc chooses no. of poles. Press ENTER', 's' );
144 %! user_freq=[-0.2:0.02:0.2]*Fs;
145 %! pburg(skewed,7,user_freq,Fs,'AKICc','semilogy');
146 %! input('User-specified frequency values. Press ENTER', 's' );