1## Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
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 3 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; see the file COPYING. If not, see
15## <https://www.gnu.org/licenses/>.
16
17## -*- texinfo -*-
18## @deftypefn  {Function File} {} ar_psd (@var{a}, @var{v})
19## @deftypefnx {Function File} {} ar_psd (@var{a}, @var{v}, @var{freq})
20## @deftypefnx {Function File} {} ar_psd (@var{a}, @var{v}, @var{freq}, @var{Fs})
21## @deftypefnx {Function File} {} ar_psd (@dots{}, @var{range})
22## @deftypefnx {Function File} {} ar_psd (@dots{}, @var{method})
23## @deftypefnx {Function File} {} ar_psd (@dots{}, @var{plottype})
24## @deftypefnx {Function File} {[@var{psd}, @var{f_out}] =} ar_psd (@dots{})
25##
26## Calculate the power spectrum of the autoregressive model
27##
28## @example
29## @group
30##                        M
31## x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
32##                       k=1
33## @end group
34## @end example
35##
36## where @math{x(n)} is the output of the model and @math{e(n)} is white noise.
37## This function is intended for use with
38## @code{[a, v, k] = arburg (x, poles, criterion)}
39## which use the Burg (1968) method to calculate a "maximum entropy"
40## autoregressive model of @var{x}.
41##
42## If the @var{freq} argument is a vector (of frequencies) the spectrum is
43## calculated using the polynomial method and the @var{method} argument is
44## ignored.  For scalar @var{freq}, an integer power of 2, or @var{method} =
45## "FFT", causes the spectrum to be calculated by FFT.  Otherwise, the spectrum
46## is calculated as a polynomial.  It may be computationally more
47## efficient to use the FFT method if length of the model is not much
48## smaller than the number of frequency values.  The spectrum is scaled so
49## that spectral energy (area under spectrum) is the same as the time-domain
50## energy (mean square of the signal).
51##
52## ARGUMENTS:
53## All but the first two arguments are optional and may be empty.
54## @itemize
55## @item
56## @var{a}
57## list of M=(order+1) autoregressive model
58## coefficients.  The first element of "ar_coeffs" is the
59## zero-lag coefficient, which always has a value of 1.
60## @item
61## @var{v}
62## square of the moving-average coefficient of the AR model.
63## @item
64## @var{freq}
65## frequencies at which power spectral density is calculated, or a scalar
66## indicating the number of uniformly distributed frequency
67## values at which spectral density is calculated.
68## (default = 256)
69## @item
70## @var{Fs}
71## sampling frequency (Hertz) (default=1)
72## @end itemize
73##
74## CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
75## Control-string arguments can be in any order after the other arguments.
76##
77## Range:
78##
79## 'half',  'onesided' : frequency range of the spectrum is
80## from zero up to but not including sample_f/2.  Power
81## from negative frequencies is added to the positive
82## side of the spectrum.
83## 'whole', 'twosided' : frequency range of the spectrum is
84## -sample_f/2 to sample_f/2, with negative frequencies
85## stored in "wrap around" order after the positive
86## frequencies; e.g. frequencies for a 10-point 'twosided'
87## spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
88## 'shift', 'centerdc' : same as 'whole' but with the first half
89## of the spectrum swapped with second half to put the
90## zero-frequency value in the middle. (See "help
91## fftshift". If "freq" is vector, 'shift' is ignored.
92## If model coefficients "ar_coeffs" are real, the default
93## range is 'half', otherwise default range is 'whole'.
94##
95## Method:
96##
97## 'fft':  use FFT to calculate power spectrum.
98## 'poly': calculate power spectrum as a polynomial of 1/z
99## N.B. this argument is ignored if the "freq" argument is a
100## vector.  The default is 'poly' unless the "freq"
101## argument is an integer power of 2.
102##
103## Plot type:
104##
105## 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
106## specifies the type of plot.  The default is 'plot', which
107## means linear-linear axes. 'squared' is the same as 'plot'.
108## 'dB' plots "10*log10(psd)".  This argument is ignored and a
109## spectrum is not plotted if the caller requires a returned
110## value.
111##
112## RETURNED VALUES:
113## If returned values are not required by the caller, the spectrum
114## is plotted and nothing is returned.
115## @itemize
116## @item
117## @var{psd}
118## estimate of power-spectral density
119## @item
120## @var{f_out}
121## frequency values
122## @end itemize
123##
124## REFERENCE
125## [1] Equation 2.28 from Steven M. Kay and Stanley Lawrence Marple Jr.:
126##   "Spectrum analysis -- a modern perspective",
127##   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
128##
129## @end deftypefn
130
131function varargout = ar_psd(a,v,varargin)
132
133  ##
134  ## Check fixed arguments
135  if ( nargin < 2 )
136    error( 'ar_psd: needs at least 2 args. Use help ar_psd.' );
137  elseif ( ~isvector(a) || length(a)<2 )
138    error( 'ar_psd: arg 1 (a) must be vector, length>=2.' );
139  elseif ( ~isscalar(v) )
140    error( 'ar_psd: arg 2 (v) must be real scalar >0.' );
141  else
142    real_model = isreal(a);
143  ##
144  ##  default values for optional arguments
145    freq = 256;
146    user_freqs = 0;    ## boolean: true for user-specified frequencies
147    Fs   = 1.0;
148    ##  FFT padding factor (is also frequency range divisor): 1=whole, 2=half.
149    pad_fact = 1 + real_model;
150    do_shift   = 0;
151    force_FFT  = 0;
152    force_poly = 0;
153    plot_type  = 1;
154  ##
155  ##  decode and check optional arguments
156  ##  end_numeric_args is boolean; becomes true at 1st string arg
157    end_numeric_args = 0;
158    for iarg = 1:length(varargin)
159      arg = varargin{iarg};
160      end_numeric_args = end_numeric_args || ischar(arg);
161      ## skip empty arguments
162      if ( isempty(arg) )
163        1;
164      ## numeric optional arguments must be first, cannot follow string args
165      ## N.B. older versions of matlab may not have the function "error" so
166      ## the user writes "function error(msg); disp(msg); end" and we need
167      ## a "return" here.
168      elseif ( ~ischar(arg) )
169        if ( end_numeric_args )
170          error( 'ar_psd: control arg must be string.' );
171        ##
172        ## first optional numeric arg is "freq"
173        elseif ( iarg == 1 )
174          user_freqs = isvector(arg) && length(arg)>1;
175          if ( ~isscalar(arg) && ~user_freqs )
176            error( 'ar_psd: arg 3 (freq) must be vector or scalar.' );
177          elseif ( ~user_freqs && ( ~isreal(arg) || ...
178                   fix(arg)~=arg || arg <= 2 || arg >= 1048576 ) )
179            error('ar_psd: arg 3 (freq) must be integer >=2, <=1048576' );
180          elseif ( user_freqs && ~isreal(arg) )
181            error( 'ar_psd: arg 3 (freq) vector must be real.' );
182          endif
183          freq = arg(:); # -> column vector
184        ##
185        ## second optional numeric arg is  "Fs" - sampling frequency
186        elseif ( iarg == 2 )
187          if ( ~isscalar(arg) || ~isreal(arg) || arg<=0 )
188            error( 'ar_psd: arg 4 (Fs) must be real positive scalar.' );
189          endif
190          Fs = arg;
191        ##
192        else
193          error( 'ar_psd: control arg must be string.' );
194        endif
195    ##
196    ## decode control-string arguments
197      elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
198        plot_type = 1;
199      elseif ( strcmp(arg,'semilogx') )
200        plot_type = 2;
201      elseif ( strcmp(arg,'semilogy') )
202        plot_type = 3;
203      elseif ( strcmp(arg,'loglog') )
204        plot_type = 4;
205      elseif ( strcmp(arg,'dB') )
206        plot_type = 5;
207      elseif ( strcmp(arg,'fft') )
208        force_FFT  = 1;
209        force_poly = 0;
210      elseif ( strcmp(arg,'poly') )
211        force_FFT  = 0;
212        force_poly = 1;
213      elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
214        pad_fact = 2;    # FFT zero-padding factor (pad FFT to double length)
215        do_shift = 0;
216      elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
217        pad_fact = 1;    # FFT zero-padding factor (do not pad)
218        do_shift = 0;
219      elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
220        pad_fact = 1;
221        do_shift = 1;
222      else
223        error( 'ar_psd: string arg: illegal value: %s', arg );
224      endif
225    endfor
226  ##  end of decoding and checking args
227  ##
228    if ( user_freqs )
229      ## user provides (column) vector of frequencies
230      if ( any(abs(freq)>Fs/2) )
231        error( 'ar_psd: arg 3 (freq) cannot exceed half sampling frequency.' );
232      elseif ( pad_fact==2 && any(freq<0) )
233        error( 'ar_psd: arg 3 (freq) must be positive in onesided spectrum' );
234      endif
235      freq_len = length(freq);
236      fft_len  = freq_len;
237      use_FFT  = 0;
238      do_shift = 0;
239    else
240      ## internally generated frequencies
241      freq_len = freq;
242      freq = (Fs/pad_fact/freq_len) * [0:freq_len-1]';
243      ## decide which method to use (poly or FFT)
244      is_power_of_2 = rem(log(freq_len),log(2))<10.*eps;
245      use_FFT = ( ~ force_poly && is_power_of_2 ) || force_FFT;
246      fft_len = freq_len * pad_fact;
247    endif
248    ##
249    ## calculate denominator of Equation 2.28, Kay and Marple, ref [1]Jr.:
250    len_coeffs = length(a);
251    if ( use_FFT )
252      ## FFT method
253      fft_out = fft( [ a(:); zeros(fft_len-len_coeffs,1) ] );
254    else
255      ## polynomial method
256      ## complex data on "half" frequency range needs -ve frequency values
257      if ( pad_fact==2 && ~real_model )
258        freq = [freq; -freq(freq_len:-1:1)];
259        fft_len = 2*freq_len;
260      endif
261      fft_out = polyval( a(len_coeffs:-1:1), exp( (-i*2*pi/Fs) * freq ) );
262    endif
263    ##
264    ## The power spectrum (PSD) is the scaled squared reciprocal of amplitude
265    ## of the FFT/polynomial. This is NOT the reciprocal of the periodogram.
266    ## The PSD is a continuous function of frequency.  For uniformly
267    ## distributed frequency values, the FFT algorithm might be the most
268    ## efficient way of calculating it.
269    ##
270    psd = ( v / Fs ) ./ ( fft_out .* conj(fft_out) );
271    ##
272    ## range='half' or 'onesided',
273    ##   add PSD at -ve frequencies to PSD at +ve frequencies
274    ## N.B. unlike periodogram, PSD at zero frequency _is_ doubled.
275    if ( pad_fact==2 )
276      freq = freq(1:freq_len);
277      if ( real_model )
278        ## real data, double the psd
279        psd = 2 * psd(1:freq_len);
280      elseif ( use_FFT )
281        ## complex data, FFT method, internally-generated frequencies
282        psd = psd(1:freq_len)+[psd(1); psd(fft_len:-1:freq_len+2)];
283      else
284        ## complex data, polynomial method
285        ##  user-defined and internally-generated frequencies
286        psd = psd(1:freq_len)+psd(fft_len:-1:freq_len+1);
287      endif
288    ##
289    ## range='shift'
290    ##   disabled for user-supplied frequencies
291    ##   Shift zero-frequency to the middle (pad_fact==1)
292    elseif ( do_shift )
293      len2 = fix((fft_len+1)/2);
294      psd  = [psd(len2+1:fft_len); psd(1:len2)];
295      freq = [freq(len2+1:fft_len)-Fs; freq(1:len2)];
296    endif
297    ##
298    ## Plot the spectrum if there are no return variables.
299    if ( nargout >= 2 )
300       varargout{1} = psd;
301       varargout{2} = freq;
302    elseif ( nargout == 1 )
303       varargout{1} = psd;
304    else
305      if ( plot_type == 1 )
306        plot(freq,psd);
307      elseif ( plot_type == 2 )
308        semilogx(freq,psd);
309      elseif ( plot_type == 3 )
310        semilogy(freq,psd);
311      elseif ( plot_type == 4 )
312        loglog(freq,psd);
313      elseif ( plot_type == 5 )
314        plot(freq,10*log10(psd));
315      endif
316    endif
317  endif
318
319endfunction
320