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