Re: Burg and Welch power spectra

Peter Vernon Lanspeary wrote:
> 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

In fact moving this mail to octave-dev makes even more sense...

I'd suggest just one modification the the above behavior. Get rid of the
global variable and use something like the following instead


function [Pxx, f] = pwelch (vargin)
  persistent _pwelch_compatibility = 'psd';

  if (nargin == 1)
    if (ischar(vargin{1})
      if (strcmp('R11-', vargin{1}))
        _pwlech_compatibility = 'R11-'
      elseif (strcmp('R12+', vargin{1}))
        _pwlech_compatibility = 'R12+'
      elseif (strcmp('psd', vargin{1}))
        _pwlech_compatibility = 'psd'
      else
        error("unrecognized compatibility string");
      endif
    else
      error("expected string argument")
    endif
    return
  endif

  # The rest of your function

Sorry, don't like global variables....

D.
if (strcmp('R12+', vargin{1}))