Burg and Welch power spectra

Hi everybody,

On Tue, Sep 26, 2006 at 11:55:11AM +0200, David Bateman wrote:
> I'd think that the best thing would be to replace the octave-forge
> functions with your new versions. Do you have a sourceforge account? If
> so then we should give you access to the CVS..

Following David's suggestion, I expect to commit replacements for arburg,
 __power, pburg and pwelch within the next week.
The replacements provide bugfixes, new features and a high level of
compatibility with matlab.  All of the replacement functions run in matlab.
Here are the usage strings and one-line comments:

[a,v,k] = arburg(x,poles,criterion)
%%               new arg, "criterion", can select the number of poles

[psd,f_out] = ar_psd(a,v,freq,Fs,range,method,plot_type)
%%               Power spectrum of AR model; replaces __power

[psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)
%%               Maximum-entropy power spectrum; wrapper for arburg and ar_psd

[psd,freq]  = pwelch(x,window,overlap,Nfft,Fs,
                       range,plot_type,detrend,sloppy)
[spectra,f] = pwelch(x,y,window,overlap,Nfft,Fs,
                       range,plot_type,detrend,sloppy,results)
[psd,ci,f]  = pwelch(x,window,overlap,Nfft,Fs,conf,
                       range,plot_type,detrend,sloppy)
%%               ci = confidence interval    
[spectra,ci,f] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
                          range,plot_type,detrend,sloppy,results)
%%               multiple "results" arg provides two-channel spectrum analyser

The upload will require changes to other functions in octave forge:

1) __power is also used by pyulear (power spectrum of Yule-Walker
   autoregressive (AR) model).  I plan to provide a new pyulear
   which takes advantage of the features in ar_psd.  However,
   I also need to fix a bug which gives an incorrectly scaled spectrum.

2) the current pwelch is a "backend" for the functions csd, tfe and cohere.
   These are compatible with Matlab R11 functions of the same names.
   Matlab R12 replaced csd, tfe and cohere with cpsd, tfestimate and mscohere,
   with different order of arguments and different defaults
   To maintain compatibility with matlab, the new pwelch recognises the
   following values of the global variable _pwelch_compatibility
     a) _pwelch_compatibility='R11-'
            produces compatibility with Matlab R11 and earlier versions
            [Pxx,f]=pwelch(x,Nfft,Fs,window,noverlap,conf,range,units);
     b) _pwelch_compatibility='R12+'
            produces compatibility with Matlab R12 and later versions
            [Pxx,f]=pwelch(x,window,noverlap,nfft,Fs,...);
     c) _pwelch_compatibility='psd.'
            gives the same default arguments as Matlab's spectrum.welch
            spectrum object and associated "psd" method.
   I have rewritten csd, tfe and cohere to use the new pwelch. Csd, tfe and
   cohere remain compatible with Matlab R11.
   While I am not entirely happy with using global variables, I have not
   been able to think of a better way... adding an extra argument won't
   work if you want old code to run unmodified.

3) I have written functions cpsd, tfestimate and mscohere. These provide
   the same capability as the corresponding Matlab sigproc functions.
   They are not currently in octaveforge.  Compatibility with
   Matlab varies according to the value of _pwelch_compatibility.

All comments and feedback will be gratefully received.

BTW. Does anyone know why Matlab spectrum functions have a default sampling
     frequency of 2*pi radians/second ?

Thanks
Peter Lanspeary