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] = ar_psd(a,v,freq,Fs,range,method,plot_type)
19 %% Calculate the power spectrum of the autoregressive model
22 %% x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
24 %% where x(n) is the output of the model and e(n) is white noise.
25 %% This function is intended for use with
26 %% [a,v,k] = arburg(x,poles,criterion)
27 %% which use the Burg (1968) method to calculate a "maximum entropy"
28 %% autoregressive model of "x". This function runs on octave and matlab.
30 %% If the "freq" argument is a vector (of frequencies) the spectrum is
31 %% calculated using the polynomial method and the "method" argument is
32 %% ignored. For scalar "freq", an integer power of 2, or "method='FFT'",
33 %% causes the spectrum to be calculated by FFT. Otherwise, the spectrum
34 %% is calculated as a polynomial. It may be computationally more
35 %% efficient to use the FFT method if length of the model is not much
36 %% smaller than the number of frequency values. The spectrum is scaled so
37 %% that spectral energy (area under spectrum) is the same as the
38 %% time-domain energy (mean square of the signal).
41 %% All but the first two arguments are optional and may be empty.
43 %% a %% [vector] list of M=(order+1) autoregressive model
44 %% %% coefficients. The first element of "ar_coeffs" is the
45 %% %% zero-lag coefficient, which always has a value of 1.
47 %% v %% [real scalar] square of the moving-average coefficient of
50 %% freq %% [real vector] frequencies at which power spectral density
52 %% %% [integer scalar] number of uniformly distributed frequency
53 %% %% values at which spectral density is calculated.
56 %% Fs %% [real scalar] sampling frequency (Hertz) [default=1]
58 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
59 %% Control-string arguments can be in any order after the other arguments.
61 %% range %% 'half', 'onesided' : frequency range of the spectrum is
62 %% %% from zero up to but not including sample_f/2. Power
63 %% %% from negative frequencies is added to the positive
64 %% %% side of the spectrum.
65 %% %% 'whole', 'twosided' : frequency range of the spectrum is
66 %% %% -sample_f/2 to sample_f/2, with negative frequencies
67 %% %% stored in "wrap around" order after the positive
68 %% %% frequencies; e.g. frequencies for a 10-point 'twosided'
69 %% %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
70 %% %% 'shift', 'centerdc' : same as 'whole' but with the first half
71 %% %% of the spectrum swapped with second half to put the
72 %% %% zero-frequency value in the middle. (See "help
73 %% %% fftshift". If "freq" is vector, 'shift' is ignored.
74 %% %% If model coefficients "ar_coeffs" are real, the default
75 %% %% range is 'half', otherwise default range is 'whole'.
77 %% method %% 'fft': use FFT to calculate power spectrum.
78 %% %% 'poly': calculate power spectrum as a polynomial of 1/z
79 %% %% N.B. this argument is ignored if the "freq" argument is a
80 %% %% vector. The default is 'poly' unless the "freq"
81 %% %% argument is an integer power of 2.
83 %% plot_type%% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
84 %% %% specifies the type of plot. The default is 'plot', which
85 %% %% means linear-linear axes. 'squared' is the same as 'plot'.
86 %% %% 'dB' plots "10*log10(psd)". This argument is ignored and a
87 %% %% spectrum is not plotted if the caller requires a returned
91 %% If returned values are not required by the caller, the spectrum
92 %% is plotted and nothing is returned.
93 %% psd %% [real vector] estimate of power-spectral density
94 %% f_out %% [real vector] frequency values
96 %% N.B. arburg runs in octave and matlab, and does not depend on octave-forge
97 %% or signal-processing-toolbox functions.
100 %% [1] Equation 2.28 from Steven M. Kay and Stanley Lawrence Marple Jr.:
101 %% "Spectrum analysis -- a modern perspective",
102 %% Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
105 function [varargout]=ar_psd(a,v,varargin)
107 %% Check fixed arguments
109 error( 'ar_psd: needs at least 2 args. Use help ar_psd.' );
110 elseif ( ~isvector(a) || length(a)<2 )
111 error( 'ar_psd: arg 1 (a) must be vector, length>=2.' );
112 elseif ( ~isscalar(v) )
113 error( 'ar_psd: arg 2 (v) must be real scalar >0.' );
115 real_model = isreal(a);
117 %% default values for optional areguments
119 user_freqs = 0; %% boolean: true for user-specified frequencies
121 %% FFT padding factor (is also frequency range divisor): 1=whole, 2=half.
122 pad_fact = 1 + real_model;
128 %% decode and check optional arguments
129 %% end_numeric_args is boolean; becomes true at 1st string arg
130 end_numeric_args = 0;
131 for iarg = 1:length(varargin)
132 arg = varargin{iarg};
133 end_numeric_args = end_numeric_args || ischar(arg);
134 %% skip empty arguments
137 %% numeric optional arguments must be first, cannot follow string args
138 %% N.B. older versions of matlab may not have the function "error" so
139 %% the user writes "function error(msg); disp(msg); end" and we need
141 elseif ( ~ischar(arg) )
142 if ( end_numeric_args )
143 error( 'ar_psd: control arg must be string.' );
145 %% first optional numeric arg is "freq"
147 user_freqs = isvector(arg) && length(arg)>1;
148 if ( ~isscalar(arg) && ~user_freqs )
149 error( 'ar_psd: arg 3 (freq) must be vector or scalar.' );
150 elseif ( ~user_freqs && ( ~isreal(arg) || ...
151 fix(arg)~=arg || arg <= 2 || arg >= 1048576 ) )
152 error('ar_psd: arg 3 (freq) must be integer >=2, <=1048576' );
153 elseif ( user_freqs && ~isreal(arg) )
154 error( 'ar_psd: arg 3 (freq) vector must be real.' );
156 freq = arg(:); % -> column vector
158 %% second optional numeric arg is "Fs" - sampling frequency
160 if ( ~isscalar(arg) || ~isreal(arg) || arg<=0 )
161 error( 'ar_psd: arg 4 (Fs) must be real positive scalar.' );
166 error( 'ar_psd: control arg must be string.' );
169 %% decode control-string arguments
170 elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
172 elseif ( strcmp(arg,'semilogx') )
174 elseif ( strcmp(arg,'semilogy') )
176 elseif ( strcmp(arg,'loglog') )
178 elseif ( strcmp(arg,'dB') )
180 elseif ( strcmp(arg,'fft') )
183 elseif ( strcmp(arg,'poly') )
186 elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
187 pad_fact = 2; % FFT zero-padding factor (pad FFT to double length)
189 elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
190 pad_fact = 1; % FFT zero-padding factor (do not pad)
192 elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
196 error( 'ar_psd: string arg: illegal value: %s', arg );
199 %% end of decoding and checking args
202 %% user provides (column) vector of frequencies
203 if ( any(abs(freq)>Fs/2) )
204 error( 'ar_psd: arg 3 (freq) cannot exceed half sampling frequency.' );
205 elseif ( pad_fact==2 && any(freq<0) )
206 error( 'ar_psd: arg 3 (freq) must be positive in onesided spectrum' );
208 freq_len = length(freq);
213 %% internally generated frequencies
215 freq = (Fs/pad_fact/freq_len) * [0:freq_len-1]';
216 %% decide which method to use (poly or FFT)
217 is_power_of_2 = rem(log(freq_len),log(2))<10.*eps;
218 use_FFT = ( ~ force_poly && is_power_of_2 ) || force_FFT;
219 fft_len = freq_len * pad_fact;
223 %% calculate denominator of Equation 2.28, Kay and Marple, ref [1]Jr.:
224 len_coeffs = length(a);
227 fft_out = fft( [ a(:); zeros(fft_len-len_coeffs,1) ] );
230 %% complex data on "half" frequency range needs -ve frequency values
231 if ( pad_fact==2 && ~real_model )
232 freq = [freq; -freq(freq_len:-1:1)];
233 fft_len = 2*freq_len;
235 fft_out = polyval( a(len_coeffs:-1:1), exp( (-i*2*pi/Fs) * freq ) );
238 %% The power spectrum (PSD) is the scaled squared reciprocal of amplitude
239 %% of the FFT/polynomial. This is NOT the reciprocal of the periodogram.
240 %% The PSD is a continuous function of frequency. For uniformly
241 %% distributed frequency values, the FFT algorithm might be the most
242 %% efficient way of calculating it.
244 psd = ( v / Fs ) ./ ( fft_out .* conj(fft_out) );
246 %% range='half' or 'onesided',
247 %% add PSD at -ve frequencies to PSD at +ve frequencies
248 %% N.B. unlike periodogram, PSD at zero frequency _is_ doubled.
250 freq = freq(1:freq_len);
252 %% real data, double the psd
253 psd = 2 * psd(1:freq_len);
255 %% complex data, FFT method, internally-generated frequencies
256 psd = psd(1:freq_len)+[psd(1); psd(fft_len:-1:freq_len+2)];
258 %% complex data, polynomial method
259 %% user-defined and internally-generated frequencies
260 psd = psd(1:freq_len)+psd(fft_len:-1:freq_len+1);
264 %% disabled for user-supplied frequencies
265 %% Shift zero-frequency to the middle (pad_fact==1)
267 len2 = fix((fft_len+1)/2);
268 psd = [psd(len2+1:fft_len); psd(1:len2)];
269 freq = [freq(len2+1:fft_len)-Fs; freq(1:len2)];
272 %% Plot the spectrum if there are no return variables.
276 elseif ( nargout == 1 )
279 if ( plot_type == 1 )
281 elseif ( plot_type == 2 )
283 elseif ( plot_type == 3 )
285 elseif ( plot_type == 4 )
287 elseif ( plot_type == 5 )
288 plot(freq,10*log10(psd));