Reply on Burg and Welch power spectra

[Peter Lanspeary]

> Date: Thu, 15 Jun 2006 18:43:00 +0200

>>1) welch_psd.m   -- (Welch 1967) averaged-periodogram method
>>2) burg_filter.m -- Burg (1968) lattice-filter algorithm
>>3) mem_psd.m     -- spectral estimation from Burg coefficients
>
> Can you please compare your implementation of the pwelch to the
> octave-forge version and to Matlab 2006a for compatiability? What would
> it take to create a fully Matlab compatiable pburg function? Once that
> is done these should be included in the signal processing toolbox of
> octave-forge...
> Regards
> David
>--
>David Bateman                                David.Bateman@...
>Motorola Labs - Paris                        [phone] (Ph)
>Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
>91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)

Hi David and John (jwe),

Initial versions of welch_psd.m, burg_filter.m and mem_psd.m are inadequate
and have been rewritten.  The last two are renamed burg_model.m and ar_psd.m
respectively, and have been upgraded to work with complex data.  I have added
burg_psd.m which is compatible with the matlab pburg.m

Octave-forge already has ".m" functions with the same names as
Matlab's pwelch, pburg and arburg, in each case for the purpose of performing
essentially the same job, but in detail there are many differences and
some bugs.

I have compared arguments and behaviour of
   a) welch_psd with pwelch in octave and matlab,
   b) burg_model with arburg in octave and matlab,
   c) burg_psd with pburg in octave and matlab.

Octave/matlab versions and O/S used in these tests are
   a) Octave 2.1.69 + octave-forge-2006.07.09 on Debian 3.1r1
      with sourceforge.net/tracker bugfix #1530283   for arburg.m
      with sourceforge.net/tracker bugfix #1530309   for __power.m
   b) Matlab 7.0.1.24704(R14)SP1 + sigproc toolbox v6.2.1, on MS-windows XP
      (this is the most recent Matlab version available to me)

Details of the comparisons and tests, and source code are attached.

Thanks,
Peter Lanspeary

COMPARISON OF WELCH and BURG POWER-SPECTRUM FUNCTIONS
P. V. Lanspeary

Contents
--------
Section 1. Comparison of arguments. Summary.
Section 2. Comparison of behaviours. Summary.
Section 3. pburg: Full compatibility with Matlab ?
Section 4. Comparison of arguments. Details.
Section 5. Comparison of behaviours. Details.


Section 1. Comparison of arguments. Summary.
--------------------------------------------
I have compared the arguments of welch_psd.m, burg_model.m and burg_psd.m
with equivalent functions in octave-forge and Matlab-sigproc.
This is a summary. See section (4) for details.

(a) welch_psd with pwelch (octave) and pwelch (matlab):
  The behaviour of welch_psd and pwelch is complicated, so details of the
  comparison are messy and complicated. In summary, the news is not
  good.  Overall functionality is similar (not the same). Only the
  output args and first input arg are compatible.  

(b) burg_model with arburg (octave) and arburg (matlab):
  Arguments are compatible except that burg_model has an extra optional
  input argument.  If each function is called with the same 2 arguments,
  behaviour should (according to the documentation) be the same.

(c) burg_psd with pburg (octave) and pburg (matlab):
  Arguments are compatible except that Octave's pburg has one an extra
  optional input argument and burg_psd has 3 extra optional input arguments.


Section 2. Comparison of behaviours. Summary.
---------------------------------------------
I have run _two_ test cases on the "Welch" functions and two more on
the "Burg" functions.  One test case has a real signal (data); the
other has a complex signal.  This is a summary.  See section (5) for
details.

(a) pwelch (octave), pwelch (matlab), welch_psd:
  Arguments are incompatible; so to do the same thing with each function
  requires different arguments in each case.  If each function is
  called with the same arguments, behaviour (except in special cases) is
  different.  Producing the same behaviour requires care. Where documentation
  requirees the same behaviour, output is the same to an accuracy of no worse
  than 10x machine precision (eps*10).

(b) arburg (octave), arburg (matlab), burg_model:
  For the real-data test case, results are exactly the same --- except that
  reflection coefficients (k) from octave arburg have the wrong sign.
  For complex data, the Matlab arburg gives the same (correct) answer as
  burg_model, but the octave arburg does not give the correct answer.

(c) pburg (octave), pburg (matlab), burg_psd:
  For 'twosided' spectra Matlab pburg gives the same answers as burg_psd, but
  for 'onesided' spectra, the first and last values in the spectrum are too
  small by a factor of 2.
  The octave pburg does not work for complex data because it calls the octave
  arburg.  The octave __power.m works correctly for both real and complex
  data.


Section 3. pburg: Full compatibility with Matlab ?
--------------------------------------------------

>> What would it take to create a fully Matlab compatiable pburg function?
Arguments are fully compatible, but with optional extra arguments.
Compatible behaviour is harder to achieve, and I am not sure that this is
desirable. For example, the Matlab pburg
a) assumes units of frequency are radians/sec unless sampling
   frequency, Fs, is specified.  If Fs is specified, the units are Hertz.
   Octave pburg uses only Hertz.
b) returns only nfft/2+1 (frequency,spectral-density) coordinates if the
   spectrum is onesided.  Octave pburg returns nfft values.
c) plots the spectral density in dB rather than correct physical units;
   the graph should be a power spectrum of a physical quantity rather
   than the response of a filter.
d) has a error in onesided spectra (the first and last values of spectral
   density are too small by a factor of 2).
A completely fixed Octave pburg would seem to have more appropriate behaviour.


Section 4. Comparison of arguments. Details.
--------------------------------------------

(a) pwelch (octave), pwelch (matlab), welch_psd:
  The number of arguments is large and behaviour of the functions is a bit
  complicated ... for the sake of brevity this comparison/description
  therefore contains some minor approximations. First we write the first line
  of each function, adjusting argument names to make them more consistent::

  1) octave [Pxx,Pci,f]=pwelch(x,nfft,Fs,window,overlap,ci,range,units,trend)
  2) matlab [Pxx,f]    =pwelch(x,window,overlap,nfft,Fs,range)
  3)        [Pxx,f]=welch_psd(x,window,Fs,overlap,pad,range,units,trend,sloppy)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave pwelch     |   matlab pwelch    |   welch_psd        |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  |  Pxx     | spectral density   | spectral density   | spectral density   |
  |          | units specified by | ["x" units]/Hertz  | ["x" units]/Hertz  |
  |          | "units" arg        | or per (rad/s)     |                    |
  +----------+--------------------+--------------------+--------------------+
  | Pci      | confidence interval| not available      | not available      |
  +----------+--------------------+--------------------+--------------------+
  | f        | frequency (Hertz)  | frequency (Hertz  | frequency          |
  |          |                    | or rad/s)          | (Hertz if given Fs)|
  |          |                    |                    | (rad/s without Fs) |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | data (signal)                                                |
  +----------+--------------------+--------------------+--------------------+
  | window   | windowing vector, or the length of data segment              |
  |          +--------------------+--------------------+--------------------+
  |          | default = Hann     | default = Hamming  | default = Hann     |
  |          | default_length=nfft| default_length=    | default_length=    |
  |          |                    | length(x)/8        | sqrt(length(x))    |
  +----------+--------------------+--------------------+--------------------+
  | Fs       | sampling frequency | sampling frequency | sampling frequency |
  |          | default=2 Hz       | default= 2*pi Hz   | default=1 Hz       |
  +----------+--------------------+--------------------+--------------------+
  | nfft     | length of FFT      | length of FFT or   | not available;     |
  |          | default = 256      | list of frequencies| uses window length |
  |          |                    | default=256 or     | + "pad" as nfft    |
  |          |                    | rounded window len |                    |
  +----------+--------------------+--------------------+--------------------+
  | overlap  | number of overlap  | number of overlap  | overlap fraction   |
  |          | samples, default=  | samples, default=  | default=0.5, i.e.  |
  |          | length(window)/2   | length(window)/2   | 50%                |
  +----------+--------------------+--------------------+--------------------+
  | pad      | not available      | not available      | length of zero     |
  |          |                    |                    | padding on FFT     |
  |          |=nfft-length(window)|=nfft-length(window)| default=0          |
  +----------+--------------------+--------------------+--------------------+
  | ci       | confidence interval| not available      | not available      |
  +----------+--------------------+--------------------+--------------------+
  | range    | 'half'             | 'onesided'         | 'half', 'onesided' |
  |          | 'whole'            | 'twosided'         | 'whole', 'twosided'|
  |          |                    |                    | 'shift'            |
  +----------+--------------------+--------------------+--------------------+
  | units    | plot scaling       | not available      | 'plot', 'semilogx' |
  |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
  |          |                    |                    | 'loglog', 'squared'|
  +----------+--------------------+--------------------+--------------------+
  | trend    | remove trend       | not available      |remove trend, 'long'|
  |          | 'mean', 'linear'   | -- does not        |'short','no-detrend'|
  |          | default=none       |    remove mean     | default='long'     |
  +----------+--------------------+--------------------+--------------------+
  | sloppy   | not available      | not available      | round nfft to 2^n  |
  +----------+--------------------+--------------------+--------------------+
    The above table does not describe the special cases where the octave pwelch
    and welch_psd are used for calculating cross-power spectra, coherence
    and transfer functions.

(b) arburg (octave), arburg (matlab), burg_model:
  (1) octave  [a,v,k]=arburg(x,p)
  (2) matlab  [a,v,k]=arburg(x,p)
  (3)         [a,v,k]=burg_model(x,p,stop)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave arburg     |   matlab arburg    |   burg_model       |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | a        | autoregression coefficients                                  |
  +----------+--------------------+--------------------+--------------------+
  | v        | mean square residual (noise)                                 |
  +----------+--------------------+--------------------+--------------------+
  | k        | reflection coefficients                                      |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | real data          | real data          | complex data       |
  +----------+--------------------+--------------------+--------------------+
  | p        | number of poles (required arg)                               |
  +----------+--------------------+--------------------+--------------------+
  | stop     | not available      | not available      | stopping criterion |
  |          |                    |                    | 'FPE','AIC'        |
  |          |                    |                    | default='none'     |
  +----------+--------------------+--------------------+--------------------+

(c) pburg (octave), pburg (matlab), burg_psd:
  (1) octave  [Pxx,f]=pburg(x,p,nfft,Fs,range,units)
  (2) matlab  [Pxx,f]=pburg(x,p,nfft,Fs,range)
  (3)         [Pxx,f]=burg_psd(x,p,nfft,Fs,range,units,method,stop)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave pburg      |   matlab pburg     |   burg_psd         |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  |  Pxx     | spectral density   | spectral density   | spectral density   |
  |          | ["x" units]/Hertzy | ["x" units]/Hertz  | ["x" units]/Hertz  |
  |          |                    | or per (rad/s)     |                    |
  +----------+--------------------+--------------------+--------------------+
  | f        | frequency (Hertz)  | frequency (Hertz   | frequency          |
  |          |                    | or rad/s)          | (Hertz if given Fs)|
  |          |                    |                    | (rad/s without Fs) |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | real data          | real data          | complex data       |
  +----------+--------------------+--------------------+--------------------+
  | p        | number of poles (required arg)                               |
  +----------+--------------------+--------------------+--------------------+
  | nfft     | number of frequncy | length of FFT or   | number of frequncy |
  |          | values             | list of frequencies| values or list of  |
  |          |                    |                    | frequencies        |
  |          | default = 256      | default=256        | default=256        |
  +----------+--------------------+--------------------+--------------------+
  | Fs       | sampling frequency | sampling frequency | sampling frequency |
  |          | default=2 Hz       | default= 2*pi Hz   | default=1 Hz       |
  +----------+--------------------+--------------------+--------------------+
  | range    | 'half'             | 'onesided'         | 'half', 'onesided' |
  |          | 'whole'            | 'twosided'         | 'whole', 'twosided'|
  |          |                    |                    | 'shift'            |
  +----------+--------------------+--------------------+--------------------+
  | units    | plot scaling       | not available      | 'plot', 'semilogx' |
  |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
  |          |                    |                    | 'loglog', 'squared'|
  +----------+--------------------+--------------------+--------------------+
  | method   | not available      | not available      | 'FFT', 'poly'      |
  +----------+--------------------+--------------------+--------------------+
  | stop     | not available      | not available      | stopping criterion |
  |          |                    |                    | 'FPE','AIC'        |
  |          |                    |                    | default='none'     |
  +----------+--------------------+--------------------+--------------------+


Section 5. Comparison of behaviours. Details.
---------------------------------------------
The test platforms are
   a) Octave 2.1.69 + octave-forge-2006.07.09 on Debian 3.1r1
      with sourceforge.net/tracker bugfix #1530283   for arburg.m
      with sourceforge.net/tracker bugfix #1530309   for __power.m
   b) Matlab 7.0.1.24704(R14)SP1 + sigproc toolbox v6.2.1, on MS-windows XP
      (this is the most recent Matlab version available to me)

(a) pwelch (octave), pwelch (matlab), welch_psd:
  The test script is

     format long
     rand('seed',2038014164);
     a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
     white = rand(1,16384);
     signal = detrend(filter(0.70181,a,white));
     skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
     %% REAL-DATA TEST
    [psd1,f1]=welch_psd(signal,hann(256),25,0.5,0,'whole','plot','no-detrend');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd2,f2]=pwelch(signal,256,25,hann(256),128,'whole','squared','none');
     else
       disp( 'MATLAB' )
       [psd2,f2]=pwelch(signal,hann(256),256*0.5,256,25,'twosided');
       end
     sum_psd1=sum(psd1)
     sum_psd2=sum(psd2)
     maxdiffr=max(abs(psd1./reshape(psd2,1,[])-1))  %% max relative difference
     %% COMPLEX-DATA TEST
    [psd3,f3]=welch_psd(skewed,hann(256),25,0.5,0,'whole','plot','no-detrend');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd4,f4]=pwelch(skewed,256,25,hann(256),128,'whole','squared','none');
     else
       disp( 'MATLAB' )
       [psd4,f4]=pwelch(skewed,hann(256),256*0.5,256,25,'twosided');
       end
     sum_psd3=sum(psd3)
     sum_psd4=sum(psd4)
     maxdiffc=max(abs(psd3./psd4'-1))  %% maximum relative difference
     eps

  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
  | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | sum_psd1 | 3.62735601946465      | 3.61040879075053          |
  | sum_psd2 | 3.62735601946465      | 3.61040879075053          |
  | maxdiffr | 1.33226762955019e-15  | 1.110223024625157e-015    |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | sum_psd3 | 3.62735601946465      | 3.61040879075054          |
  | sum_psd4 | 3.62735601946465      | 3.61040879075054          |
  | maxdiffc | 1.33226762955019e-15  | 1.332267629550188e-015    |
  +----------+-----------------------+---------------------------+
    We should have sum_psd1==sum_psd2, maxdiffr==0, sum_psd3==sum_psd4,
    maxdiffc==0.


(b) arburg (octave), arburg (matlab), burg_model:
  The additional test script (carrying on from previous test) is

     %% REAL-DATA TEST
     [ar1,vr1,kr1]=burg_model(signal,5);
     [ar2,vr2,kr2]=arburg(signal,5);
     maxdifar=max(abs(ar2./ar1-1))  %% maximum relative difference
     vr1/vr2-1
     maxdifkr=max(abs(reshape(kr2,1,[])./kr1-1)) %% kr2 has inconsistent shape
     %% COMPLEX-DATA TEST
     [ac1,vc1,kc1]=burg_model(skewed,5);
     [ac2,vc2,kc2]=arburg(skewed,5);
     maxdifac=max(abs(ac2./ac1-1))  %% maximum relative difference
     vc2/vc1-1
     maxdifkc=max(abs(reshape(kc2,1,[])./kc1-1))

  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | burg_model vs. arburg | burg_model vs. arburg     |
  | eps      | 2.22044604925031e-16  |  2.220446049250313e-016   |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | maxdifar | 8.88178419700125e-16  |  0                        |
  | vr1/vr2-1| -5.55111512312578e-16 |  0                        |
  | maxdifkr | 2.00000000000000      |  0                        |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | maxdifac | 2.15561046016939      |  0                        |
  | vc2/vc1-1| 7.679016 - 0.0026252i |  0                        |
  | maxdifkc | 1.68536747392439      |  5.590172498962309e-017   |
  +----------+-----------------------+---------------------------+
     We should have maxdifar==0, vr1/vr2-1==0, maxdifkr==0, maxdifac==0,
     vc2/vc1-1==0, maxdifkc==0.
     ... Octave arburg gives the correct autoregression coefficients and
     variance (ar1,vr1) for real data but the reflection coefficients have
     the wrong sign; for complex data, the octave arburg gives the wrong
     answers.  Matlab pburg and burg_model give exactly the same answer
     for both real and complex data.
     Now compare __power with ar_psd using coefficients (ac1,vc1)
     calculated by burg_model  ...

     [psd3,f3]=ar_psd(ac1,vc1,256,25,'whole','squared');
     [psd4,f4]=__power(sqrt(vc1),ac1,256,25,'whole','squared');
     sum_psd3=sum(psd3)
     %%     sum_psd3 = 3.63449149426036e+00 + 4.89110705062523e-19i
     sum_psd4=sum(psd4)
     %%     sum_psd4 = 3.63449149426036
     maxdiffc=max(abs(psd3./psd4-1))
     %%     maxdiffc =  3.33066934574812e-16

     ... we conclude that __power works OK for complex filter coefficients

(b) pburg (octave), pburg (matlab), burg_psd:
  The additional test script (carrying on from previous test) is

     %% REAL-DATA TEST
     [psd5,f5]=burg_psd(signal,5,256,25,'whole');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd6,f6]=pburg(signal,5,256,25,'whole');
     else
       disp( 'MATLAB' )
       [psd6,f6]=pburg(signal,5,256,25,'twosided');
       end
     sum_psd5=sum(psd5)
     sum_psd6=sum(psd6)
     maxdiffr=max(abs(psd5./reshape(psd6,1,[])-1))  %% max relative difference
     %% COMPLEX-DATA TEST
     [psd7,f7]=burg_psd(skewed,5,256,25,'whole');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd8,f8]=pburg(skewed,5,256,25,'whole');
     else
       disp( 'MATLAB' )
       [psd8,f8]=pburg(skewed,5,256,25,'twosided');
       end
     sum_psd7=sum(psd7)
     sum_psd8=sum(psd8)
     maxdiffc=max(abs(psd7./reshape(psd8,1,[])-1))  %% max relative difference
 
  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
  | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | sum_psd5 | 3.63449149426037      | 3.62433951286238          |
  | sum_psd6 | 3.63449149426037      | 3.62433951286238          |
  | maxdiffr | 1.24344978758018e-14  | 5.551115123125783e-016    |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | sum_psd7 | 3.634491494 +4.9e-19i | 3.62433951286238          |
  | sum_psd8 | 3.64552153406807      | 3.62433951286238          |
  | maxdiffc | 6.15692319331572      | 6.661338147750939e-016    |
  +----------+-----------------------+---------------------------+
     We should have sum_psd5==sum_psd6, maxdiffr==0, sum_psd7==sum_psd8,
     maxdiffc==0.
     Octave pburg doesn't work for complex data because it calls  arburg.
     Twosided spectra from matlab pburg are OK, but for "onesided" spectrum
      [psd8,f8]=pburg(skewed,5,256,25,'onesided');
    the first and last values in the spectrum from matlab are too small by
    a factor of 2.

%% Copyright (C) 2006 Peter V. Lanspeary
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
%% as published by the Free Software Foundation; either version 2,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software Foundation,
%% Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% Usage:
%%   [psd,freq_out] = ar_psd(ar_coeffs,residual,freq,sample_f,
%%                           range,method,plot_type)
%%
%%  Calculate the power spectrum of the autoregressive model
%%
%%                                M
%%  x(n) = sqrt(residual).e(n) + SUM ar_coeffs(k).x(n-k)
%%                               k=1
%%  where x(n) is the output of the model and e(n) is white noise.
%%  This function is intended for use with
%%    [ar_coeffs,residual,lattice] = burg_model(data,poles,stop_crit)
%%    [ar_coeffs,residual,lattice] = arburg(data,poles)
%%  which use the Burg (1968) method to calculate a "maximum entropy"
%%  autoregressive model of "data".  This function runs on octave and matlab.
%%  
%%  If the "freq" argument is a vector (of frequencies) the spectrum is
%%  calculated using the polynomial method and the "method" argument is
%%  ignored.  For scalar "freq", an integer power of 2, or "method='FFT'",
%%  causes the spectrum to be calculated by FFT.  Otherwise, the spectrum
%%  is calculated as a polynomial.  It may be computationally more
%%  efficient to use the FFT method if length of the model is not much
%%  smaller than the number of frequency values. The spectrum is scaled so
%%  that spectral energy (area under spectrum) is the same as the
%%  time-domain energy (mean square of the signal).
%%
%% ARGUMENTS:
%%     All but the first two arguments are optional and may be empty.
%%
%%   ar_coeffs %% [vector] list of M=(order+1) autoregressive model
%%             %%      coefficients.  The first element of "ar_coeffs" is the
%%             %%       zero-lag coefficient, which always has a value of 1.
%%
%%   residual  %% [real scalar] square of the moving-average coefficient of
%%             %%               the AR model.
%%
%%   freq      %% [real vector] frequencies at which power spectral density
%%             %%               is calculated
%%             %% [integer scalar] number of uniformly distributed frequency
%%             %%          values at which spectral density is calculated.
%%             %%          [default=256]
%%
%%   sample_f  %% [real scalar] sampling frequency (Hertz) [default=1]
%%
%% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
%%   Control-string arguments can be in any order after the other arguments.
%%
%%   range     %% 'half',  'onesided' : frequency range of the spectrum is
%%             %%       zero up to but not including sample_f/2.  Power from
%%             %%       negative frequencies is added to the positive side of
%%             %%       the spectrum.
%%             %% 'whole', 'twosided' : frequency range of the spectrum is
%%             %%       -sample_f/2 to sample_f/2, with negative frequencies
%%             %%       stored in "wrap around" order after the positive
%%             %%       frequencies; e.g. frequencies for a 10-point 'twosided'
%%             %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
%%             %% 'shift' : same as 'whole' but with the first half of the
%%             %%       spectrum swapped with second half to put the zero-
%%             %%       frequency value in the middle. (See "help fftshift".)
%%             %%       If "freq" is vector, 'shift' is ignored.
%%             %% If model coefficients "ar_coeffs" and "residual" are real,
%%             %% default range is 'half', otherwise default range is 'whole'.
%%
%%   method    %% 'FFT':  use FFT to calculate power spectrum.
%%             %% 'poly': calculate power spectrum as a polynomial of 1/z
%%             %% N.B. this argument is ignored if the "freq" argument is a
%%             %%      vector.  The default is 'poly' unless the "freq"
%%             %%      argument is an integer power of 2.
%%  
%%   plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
%%             %%       specifies the type of plot.  The default is 'plot',
%%             %%       which means linear-linear axes. 'squared' is the same
%%             %%       as 'plot'. 'db' plots "10*log10(psd)".
%%
%% RETURNED VALUES:
%%     If return values are not required by the caller, the spectrum
%%     is plotted and nothing is returned.
%%   psd       %% [real vector] estimate of power-spectral density
%%   freq_out  %% [real vector] frequency values
%%
%%
%% REFERENCES
%% John Parker Burg (1968):
%%   "A new analysis technique for time series data",
%%   NATO advanced study Institute on Signal Processing with Emphasis on
%%   Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.
%%
%% William H. Press and Saul A. Teukolsky and William T. Vetterling and
%%               Brian P. Flannery",
%% "Numerical recipes in C, The art of scientific computing", 2nd edition,
%%    Cambridge University Press, 2002 --- Section 13.7.
%% N.B. The algorithm in Press et al. expects prediction-filter coefficients.
%%      ar_psd requires the whitening-filter coefficients.
%%

function [varargout]=ar_psd(ar_coeffs,residual,varargin)
%%
%% Check fixed arguments
if ( nargin < 2 )
  error( 'ar_psd: error: needs >=2 args. Use "help ar_psd"');
elseif ( ~isvector(ar_coeffs) || length(ar_coeffs)<2 )
  error( 'ar_psd:error: arg 1 (ar_coeffs) must be vector, length>=2.' );
elseif ( ~isscalar(residual) )
  error( 'ar_psd: error: arg 2 (residual) must be real scalar >0' );
else
  real_model = isreal(ar_coeffs);
%%
%%  default values for optional areguments
  freq = 256;
  user_freqs = 0;    %% boolean: true for user-specified frequencies
  sample_f   = 1.0;
  %%  FFT padding factor (is also frequency range divisor): 1=whole, 2=half.
  pad_fact = 1 + real_model;
  do_shift   = 0;
  force_FFT  = 0;
  force_poly = 0;
  plot_type  = 1;
%%
%%  decode and check optional arguments
%%  end_numeric_args is boolean; becomes true at 1st string arg
  end_numeric_args = 0;
  for iarg = 1:length(varargin)
    arg = varargin{iarg};
    end_numeric_args = end_numeric_args || ischar(arg);
    %% skip empty arguments
    if ( isempty(arg) )
      1;
    %% numeric optional arguments must be first, cannot follow string args
    elseif ( ~ischar(arg) )
      if ( end_numeric_args )
        error( 'ar_psd: error: control arg must be string' );
        return;
      %%
      %% first optional numeric arg is "freq"
      elseif ( iarg == 1 )
        user_freqs = isvector(arg) && length(arg)>1;
        if ( ~isscalar(arg) && ~user_freqs )
          error( 'ar_psd: error: arg 3 (freq) must be vector or scalar.');
          return;
        elseif ( ~user_freqs && ( ~isreal(arg) || ...
                 fix(arg)~=arg || arg <= 2 || arg >= 1048576 ) )
          error('ar_psd: error: arg 3 (freq) is not integer >=2, <=1048576');
          return;
        elseif ( user_freqs && ~isreal(arg) )
          error( 'ar_psd: error: arg 3 (freq) vector must be real' );
          return;
          end
        freq = arg;
      %%
      %% second optional numeric arg is  "sample_f" - sampling frequency
      elseif ( iarg == 2 )
        if ( ~isscalar(arg) || ~isreal(arg) || arg<=0 )
          error( 'ar_psd: error: arg 4 (sample_f) must be real scalar >0.' );
          return;
          end
        sample_f = arg;
      %%
      else
        error( 'ar_psd: error: control arg must be string' );
        return;
        end
  %%
  %% decode control-string arguments
    elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
      plot_type = 1;
    elseif ( strcmp(arg,'semilogx') )
      plot_type = 2;
    elseif ( strcmp(arg,'semilogy') )
      plot_type = 3;
    elseif ( strcmp(arg,'loglog') )
      plot_type = 4;
    elseif ( strcmp(arg,'db') )
      plot_type = 5;
    elseif ( strcmp(arg,'FFT') )
        force_FFT  = 1;
        force_poly = 0;
    elseif ( strcmp(arg,'poly') )
        force_FFT  = 0;
        force_poly = 1;
    elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
        pad_fact = 2;
        do_shift = 0;
    elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
      pad_fact = 1;
      do_shift = 0;
    elseif ( strcmp(arg,'shift') )
      pad_fact = 1;
      do_shift = 1;
    else
      error( strcat( 'ar_psd: error: string arg: illegal value: ', arg ) );
      return;
      end
    end
%%  end of decoding and checking args
%%
%% define the frequencies
  if ( user_freqs )
    %% user-supplied frequencies
    if ( any(abs(freq))>sample_f/2 )
      error( 'ar_psd: error: arg 3 (freq) must not exceed sampling freq/2' );
      return;
      end
    len_freq = length(freq);
    use_FFT = 0;
    do_shift = 0;
    %% "half" frequency range with complex data needs -ve frequency values
    if ( pad_fact==2 && ~real_model )
      freq = [freq -freq(len_freq:-1:1)];
      end
    fft_len = length(freq);
  else
    %% internally generated frequencies
    len_freq = freq;
    freq = [0:len_freq-1] * sample_f /pad_fact / len_freq;
    %% decide which method to use (poly or FFT)
    is_power_of_2 = rem(log(len_freq),log(2))<10.*eps;
    use_FFT = ( ~ force_poly && is_power_of_2 ) || force_FFT;
    fft_len = len_freq*pad_fact;
    end
%%
%% calculate FFT or polynomial
  len_coeffs = length(ar_coeffs);
  if ( use_FFT )
    fft_out = fft( [ ar_coeffs zeros(1,fft_len-len_coeffs) ] );
  else %% use polynomial
    two_pi_i = (0+i) * 2 * pi / sample_f;
    fft_out = polyval( ar_coeffs(len_coeffs:-1:1), exp( two_pi_i * freq ) );
    end
%%
%% The power spectrum (PSD) is the scaled squared reciprocal of amplitude of
%% the FFT/polynomial. This is NOT the reciprocal of the periodogram.
%% The PSD is a continuous function of frequency.  For carefully chosen
%% frequency values, the FFT algorithm might be the most efficient way of
%%    calculating it.
%%
  psd = (residual/sample_f) ./ ( fft_out .* conj(fft_out) );
%%
%% for range='half', add PSD at -ve frequencies to PSD at +ve frequencies
%% N.B. unlike periodogram, PSD at zero frequency is doubled.
  if ( pad_fact==2 )
    freq = freq(1:len_freq);
    if ( real_model )
      psd = 2 * psd(1:len_freq); % might be 'poly' or 'fft'
    elseif ( user_freq ) % complex model
      psd = psd(1:len_freq)+psd(fft_len:-1:len_freq+1);
    else  % internally-generated frequencies, complex model
      psd = psd(1:len_freq)+[psd(1) psd(fft_len:-1:len_freq+2)];
      end
%%
%% range='shift': Shift zero-frequency to the middle (pad_fact==1)
  elseif ( do_shift )
    len2 = fix((fft_len+1)/2);
    psd  = [psd(len2+1:fft_len), psd(1:len2)];
    freq = [freq(len2+1:fft_len)-sample_f, freq(1:len2)];
    end
%%
%%
%% Plot the spectrum if there are no return variables.
  if ( nargout >= 2 )
     varargout{1} = psd;
     varargout{2} = freq;
  elseif ( nargout == 1 )
     varargout{1} = psd;
  else
    if ( plot_type == 1 )
      plot(freq,psd);
    elseif ( plot_type == 2 )
      semilogx(freq,psd);
    elseif ( plot_type == 3 )
      semilogy(freq,psd);
    elseif ( plot_type == 4 )
      loglog(freq,psd);
    elseif ( plot_type == 5 )
      plot(freq,10*log10(psd));
      end
    end
  end
end

%% Copyright (C) 2006 Peter V. Lanspeary
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
%% as published by the Free Software Foundation; either version 2,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software Foundation,
%% Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% [coeffs,residual,lattice] = burg_model(data,poles,stop_crit)
%%
%% Calculate coefficients of an autoregressive (AR) model of "data" using the
%% lattice-filter method of Burg (1968).  The inverse of the model is a
%% moving-average filter which reduces the "data" to white noise.
%% The algorithm is adapted for complex data (Kay & Marple, 1981).
%% The power spectrum of the AR model is an estimate of the maximum
%% entropy power spectrum of the data.  The function "ar_psd" calculates the
%% power spectrum of the AR model.
%%
%% burg_model does not remove the mean from the data.  You should remove the
%% mean from the data if you want a power spectrum.  A non-zero mean can
%% produce large errors in a power-spectrum estimate.  See "help detrend".
%%
%% ARGUMENTS:
%%   data      %% [vector] sampled data
%%
%%   poles     %% [integer scalar] required number of poles of AR model
%%
%%   stop_crit %% [optional string arg]
%%             %% 'FPE' -- apply "final prediction error" criterion
%%             %% 'AIC' -- apply "Akaike information criterion"
%%             %%    (Kay & Marple, 1981) to limit the number of poles so that
%%             %%    spurious poles are not added when the whitened data has
%%             %%    no more information in it. The default is to NOT apply
%%             %%    the AIC or FPE.
%%
%% RETURNED VALUES:
%%   coeffs    %% [real vector] list of M=(poles+1) moving-average model
%%             %%               coefficients; for data input x(n) and
%%             %%               white noise output w(n), the model is
%%             %%                                     M
%%             %%       x(n) = sqrt(residual).e(n) + SUM ar_coeffs(k).x(n-k)
%%             %%                                    k=1
%%
%%   residual  %% mean square of residual (white) noise from the whitening
%%             %% operation of the Burg lattice filter.
%%
%%   lattice   %% reflection coefficients defining the lattice-filter
%%             %% embodiment of the model
%%
%% REFERENCES
%% John Parker Burg (1968)
%%   "A new analysis technique for time series data",
%%   NATO advanced study Institute on Signal Processing with Emphasis on
%%   Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.
%%
%% Steven M. Kay and Stanley Lawrence Marple Jr. (1981):
%%   "Spectrum analysis -- a modern perspective",
%%   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
%%
%% William H. Press and Saul A. Teukolsky and William T. Vetterling and
%%               Brian P. Flannery",
%% "Numerical recipes in C, The art of scientific computing", 2nd edition,
%%    Cambridge University Press, 2002 --- Section 13.7.

function [varargout] = burg_model( data, poles, stop_crit )
%%
%% Check arguments
if ( nargin < 2 )
  error( 'burg_model(data,poles): Need at 2 least args.' );
elseif ( ~isvector(data) || length(data) < 3 )
  error( 'burg_model:error: arg 1 (data) must be vector of length >3.' );
elseif ( ~isscalar(poles) || ~isreal(poles) || fix(poles)~=poles || poles<=0.5)
  error( 'burg_model:error: arg 2 (poles) must be positive integer.' );
elseif ( floor(poles+0.5) > length(data)-2 )
  error( 'burg_model:error: arg 2 (poles)+2 must be less than data length' );
elseif ( nargin>2 && ~isempty(stop_crit) && ~ischar(stop_crit) )
  error( 'burg_model:error: arg 3 (stop_crit) must be string' );
else
%%
%%
%% Decide if FPE or AIC criteria will be applied to stop before
%% getting to the specified number of model poles.
  if ( nargin > 2 && ~isempty(stop_crit) )
    use_FPE = strcmp(stop_crit,'FPE');
    use_AIC = strcmp(stop_crit,'AIC');
  else
    use_FPE = 0;
    use_AIC = 0;
    end
%%
%% Storage of forward and backward prediction errors is a little tricky.
%% Because the forward error f(i) is always combined with the lagged
%% backward error b(i-1), f(1) and b(n) are never used, and therefore are
%% never stored.  Not storing unused data makes the calculation of the
%% reflection coefficient look much cleaner :)
%% N.B. {initial residual} = {error for zero-order model} =
%%      {zero-lag autocorrelation} =  E(x*conj(x)) = x*x'
%% Also there is no sumsq or meansq function in Matlab (for real data)
  N = length(data);
  lattice = [];
  if ( size(data,1) > 1 ) %% if data is row vector
    forw_err = data(2:N).';
    back_err = data(1:N-1).';
    residual = data' * data / N;
  else %% if data is column vector
    forw_err = data(2:N);
    back_err = data(1:N-1);
    residual = data * data' / N;
    end
  %% new_criterion/old_criterion are either FPE or AIC
  new_criterion = abs(residual);
  old_criterion = 2 * new_criterion;
  for k = 1:poles
    %%
    %% reflection_coeff = -2* E(f(i)*conj(b(i-1))) / ( E(f(i)^2)+E(b(i-1)^2) )
    refl_coeff= -2 * forw_err * back_err' / ...
                ( forw_err * forw_err' + back_err * back_err');
    %%  Levinson-Durbin recursion for residual
    new_residual = residual * ( 1.0 - refl_coeff * conj(refl_coeff) );
    if ( k > 1 )
      %%
      %% Apply the FPE or AIC criterion and stop if the FPE or AIC is
      %% increasing rather than decreasing. For complex data, use abs(residual)
      %% Do it before we update the old model "coeffs" and "residual".
      if ( use_FPE )
        old_criterion = new_criterion;
        new_criterion = abs(new_residual) * ( N + k + 1 ) / ( N - k - 1 );
        if ( new_criterion > old_criterion )
           break;
           end
      elseif ( use_AIC )
        old_criterion = new_criterion;
        new_criterion = log(abs(new_residual)) + 2 * ( k + 1 ) / N;
        if ( new_criterion > old_criterion )
          break;
           end
        end
      %% Update model "coeffs" and "residual".
      %% Use Levinson-Durbin recursion formula (for complex data).
      coeffs = [ prev_coeffs + refl_coeff .* conj(prev_coeffs(k-1:-1:1)),...
                 refl_coeff ];
    else %% if( k==1 )
      coeffs = refl_coeff;
      end
    lattice = [ lattice, refl_coeff ];
    residual = new_residual;
    if ( k < poles )
      prev_coeffs = coeffs;
      %%  calculate new prediction errors (by recursion):
      %%  f(p,i) = f(p-1,i)   + k * b(p-1,i-1)  i=2,3,...n
      %%  b(p,i) = b(p-1,i-1) + conj(k) * f(p-1,i)    i=2,3,...n
      %%  remember f(p,1) is not stored, so don't calculate it; make f(p,2)
      %%  the first element in forw_err.  b(p,n) isn't calculated either.
      nn = length(forw_err);
      forw_new = forw_err(2:nn)   + refl_coeff .* back_err(2:nn);
      back_err = back_err(1:nn-1) + conj(refl_coeff) .* forw_err(1:nn-1);
      forw_err = forw_new;
      end
    end
  %% end of for loop
  %%
  varargout{1} = [1 coeffs];
  varargout{2} = residual;
  if ( nargout>=3 )
    varargout{3} = lattice;
    end
  end
end

%% Copyright (C) 2006 Peter V. Lanspeary
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
%% as published by the Free Software Foundation; either version 2,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software Foundation,
%% Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% usage:
%%    [psd,freq_out] = burg_psd(data,poles,freq,sample_f,
%%                              range,method,plot_type,stop_crit)
%%
%% Calculate Burg maximum-entropy power spectrum.  This is a wrapper for
%% functions "burg_model" and "ar_psd" which perform the argument checking.
%% See "help burg_model" and "help ar_psd" for further details.
%%
%% ARGUMENTS:
%%     All but the first two arguments are optional and may be empty.
%%   data      %% [vector] sampled data
%%
%%   poles     %% [integer scalar] required number of poles of the AR model
%%
%%   freq      %% [real vector] frequencies at which power spectral density
%%             %%               is calculated
%%             %% [integer scalar] number of uniformly distributed frequency
%%             %%          values at which spectral density is calculated.
%%             %%          [default=256]
%%
%%   sample_f  %% [real scalar] sampling frequency (Hertz) [default=1]
%%
%%
%% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
%%   Control-string arguments can be in any order after the other arguments.
%%
%%
%%   range     %% 'half',  'onesided' : frequency range of the spectrum is
%%             %%       zero up to but not including sample_f/2.  Power from
%%             %%       negative frequencies is added to the positive side of
%%             %%       the spectrum.
%%             %% 'whole', 'twosided' : frequency range of the spectrum is
%%             %%       -sample_f/2 to sample_f/2, with negative frequencies
%%             %%       stored in "wrap around" order after the positive
%%             %%       frequencies; e.g. frequencies for a 10-point 'twosided'
%%             %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
%%             %% 'shift' : same as 'whole' but with the first half of the
%%             %%       spectrum swapped with second half to put the zero-
%%             %%       frequency value in the middle. (See "help fftshift".)
%%             %%       If "freq" is vector, 'shift' is ignored.
%%             %% If model coefficients "ar_coeffs" and "residual" are real,
%%             %% default range is 'half', otherwise default range is 'whole'.
%%
%%   method    %% 'FFT':  use FFT to calculate power spectrum.
%%             %% 'poly': calculate power spectrum as a polynomial of 1/z
%%             %% N.B. this argument is ignored unless "freq" is scalar.
%%             %%      The default is 'poly' unless the number of frequency
%%             %%      values -1 (freq) is an integer power of 2.
%%  
%%   plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
%%             %%       specifies the type of plot.  The default is 'plot',
%%             %%       which means linear-linear axes. 'squared' is the same
%%             %%       as 'plot'. 'db' plots "10*log10(psd)".
%%
%%   stop_crit %% 'FPE' -- apply "final prediction error" criterion
%%             %% 'AIC' -- apply "Akaike information criterion"
%%             %%    (Kay & Marple, 1981) to limit the number of poles so that
%%             %%    spurious poles are not added when the whitened data has
%%             %%    no more information in it. The default is to NOT apply
%%             %%    the AIC or FPE.
%%
%% RETURNED VALUES:
%%     If return values are not required by the caller, the spectrum
%%     is plotted and nothing is returned.
%%   psd       %% [real vector] power-spectrum estimate
%%   freq      %% [real vector] frequency values


function [psd,freq]=burg_psd(data,poles,varargin)
%%
nvarargin=length(varargin);
stop_crit=[];
%%
%% Search for a "stop_crit" arg. If found, remove it
%% from "varargin" list and feed it to burg_model instead.
for iarg = 1: nvarargin
  arrgh = varargin{iarg};
  if ( ischar(arrgh) && (strcmp(arrgh,'FPE') || strcmp(arrgh,'AIC')) )
    stop_crit=arrgh;
    if ( nvarargin>1 )
      varargin{iarg}= [];
    else
      varargin={};
      end
    end
  end
%%
[ar_coeffs,residual]=burg_model(data,poles,stop_crit);
if ( nargout==0 )
  ar_psd(ar_coeffs,residual,varargin{:});
elseif ( nargout==1 )
  [psd,freq]=ar_psd(ar_coeffs,residual,varargin{:});
elseif ( nargout==2 )
  [psd,freq]=ar_psd(ar_coeffs,residual,varargin{:});
  end
end

%% Copyright (C) 2006 Peter V. Lanspeary
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
%% as published by the Free Software Foundation; either version 2,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software Foundation,
%% Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% Usage:
%%   [spectra,freq] = welch_psd( x, window, sample_f, overlap, padding,...
%%                               range, plot_type, detrend, sloppy)
%%        Estimate power spectrum of time-series data "x" by the Welch (1967)
%%        (periodogram/FFT) method.  All arguments except "x" are optional.
%%
%%   [spectra,freq] = welch_psd( x, y, window, sample_f, overlap, padding,...
%%                               range, plot_type, detrend, sloppy, results)
%%        Estimate cross-spectral density, transfer function and/or coherence
%%        functions  of time-series input data "x" and output data "y" by the
%%        Welch (1967) (periodogram/FFT) method.  Requires arguments "x", "y",
%%        "results" = 'cross', 'trans' or 'coher'.  All other arguments are
%%        optional.  All spectra are returned in matrix "spectra".
%%
%% The data is divided into segments.  If "window" is a vector,
%% each segment has the same length as "window" and is multiplied by "window"
%% before (optional) zero-padding and calculation of its periodogram.  If
%% "window" is a scalar, each segment has a length of "window" and a Hann
%% window is used.  Hann (hanning), hamming, bartlett, blackman, flattopwin
%% etc are available as separate Matlab sigproc or Octave functions.
%%
%% The power spectrum is the mean of the periodograms, scaled so that area
%% under the power spectrum is the same as the mean square of the data.
%% NOTE that this equivalence is supposed to be exact, but in practice I
%%      have found a mismatch of up to 0.5% when comparing area under a
%%      periodogram with the mean square of the data.
%%
%% ARGUMENTS
%% All but the first argument are optional and may be empty, except that
%% the "results" argument may require the second argument to be "y".
%%
%% x           %% [non-empty vector] system-input time-series data
%% y           %% [non-empty vector] system-output time-series data
%%
%% window      %% [real vector] of window-function values between 0 and 1; the
%%             %%       data segment has the same length as the window.
%%             %% [integer scalar] length of each data segment (and
%%             %%       unpadded FFT); default value is window=sqrt(length(x))
%%             %%       rounded down to the nearest integer power of 2; see
%%             %%    'sloppy' argument. Default window shape is Hann (hanning).
%%
%% sample_f    %% [real scalar] sampling frequency (Hertz); default=1.0
%%
%% overlap     %% [real scalar] segment overlap factor  0 <= overlap < 1,
%%             %%       The default is 0.5.
%%
%% padding     %% [integer scalar] number of samples of zero padding (per FFT)
%%             %%       default is zero
%%
%% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
%%   Control-string arguments must be after the other arguments but can be in
%%   any order.
%%  
%%   range     %% 'half',  'onesided' : frequency range of the spectrum is
%%             %%       zero up to but not including sample_f/2.  Power from
%%             %%       negative frequencies is added to the positive side of
%%             %%       the spectrum, but not at zero or Nyquist (sample_f/2)
%%             %%       frequencies.  This keeps power equal in time and
%%             %%       spectral domains.  See reference [2].
%%             %% 'whole', 'twosided' : frequency range of the spectrum is
%%             %%       -sample_f/2 to sample_f/2, with negative frequencies
%%             %%       stored in "wrap around" order after the positive
%%             %%       frequencies; e.g. frequencies for a 10-point 'twosided'
%%             %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
%%             %% 'shift' : same as 'whole' but with the first half of the
%%             %%       spectrum swapped with second half to put the zero-
%%             %%       frequency value in the middle. (See "help fftshift".)
%%             %%  Default is 'half' for real data, 'whole' for complex data.
%%
%% plot_type   %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
%%             %%       specifies the type of plot.  The default is 'plot',
%%             %%       which means linear-linear axes. 'squared' is the same
%%             %%       as 'plot'. 'db' plots "10*log10(psd)".
%%
%% detrend     %% 'no-detrend' -- do NOT remove a mean value from the data
%%             %% 'short' -- remove the mean value of each segment from each
%%             %%            segment of the data
%%             %% 'long'  -- remove the mean value from the data (before
%%             %%            splitting into segments). This is the default.
%%
%% sloppy      %% 'sloppy': FFT length is rounded up to the nearest integer
%%             %%       power of 2, segment length is rounded up unless the
%%             %%       "window" argument is specified as a vector. FFT length
%%             %%       is adjusted after addition of any padding.  The
%%             %%       default is to use exactly the segment lengths and
%%             %%       padding lengths (hence FFT length) specified in
%%             %%       argument list.
%%
%% results     %% specifies what results to return (in the order specified
%%             %%   and as many as desired).
%%             %% 'power' calculate power spectral density of "x"
%%             %% 'cross' calculate cross spectral density of "x" and "y"
%%             %% 'trans' calculate transfer function of a system with
%%             %%         input "x" and output "y"
%%             %% 'coher' calculate coherence function of "x" and "y"
%%             %% 'ypower' calculate power spectral density of "y"
%%             %%  The default is 'power'.
%%
%% RETURNED VALUES:
%%   If return values are not required by the caller, the results are
%%     plotted and nothing is returned.
%%
%% spectra     %% [real matrix] rows of the matrix contain results in the
%%             %%               same order as specified by "results" arguments.
%%             %%               Each row contains one of the result vectors.
%%
%% freq        %% [real vector] frequency values
%%
%% REFERENCES
%%  [1] Peter D. Welch (June 1967):
%%   "The use of fast Fourier transform for the estimation of power spectra:
%%   a method based on time averaging over short, modified periodograms."
%%   IEEE Transactions on Audio Electroacoustics, Vol AU-15(6), pp 70-73
%%
%%  [2] William H. Press and Saul A. Teukolsky and William T. Vetterling and
%%               Brian P. Flannery",
%%   "Numerical recipes in C, The art of scientific computing", 2nd edition,
%%      Cambridge University Press, 2002 --- Section 13.7.
%%  [3] Paul Kienzle (1999-2001): "pwelch", http://octave.sourceforge.net/


function [varargout] = welch_psd(x,varargin)
%%
%% Check fixed arguments
if ( nargin <= 0 )
  error( 'welch_psd: error: Need at least 1 arg. Use "help welch_psd"' );
elseif ( isempty(x) || ~isvector(x) )
  error( 'welch_psd: error: arg 1 (x) must be vector.' );
else
%%  force x to be ROW vector
  if ( size(x,2)==1 )
    x=reshape(x,1,[]);
    end
%%
%% get the second data vector if it is required
  need_y = 0;
  x_len = length(x);
  nvarargin = length(varargin);
  for iarg=1:nvarargin
    arg = varargin{iarg};
    if ( ~isempty(arg) && ischar(arg) && ( strcmp(arg,'cross') || ...
             strcmp(arg,'trans') || strcmp(arg,'coher') ))
      need_y = 1;
      if ( nargin<2 || isempty(varargin{1}) || ~isvector(varargin{1}) || ...
           length(varargin{1}) ~= x_len )
        error( 'welch_psd: arg error: y must be vector & same length as x.' );
        return;
        end
      y = varargin{1};
      %% force ROW vector
      if ( size(x,2)==1 )
        y = reshape(y,1,[]);
        end
      break;
      end
    end
%%
%%  default values for optional areguments
  is_win    = 0;
  sample_f  = 1;
  overlap   = 0.5;
  padding   = 0;
  range     = ~isreal(x) || ( need_y && ~isreal(y) );
  is_sloppy = 1;
  plot_type = 1;
  rm_mean   = 2; %% remove mean
  n_results = 0;
  do_power  = 0;
  do_cross  = 0;
  do_trans  = 0;
  do_coher  = 0;
  do_ypower = 0;
%%
%%  decode and check optional arguments
  end_numeric_args = 0;
  for iarg = 1+need_y:nvarargin
    arg = varargin{iarg};
    %% skip empty arguments
    end_numeric_args = end_numeric_args || ischar(arg);
    if ( isempty(arg) )
      1;
    %% first 4 optional arguments are numeric -- in fixed order
    elseif ( ~ischar(arg) )
      if ( end_numeric_args )
        error( 'welch_psd: error: control arg must be string' );
        return;
      %%
      elseif ( iarg-need_y == 1 )
        window = arg;
        is_sloppy = 0;
        if ( isscalar(window) )
          is_win = 1;
        elseif ( isvector(window) )
          is_win = length(window);
          if ( size(window,1)>1 )
            window = reshape(window,1,[]);
            end
        else
          is_win = 0;
          end
        if ( ~is_win )
        error('welch_psd: arg error: window must be vector or segment length');
          return;
        elseif ( is_win==1 && ( ~isreal(window) || ...
                 fix(window)~=window || window<3 || x_len<window ) )
       error('welch_psd: error: arg window must be integer, >3 & <=length(x)');
          return;
        elseif ( is_win>1 && ( ~isreal(window) || any(window<0) ) )
          error( 'welch_psd: arg error: window vector must be real and >=0');
          return;
          end
      %%
      elseif ( iarg-need_y == 2 )
        sample_f = arg;
        if ( ~isscalar(sample_f) || ~isreal(sample_f) || sample_f<0 )
          error( 'welch_psd: arg error: sample_f must be real scalar >0.' );
          return;
          end
      %%
      elseif ( iarg-need_y == 3 )
        overlap = arg;
        if (~isscalar(overlap) || ~isreal(overlap) || overlap<0 || overlap>=1)
          error('welch_psd: arg error: overlap must be between 0 and 1.0');
          return;
          end
      elseif ( iarg-need_y == 4 )
      %%
        padding = arg;
        if ( ~isscalar(padding) || ...
             ~isreal(padding) || fix(padding)~=padding || padding<0 )
          error( 'welch_psd: arg error: padding must be integer >=0' );
          return;
          end
      else
      %%
        error( 'welch_psd: error: control arg must be string' );
        return;
        end
  %%
  %% decode control-string arguments
    elseif ( strcmp(arg,'sloppy') )
      is_sloppy = ~is_win || is_win==1;
    elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
      plot_type = 1;
    elseif ( strcmp(arg,'semilogx') )
      plot_type = 2;
    elseif ( strcmp(arg,'semilogy') )
      plot_type = 3;
    elseif ( strcmp(arg,'loglog') )
      plot_type = 4;
    elseif ( strcmp(arg,'db') )
      plot_type = 5;
    elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
      range = 0;
    elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
      range = 1;
    elseif ( strcmp(arg,'shift') )
      range = 2;
    elseif ( strcmp(arg,'long') )
      rm_mean = 2;
    elseif ( strcmp(arg,'short') )
      rm_mean = 1;
    elseif ( strcmp(arg,'no-detrend') )
      rm_mean = 0;
    elseif ( strcmp(arg, 'power' ) )
      if ( ~do_power )
        n_results = n_results+1;
        do_power = n_results;
        end
    elseif ( strcmp(arg, 'cross' ) )
      if ( ~do_cross )
        n_results = n_results+1;
        do_cross = n_results;
        end
    elseif ( strcmp(arg, 'trans' ) )
      if ( ~do_trans )
        n_results = n_results+1;
        do_trans = n_results;
        end

[David Bateman]

Peter,

Wow, the shear quantity of the comparison data is awesome... Frankly,
I'm not sure I have the time to fully understand, but its seems to me
from a brief reading of your message that

1) You don't think full matlab compatibility is desireable
2) Your code is more compatible that the existing code
3) Your have additional optional arguments

The additional arguments are not a problem as they don't affect
compatibility. Paul stated that its a pity to break compatibility with
existing octave code, but I thinks its more important to adapt to the
matlab interface to attract new users. As for the first point, its seems
to me your the expert.

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..

D.

--
David Bateman                                David.Bateman@...
Motorola Labs - Paris                        [phone] (Ph)
Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[pkienzle-2]

On Sep 26, 2006, at 5:55 AM, David Bateman wrote:

> Peter,
>
> Wow, the shear quantity of the comparison data is
> awesome... Frankly, I'm not sure I have the time
> to fully understand, but its seems to me
> from a brief reading of your message that
>
> 1) You don't think full matlab compatibility is desireable
> 2) Your code is more compatible that the existing code
> 3) Your have additional optional arguments
>
> The additional arguments are not a problem as they
> don't affect compatibility. Paul stated that its a
> pity to break compatibility with existing octave code,
> but I thinks its more important to adapt to the
> matlab interface to attract new users. As for the
> first point, its seems to me your the expert.
>
> 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..

I agree with David that we should provide the new
versions so that new code is compatible.  However,
it is a pity to break compatibility with existing
_Matlab_ code, even though it is from an older version.

Do we want to save code compatible with R11?

- Paul

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[Peter Lanspeary]

Hi David,

On Tue, Sep 26, 2006 at 11:55:11AM +0200, David Bateman wrote:
> Wow, the shear quantity of the comparison data ...

I am sorry to report that my comparison contains a number of mistakes.
I attach (first attachment) a new, more accurate version of the comparison.

> ... its seems to me from a brief reading of your message that
> 1) You don't think full matlab compatibility is desireable
> 2) Your code is more compatible that the existing code
> 3) Your have additional optional arguments
> The additional arguments are not a problem as they don't affect
> compatibility. Paul stated that its a pity to break compatibility with
> existing octave code, but I thinks its more important to adapt to the
> matlab interface to attract new users...

I agree that it would be nice to be compatible with matlab, but I
would prefer to avoid duplicating matlab bugs and badly-designed parts of the
interface (see the second attachment).  For simple methods with few parameters,
such as the Burg method, the best choice of interface is obvious and there is
no problem.  The Welch method is much more complicated and tricky, and a good
pwelch interface is more difficult to design.  Even Mathworks have not set the
interface in concrete; their documentation (e.g.
http://www.mathworks.com/access/helpdesk/help/toolbox/signal/index.html?/
access/helpdesk/help/toolbox/signal/f12-6587.html)
seems to steer users to the newer "psd(spectrum.welch,...)", which has
differently structured arguments and different defaults.  A colleague says
that it is much easier to use than pwelch. Pwelch might not seem so attractive.

> 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..

I can do the following without too much trouble.

1) The Burg-function arguments (both my functions and octave's) are
   compatible with matlab, so replacing the octave arburg and octave
   pburg gives
   a) an interface compatible with matlab
   b) correct results with the following differences in behaviour
      a) returned frequencies are always Hz, never rad/sec
      b) does not duplicate the matlab R14 periodogram bug
      c) for onesided spectra, returns nfft or length(nfft) PSD values
      d) additional arguments give extra features.
   However, fixing the octave arburg (see the third attachment) delivers all
   of the above except 2)d).

2) The arguments of welch_psd would be made more compatible with matlab pwelch,
   but not __completely__ compatible.  If we replace "pad" with "nfft", where
   nfft=length(window)+pad, and make "Fs" the 5th argument,
      [Pxx,f]=welch_psd(x,window,Fs,overlap,pad,range,units,trend,sloppy)
   becomes
      [Pxx,f]=welch_psd(x,window,overlap,nfft,Fs,range,units,trend,sloppy)
   Matlab pwelch is
         [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,range)

   I would then change the default window function to Hamming and change the
   overlap from a fraction (e.g. 0.5) to a percentage (e.g. 50) so that it
   becomes compatible with matlab's "psd(spectrum.welch,...)".
   The remaining 3 differences between corresponding input arguments are:
   +----------+--------------------+--------------------+--------------------+
   | Argument |   matlab pwelch    |   welch_psd        |Is compatibility with
   +----------+--------------------+--------------------+ pwelch possible ??
   +----------+--------------------+--------------------+--------------------+
   | window   | default_length=    | default_length=    | yes, but would     |
   |          |    length(x)/8     |    sqrt(length(x)) | make bad spectra   |
   |          |Is incompatible with| rounded up to 2^n. |                    |
   |          | psd(spectrum.welch)|                    |                    |
   +----------+--------------------+--------------------+--------------------+
   | nfft     | length of FFT or   | =length(window)+pad| is compatible      |
   |          +--------------------+--------------------+--------------------+
   |          | list of frequencies|                    | NO - not available |
   +----------+--------------------+--------------------+--------------------+
   | overlap  | overlap samples    | overlap percentage |would be undesirable|
   +----------+--------------------+--------------------+--------------------+
   Also, returned frequencies would always be in Hz, never rad/sec

I would be happy build new pwelch, arburg, ar_psd and pburg
  a) if that is what you propose,
  b) if the above details are acceptable to you, and
  c) if I have access to the CVS.
I have a sourceforge account, but I am not registered as an  octave-forge
developer. I guess I'll need to "read the instructions" first...

Thanks,
Peter Lanspeary


COMPARISON OF WELCH and BURG POWER-SPECTRUM FUNCTIONS
P. V. Lanspeary

Contents
--------
Section 1. Comparison of arguments. Summary.
Section 2. Comparison of behaviours. Summary.
Section 3. pburg: Full compatibility with Matlab ?
Section 4. Comparison of arguments. Details.
Section 5. Comparison of behaviours. Details.


Section 1. Comparison of arguments. Summary.
--------------------------------------------
I have compared the arguments of welch_psd.m, burg_model.m and burg_psd.m
with equivalent functions in octave-forge and Matlab-sigproc.
This is a summary. See section (4) for details.

(a) welch_psd with pwelch (octave) and pwelch (matlab):
  The behaviour of welch_psd and pwelch is complicated, so details of the
  comparison are messy and complicated. In summary, the news is not
  good.  Overall functionality is similar (not the same). Only the
  output args and first input arg are compatible.  

(b) burg_model with arburg (octave) and arburg (matlab):
  Arguments are compatible except that burg_model has an extra optional
  input argument.  If each function is called with the same 2 arguments,
  behaviour should (according to the documentation) be the same.

(c) burg_psd with pburg (octave) and pburg (matlab):
  Arguments are compatible except that Octave's pburg has one an extra
  optional input argument and burg_psd has 3 extra optional input arguments.


Section 2. Comparison of behaviours. Summary.
---------------------------------------------
I have run _two_ test cases on the "Welch" functions and two more on
the "Burg" functions.  One test case has a real signal (data); the
other has a complex signal.  This is a summary.  See section (5) for
details.

(a) pwelch (octave), pwelch (matlab), welch_psd:
  Arguments are incompatible; so to do the same thing with each function
  requires different arguments in each case.  If each function is
  called with the same arguments, behaviour (except in special cases) is
  different.  Producing the same behaviour requires care. Where documentation
  requirees the same behaviour, output is the same to an accuracy of no worse
  than 10x machine precision (eps*10).

(b) arburg (octave), arburg (matlab), burg_model:
  For the real-data test case, results are exactly the same --- except that
  reflection coefficients (k) from octave arburg have the wrong sign.
  For complex data, the Matlab arburg gives the same (correct) answer as
  burg_model, but the octave arburg does not give the correct answer.

(c) pburg (octave), pburg (matlab), burg_psd:
  For 'twosided' spectra Matlab pburg gives the same answers as burg_psd, but
  for 'onesided' spectra, the first and last values in the spectrum are too
  small by a factor of 2.
  The octave pburg does not work for complex data because it calls the octave
  arburg.  The octave __power.m works correctly for both real and complex
  data.


Section 3. pburg: Full compatibility with Matlab ?
--------------------------------------------------

>> What would it take to create a fully Matlab compatiable pburg function?
Arguments are fully compatible, but with optional extra arguments.
Compatible behaviour is harder to achieve, and I am not sure that this is
desirable. For example, the Matlab pburg
a) assumes units of frequency are radians/sec unless sampling
   frequency, Fs, is specified.  If Fs is specified, the units are Hertz.
   Octave pburg uses only Hertz.
b) returns only nfft/2+1 (frequency,spectral-density) coordinates if the
   spectrum is onesided.  Octave pburg returns nfft values.
c) plots the spectral density in dB rather than correct physical units;
   the graph should be a power spectrum of a physical quantity rather
   than the response of a filter.
d) has a error in onesided spectra (the first and last values of spectral
   density are too small by a factor of 2).
A completely fixed Octave pburg would seem to have more appropriate behaviour.


Section 4. Comparison of arguments. Details.
--------------------------------------------

(a) pwelch (octave), pwelch (matlab), welch_psd:
  The number of arguments is large and behaviour of the functions is a bit
  complicated ... for the sake of brevity this comparison/description
  therefore contains some minor approximations. First we write the first line
  of each function, adjusting argument names to make them more consistent::

  1) octave [Pxx,Pci,f]=pwelch(x,nfft,Fs,window,overlap,ci,range,units,trend)
  2) matlab [Pxx,f]    =pwelch(x,window,overlap,nfft,Fs,range)
  3)        [Pxx,f]=welch_psd(x,window,Fs,overlap,pad,range,units,trend,sloppy)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave pwelch     |   matlab pwelch    |   welch_psd        |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  |  Pxx     | spectral density   | spectral density   | spectral density   |
  |          | units specified by | ["x" units]/Hertz  | ["x" units]/Hertz  |
  |          | "units" arg        | or per (rad/s)     |                    |
  +----------+--------------------+--------------------+--------------------+
  | Pci      | confidence interval| not available      | not available      |
  +----------+--------------------+--------------------+--------------------+
  | f        | frequency (Hertz)  | frequency          | frequency (Hertz)  |
  |          |                    | (Hertz if given Fs)|                    |
  |          |                    | (rad/s without Fs) |                    |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | complex data                                                 |
  +----------+--------------------+--------------------+--------------------+
  | window   | windowing vector, or the length of data segment              |
  |          +--------------------+--------------------+--------------------+
  |          | default = Hann     | default = Hamming  | default = Hann     |
  |          | default_length=nfft| default_length=    | default_length=    |
  |          |                    | length(x)/8        |    sqrt(length(x)) |
  |          |                    |                    |   rounded up to 2^n|
  +----------+--------------------+--------------------+--------------------+
  | Fs       | sampling frequency | sampling frequency | sampling frequency |
  |          | default=2 Hz       | default= 1 Hz      | default=1 Hz       |
  +----------+--------------------+--------------------+--------------------+
  | nfft     | length of FFT      | length of FFT or   | not available; uses|
  |          | default = 256      | list of frequencies| length(window)+pad |
  |          |                    | default=256 or 2^n |    as nfft         |
  |          |                    |   >length(window)  |                    |
  +----------+--------------------+--------------------+--------------------+
  | overlap  | number of overlap  | number of overlap  | overlap fraction   |
  |          | samples, default=  | samples, default=  | default=0.5, i.e.  |
  |          | length(window)/2   | length(window)/2   | 50%                |
  +----------+--------------------+--------------------+--------------------+
  | pad      | not available      | not available      | length of zero     |
  |          |                    |                    | padding on FFT     |
  |          |=nfft-length(window)|=nfft-length(window)| default=0          |
  +----------+--------------------+--------------------+--------------------+
  | ci       | confidence interval| not available      | not available      |
  +----------+--------------------+--------------------+--------------------+
  | range    | 'half'             | 'half', 'onesided' | 'half', 'onesided' |
  |          | 'whole'            | 'whole', 'twosided'| 'whole', 'twosided'|
  |          |                    |                    | 'shift'            |
  +----------+--------------------+--------------------+--------------------+
  | units    | plot scaling       | not available      | 'plot', 'semilogx' |
  |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
  |          |                    |                    | 'loglog', 'squared'|
  +----------+--------------------+--------------------+--------------------+
  | trend    | remove trend       | not available      |remove trend, 'long'|
  |          | 'mean', 'linear'   | -- does not        |'short','no-detrend'|
  |          | default=none       |    remove mean     | default='long'     |
  +----------+--------------------+--------------------+--------------------+
  | sloppy   | not available      | not available      | round nfft to 2^n  |
  +----------+--------------------+--------------------+--------------------+
    The above table does not describe the special cases where the octave pwelch
    and welch_psd are used for calculating cross-power spectra, coherence
    and/or transfer functions.

(b) arburg (octave), arburg (matlab), burg_model:
  (1) octave  [a,v,k]=arburg(x,p)
  (2) matlab  [a,v,k]=arburg(x,p)
  (3)         [a,v,k]=burg_model(x,p,stop)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave arburg     |   matlab arburg    |   burg_model       |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | a        | autoregression coefficients                                  |
  +----------+--------------------+--------------------+--------------------+
  | v        | mean square residual (noise)                                 |
  +----------+--------------------+--------------------+--------------------+
  | k        | reflection coefficients                                      |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | real data          | complex data       | complex data       |
  +----------+--------------------+--------------------+--------------------+
  | p        | number of poles (required arg)                               |
  +----------+--------------------+--------------------+--------------------+
  | stop     | not available      | not available      | stopping criterion |
  |          |                    |                    | 'FPE','AIC'        |
  |          |                    |                    | default='none'     |
  +----------+--------------------+--------------------+--------------------+

(c) pburg (octave), pburg (matlab), burg_psd:
  (1) octave  [Pxx,f]=pburg(x,p,nfft,Fs,range,units)
  (2) matlab  [Pxx,f]=pburg(x,p,nfft,Fs,range)
  (3)         [Pxx,f]=burg_psd(x,p,nfft,Fs,range,units,method,stop)

  +----------+--------------------+--------------------+--------------------+
  | Argument |  octave pburg      |   matlab pburg     |   burg_psd         |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  |  Pxx     | spectral density   | spectral density   | spectral density   |
  |          | ["x" units]/Hertz  | ["x" units]/Hertz  | ["x" units]/Hertz  |
  |          |                    | or per (rad/s)     |                    |
  +----------+--------------------+--------------------+--------------------+
  | f        | frequency (Hertz)  | frequency          | frequency (Hertz)  |
  |          |                    | (Hertz if given Fs)|                    |
  |          |                    | (rad/s without Fs) |                    |
  +----------+--------------------+--------------------+--------------------+
  +----------+--------------------+--------------------+--------------------+
  | x        | real data          | complex data       | complex data       |
  +----------+--------------------+--------------------+--------------------+
  | p        | number of poles (required arg)                               |
  +----------+--------------------+--------------------+--------------------+
  | nfft     | number of frequency| length of FFT or   | number of frequency|
  |          | values             | list of frequencies| values or list of  |
  |          |                    |                    | frequencies        |
  |          | default = 256      | default=256        | default=256        |
  +----------+--------------------+--------------------+--------------------+
  | Fs       | sampling frequency | sampling frequency | sampling frequency |
  |          | default=2 Hz       | default= 1 Hz      | default=1 Hz       |
  +----------+--------------------+--------------------+--------------------+
  | range    | 'half'             | 'half', 'onesided' | 'half', 'onesided' |
  |          | 'whole'            | 'whole', 'twosided'| 'whole', 'twosided'|
  |          |                    |                    | 'shift'            |
  +----------+--------------------+--------------------+--------------------+
  | units    | plot scaling       | not available      | 'plot', 'semilogx' |
  |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
  |          |                    |                    | 'loglog', 'squared'|
  +----------+--------------------+--------------------+--------------------+
  | method   | not available      | not available      | 'FFT', 'poly'      |
  +----------+--------------------+--------------------+--------------------+
  | stop     | not available      | not available      | stopping criterion |
  |          |                    |                    | 'FPE','AIC'        |
  |          |                    |                    | default='none'     |
  +----------+--------------------+--------------------+--------------------+


Section 5. Comparison of behaviours. Details.
---------------------------------------------
The test platforms are
   a) Octave 2.1.69 + octave-forge-2006.07.09 on Debian 3.1r1
      with sourceforge.net/tracker bugfix #1530283   for arburg.m
      with sourceforge.net/tracker bugfix #1530309   for __power.m
   b) Matlab 7.0.1.24704(R14)SP1 + sigproc toolbox v6.2.1, on MS-windows XP
      (this is the most recent Matlab version available to me)

(a) pwelch (octave), pwelch (matlab), welch_psd:
  The test script is

     format long
     rand('seed',2038014164);
     a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
     white = rand(1,16384);
     signal = detrend(filter(0.70181,a,white));
     skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
     %% REAL-DATA TEST
    [psd1,f1]=welch_psd(signal,hann(256),25,0.5,0,'whole','plot','no-detrend');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd2,f2]=pwelch(signal,256,25,hann(256),128,'whole','squared','none');
     else
       disp( 'MATLAB' )
       [psd2,f2]=pwelch(signal,hann(256),256*0.5,256,25,'twosided');
       end
     sum_psd1=sum(psd1)
     sum_psd2=sum(psd2)
     maxdiffr=max(abs(psd1./reshape(psd2,1,[])-1))  %% max relative difference
     %% COMPLEX-DATA TEST
    [psd3,f3]=welch_psd(skewed,hann(256),25,0.5,0,'whole','plot','no-detrend');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd4,f4]=pwelch(skewed,256,25,hann(256),128,'whole','squared','none');
     else
       disp( 'MATLAB' )
       [psd4,f4]=pwelch(skewed,hann(256),256*0.5,256,25,'twosided');
       end
     sum_psd3=sum(psd3)
     sum_psd4=sum(psd4)
     maxdiffc=max(abs(psd3./psd4'-1))  %% maximum relative difference
     eps

  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
  | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | sum_psd1 | 3.62735601946465      | 3.61040879075053          |
  | sum_psd2 | 3.62735601946465      | 3.61040879075053          |
  | maxdiffr | 1.33226762955019e-15  | 1.110223024625157e-015    |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | sum_psd3 | 3.62735601946465      | 3.61040879075054          |
  | sum_psd4 | 3.62735601946465      | 3.61040879075054          |
  | maxdiffc | 1.33226762955019e-15  | 1.332267629550188e-015    |
  +----------+-----------------------+---------------------------+
    We should have sum_psd1==sum_psd2, maxdiffr==0, sum_psd3==sum_psd4,
    maxdiffc==0.


(b) arburg (octave), arburg (matlab), burg_model:
  The additional test script (carrying on from previous test) is

     %% REAL-DATA TEST
     [ar1,vr1,kr1]=burg_model(signal,5);
     [ar2,vr2,kr2]=arburg(signal,5);
     maxdifar=max(abs(ar2./ar1-1))  %% maximum relative difference
     vr1/vr2-1
     maxdifkr=max(abs(reshape(kr2,1,[])./kr1-1)) %% kr2 has inconsistent shape
     %% COMPLEX-DATA TEST
     [ac1,vc1,kc1]=burg_model(skewed,5);
     [ac2,vc2,kc2]=arburg(skewed,5);
     maxdifac=max(abs(ac2./ac1-1))  %% maximum relative difference
     vc2/vc1-1
     maxdifkc=max(abs(reshape(kc2,1,[])./kc1-1))

  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | burg_model vs. arburg | burg_model vs. arburg     |
  | eps      | 2.22044604925031e-16  |  2.220446049250313e-016   |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | maxdifar | 8.88178419700125e-16  |  0                        |
  | vr1/vr2-1| -5.55111512312578e-16 |  0                        |
  | maxdifkr | 2.00000000000000      |  0                        |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | maxdifac | 2.15561046016939      |  0                        |
  | vc2/vc1-1| 7.679016 - 0.0026252i |  0                        |
  | maxdifkc | 1.68536747392439      |  5.590172498962309e-017   |
  +----------+-----------------------+---------------------------+
     We should have maxdifar==0, vr1/vr2-1==0, maxdifkr==0, maxdifac==0,
     vc2/vc1-1==0, maxdifkc==0.
     ... Octave arburg gives the correct autoregression coefficients and
     variance (ar1,vr1) for real data but the reflection coefficients have
     the wrong sign; for complex data, the octave arburg gives the wrong
     answers.  Matlab pburg and burg_model give exactly the same answer
     for both real and complex data.
     Now compare __power with ar_psd using coefficients (ac1,vc1)
     calculated by burg_model  ...

     [psd3,f3]=ar_psd(ac1,vc1,256,25,'whole','squared');
     [psd4,f4]=__power(sqrt(vc1),ac1,256,25,'whole','squared');
     sum_psd3=sum(psd3)
     %%     sum_psd3 = 3.63449149426036e+00 + 4.89110705062523e-19i
     sum_psd4=sum(psd4)
     %%     sum_psd4 = 3.63449149426036
     maxdiffc=max(abs(psd3./psd4-1))
     %%     maxdiffc =  3.33066934574812e-16

     ... we conclude that __power works OK for complex filter coefficients

(b) pburg (octave), pburg (matlab), burg_psd:
  The additional test script (carrying on from previous test) is

     %% REAL-DATA TEST
     [psd5,f5]=burg_psd(signal,5,256,25,'whole');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd6,f6]=pburg(signal,5,256,25,'whole');
     else
       disp( 'MATLAB' )
       [psd6,f6]=pburg(signal,5,256,25,'twosided');
       end
     sum_psd5=sum(psd5)
     sum_psd6=sum(psd6)
     maxdiffr=max(abs(psd5./reshape(psd6,1,[])-1))  %% max relative difference
     %% COMPLEX-DATA TEST
     [psd7,f7]=burg_psd(skewed,5,256,25,'whole');
     if ( strcmp(version,'2.1.69') )
       disp( 'OCTAVE 2.1.69' )
       [psd8,f8]=pburg(skewed,5,256,25,'whole');
     else
       disp( 'MATLAB' )
       [psd8,f8]=pburg(skewed,5,256,25,'twosided');
       end
     sum_psd7=sum(psd7)
     sum_psd8=sum(psd8)
     maxdiffc=max(abs(psd7./reshape(psd8,1,[])-1))  %% max relative difference
 
  +----------+-----------------------+---------------------------+
  |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
  |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
  +----------+-----------------------+---------------------------+
  |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
  | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
  +----------+-----------------------+---------------------------+
  | real data                                                    |
  | sum_psd5 | 3.63449149426037      | 3.62433951286238          |
  | sum_psd6 | 3.63449149426037      | 3.62433951286238          |
  | maxdiffr | 1.24344978758018e-14  | 5.551115123125783e-016    |
  +----------+-----------------------+---------------------------+
  | complex data                                                 |
  | sum_psd7 | 3.634491494 +4.9e-19i | 3.62433951286238          |
  | sum_psd8 | 3.64552153406807      | 3.62433951286238          |
  | maxdiffc | 6.15692319331572      | 6.661338147750939e-016    |
  +----------+-----------------------+---------------------------+
     We should have sum_psd5==sum_psd6, maxdiffr==0, sum_psd7==sum_psd8,
     maxdiffc==0.
     Octave pburg doesn't work for complex data because it calls  arburg.
     Twosided spectra from matlab pburg are OK, but for "onesided" spectrum
      [psd8,f8]=pburg(skewed,5,256,25,'onesided');
    the first and last values in the spectrum from matlab are too small by
    a factor of 2.

Bugs and Design Problems with Burg and Welch functions
------------------------------------------------------
Peter Lanspeary, 28 Sept 2006

1) The Burg functions.

If we ignore additional arguments, both octave (arburg, pburg) and my efforts
(burg_model, purg_psd) have interfaces which are compatible with matlab.
However, behaviour is different.
 (a) The octave arburg does not work with complex data, and with real data
     the returned reflection coefficients have the wrong sign.  I have
     modified the octave arburg code to fix these problems so that its
     behaviour is identical to the matlab (R14) arburg.  The modified arburg
     source is in the third attachment to this email.
     The "burg_model" function offers the addition of a "stopping criterion"
     argument which is intended to prevent addition of unecessary poles to the
     autoregressive model.  
 (b) The matlab (R14) pburg has a more subtle problem. It is OK for "whole" or
     "twosided" spectra, but when asked to do a "onesided" spectrum the
     first and last spectral-density values are only half what they ought
     to be.  It seems likely that the matlab code uses a periodogram
     function where it __should__ have calculated
     2.0*[squared amplitude of FFT].  According to the theory, the Burg
     spectrum is a continuous function of frequency (Equation 2.28,
     Kay & Marple, Proc IEEE, 69(11), 1981) and it does not have the
     end-discontinuities produced by the one-sided periodogram function.
     I conclude that, for one-sided spectra, the octave __power.m is
     mathematically correct and the matlab pburg.m is mathematically wrong.
     This is the most important reason for not exactly reproducing the
     behaviour of matlab pburg.
     When calculating a Burg spectrum, the FFT/DFT is simply a convenient
     numerical method.  There is no requirement in the theory to base the
     interface on the FFT.  The function "ar_psd" offers the option of
     achieving the same result by using the function "polyval".  Freqz
     also uses polyval.

2) The Welch functions

When I first wrote "welch_psd", I didn't know about "pwelch", so I just did
what seemed to make most sense.  I wanted it to be easy to use without
making mistakes, and I wanted the default arguments to give the best possible
spectra.  The Welch method is a bit difficult and it's easy to make mistakes.
When I first saw the current matlab interface to pwelch, it seemed to have
several design flaws:
  (a) the default number of data segments is 8.  This is far too small if we
      need a reliable spectrum of noisy data.  Matlab's new
      "psd(spectrum.welch,...)" doesn't do this; it sets a default window
      length of 64 samples and pads the FFT to 256 samples, giving 156
      segments per 10,000 samples, which is much better.
  (b) segment overlap should be expressed as a proportion or percentage
      of segment length (like an H/P spectrum analyser) rather than as a
      number of samples of overlap.  To specify the correct number of overlap
      samples the user must first know the segment length, which may not
      be available.  In "psd(spectrum.welch,...)" overlap is expressed
      as a percentage of segment length.
  (c) the sampling frequency argument (Fs) should be placed before the
      overlap and FFT-length arguments because Fs is nearly always used,
      but the overlap and FFT-length are rarely used.
      (octave pwelch has Fs in the correct position)
  (d) If the sampling frequency (Fs) argument is not given, units of
      returned frequency are radians/sec.  The observed effect is
      indistinguishable from a default Fs of 2*pi Hz, I had difficulty
      understanding this feature from the matlab documentation; it is
      a source of confusion that I would prefer to avoid.
Mathworks have addressed all of these problems except (d) in their newer
"psd(spectrum.welch,...)" function.

## Copyright (C) 1999 Paul Kienzle
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
##
## Modified 28 Sept 2006 by Peter Lanspeary for complex data.

## usage:  [a, v, k] = arburg (x, p)
##
## fits an AR (p)-model using Burg method (a so called maximum entropy model).
## x = data vector to estimate
## a: AR coefficients
## v: variance of white noise
## k: reflection coeffients for use in lattice filter
##
## The power spectrum of the resulting filter can be plotted with
## pburg(x, p), or you can plot it directly with power(sqrt(v), a).
##
## Example
##   ## Target system
##   pw=[0.2, 0.4, 0.45, 0.95];   # pole angle (nyquist freq. is 1.0)
##   pr=[0.98, 0.98, 0.98, 0.96]; # pole distance (0.0<=x<1.0)
##   sys_a = real(poly([pr, pr].*exp(1i*pi*[pw, -pw])));
##   order = 2*length(pw);
##   ## Filter impulse+random gaussian noise to produce signal
##   n = 1024;
##   s = [1 ; 0.1*randn(n-1,1)];
##   x = filter(1,sys_a,s); % AR system output
##   ## Determine system from signal
##   [a, v] = arburg(x, order);
##   ## Plot magnitude response of signal and matched system
##   figure(0);
##   mag = abs(fft(x))/sqrt(n);
##   [h, w] = freqz(sqrt(v), a, [], 2);
##   semilogy(2*[0:n/2-1]/n,mag(1:(n/2)),'1;spectrum;');
##   hold on;
##   semilogy(w,abs(h),sprintf('2;order %d burg;', order));
##   hold off;
##   ## Plot zero-pole graph of target system and matched system
##   figure(1);
##   axis("square"); __gnuplot_set__ pointsize 2; grid;
##   r = exp(2i*pi*[0:100]/100); plot(real(r), imag(r), "0;;");
##   hold on;
##   r = roots(sys_a); plot(real(r), imag(r), "1x;system;");
##   r = roots(a); plot(real(r), imag(r), "2x;arburg;");
##   hold off;
##   axis("normal"); __gnuplot_set__ pointsize 1; grid('off');
##
## See also:
## pburg, power, freqz, impz for measuring the characteristics
##    of the resulting filter
## aryule for alternative spectral estimators
##
## Note: Orphanidis '85 claims lattice filters are more tolerant of
## truncation errors, which is why you might want to use them.  However,
## lacking a lattice filter processor, I haven't tested that the lattice
## filter coefficients are reasonable.
##
## Algorithm derived from:
##    Sophocles J. Orfanidis (1985).
##    Optimum signal processing: An introduction.
##    New York: Macmillan.

function [a, v, k] = arburg (x, p)

  if (nargin != 2) usage("[a, v, k] = arburg(x,p)"); end

  k = zeros(1,p);
  n = length(x);
  x = reshape(x,1,[]);
  v = (x*x')/n;

  ## f and b are the forward and backward error sequences
  f = x(2:n);
  b = x(1:n-1);

  ## remaining stages i=2 to p
  for i=1:p

    ## get the i-th reflection coefficient
    g = (-2*f*b')/(f*f'+b*b');
    k(i) = g;

    ## generate next filter order
    if i==1
      a = [ g ] ;
    else
      a = [ g, a+g*conj(a(i-1:-1:1)) ];
    endif

    ## keep track of the error
    v = v*(1-g*conj(g));

    ## update the prediction error sequences
    oldf = f;
    f = oldf(2:n-i) + g*b(2:n-i);
    b = b(1:n-i-1) + conj(g)*oldf(1:n-i-1);

  endfor
  a = [ 1, a(p:-1:1) ] ;

endfunction

%!demo
%! % use demo('pburg');

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[dbateman3]

Peter Vernon Lanspeary wrote:

>
> I agree that it would be nice to be compatible with matlab, but I
> would prefer to avoid duplicating matlab bugs and badly-designed parts of the
> interface (see the second attachment).  For simple methods with few parameters,
> such as the Burg method, the best choice of interface is obvious and there is
> no problem.  The Welch method is much more complicated and tricky, and a good
> pwelch interface is more difficult to design.  Even Mathworks have not set the
> interface in concrete; their documentation (e.g.
> http://www.mathworks.com/access/helpdesk/help/toolbox/signal/index.html?/
> access/helpdesk/help/toolbox/signal/f12-6587.html)
> seems to steer users to the newer "psd(spectrum.welch,...)", which has
> differently structured arguments and different defaults.  A colleague says
> that it is much easier to use than pwelch. Pwelch might not seem so attractive.

Then its more important to try for psd compatibility than pwelch.

> I can do the following without too much trouble.
>
> 1) The Burg-function arguments (both my functions and octave's) are
>    compatible with matlab, so replacing the octave arburg and octave
>    pburg gives
>    a) an interface compatible with matlab
>    b) correct results with the following differences in behaviour
>       a) returned frequencies are always Hz, never rad/sec
>       b) does not duplicate the matlab R14 periodogram bug
>       c) for onesided spectra, returns nfft or length(nfft) PSD values
>       d) additional arguments give extra features.
>    However, fixing the octave arburg (see the third attachment) delivers all
>    of the above except 2)d).

If you consider the extra features are worth it, I'd say you are in the
best position to decide. Make a choice.

> 2) The arguments of welch_psd would be made more compatible with matlab pwelch,
>    but not __completely__ compatible.  If we replace "pad" with "nfft", where
>    nfft=length(window)+pad, and make "Fs" the 5th argument,
>       [Pxx,f]=welch_psd(x,window,Fs,overlap,pad,range,units,trend,sloppy)
>    becomes
>       [Pxx,f]=welch_psd(x,window,overlap,nfft,Fs,range,units,trend,sloppy)
>    Matlab pwelch is
>          [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,range)
>
>    I would then change the default window function to Hamming and change the
>    overlap from a fraction (e.g. 0.5) to a percentage (e.g. 50) so that it
>    becomes compatible with matlab's "psd(spectrum.welch,...)".
>    The remaining 3 differences between corresponding input arguments are:
>    +----------+--------------------+--------------------+--------------------+
>    | Argument |   matlab pwelch    |   welch_psd        |Is compatibility with
>    +----------+--------------------+--------------------+ pwelch possible ??
>    +----------+--------------------+--------------------+--------------------+
>    | window   | default_length=    | default_length=    | yes, but would     |
>    |          |    length(x)/8     |    sqrt(length(x)) | make bad spectra   |
>    |          |Is incompatible with| rounded up to 2^n. |                    |
>    |          | psd(spectrum.welch)|                    |                    |
>    +----------+--------------------+--------------------+--------------------+
>    | nfft     | length of FFT or   | =length(window)+pad| is compatible      |
>    |          +--------------------+--------------------+--------------------+
>    |          | list of frequencies|                    | NO - not available |
>    +----------+--------------------+--------------------+--------------------+
>    | overlap  | overlap samples    | overlap percentage |would be undesirable|
>    +----------+--------------------+--------------------+--------------------+
>    Also, returned frequencies would always be in Hz, never rad/sec
>
> I would be happy build new pwelch, arburg, ar_psd and pburg
>   a) if that is what you propose,
>   b) if the above details are acceptable to you, and
>   c) if I have access to the CVS.
> I have a sourceforge account, but I am not registered as an  octave-forge
> developer. I guess I'll need to "read the instructions" first...

Yes this is what I propose. What is your sourceforge ID and I'll add it.
As for the developer guide it has been completely written for the
package manager issues at the moment. So if you want the developers
guide I'd recommend you read at the moment, you'll have to get the CVS
and run

./configure; make; cd www; make

which will take a while for the last step due to the building of the
function references for the website. The file to look at is then
developers.html

D.

>
> Thanks,
> Peter Lanspeary
>
>
>
> ------------------------------------------------------------------------
>
> COMPARISON OF WELCH and BURG POWER-SPECTRUM FUNCTIONS
> P. V. Lanspeary
>
> Contents
> --------
> Section 1. Comparison of arguments. Summary.
> Section 2. Comparison of behaviours. Summary.
> Section 3. pburg: Full compatibility with Matlab ?
> Section 4. Comparison of arguments. Details.
> Section 5. Comparison of behaviours. Details.
>
>
> Section 1. Comparison of arguments. Summary.
> --------------------------------------------
> I have compared the arguments of welch_psd.m, burg_model.m and burg_psd.m
> with equivalent functions in octave-forge and Matlab-sigproc.
> This is a summary. See section (4) for details.
>
> (a) welch_psd with pwelch (octave) and pwelch (matlab):
>   The behaviour of welch_psd and pwelch is complicated, so details of the
>   comparison are messy and complicated. In summary, the news is not
>   good.  Overall functionality is similar (not the same). Only the
>   output args and first input arg are compatible.  
>
> (b) burg_model with arburg (octave) and arburg (matlab):
>   Arguments are compatible except that burg_model has an extra optional
>   input argument.  If each function is called with the same 2 arguments,
>   behaviour should (according to the documentation) be the same.
>
> (c) burg_psd with pburg (octave) and pburg (matlab):
>   Arguments are compatible except that Octave's pburg has one an extra
>   optional input argument and burg_psd has 3 extra optional input arguments.
>
>
> Section 2. Comparison of behaviours. Summary.
> ---------------------------------------------
> I have run _two_ test cases on the "Welch" functions and two more on
> the "Burg" functions.  One test case has a real signal (data); the
> other has a complex signal.  This is a summary.  See section (5) for
> details.
>
> (a) pwelch (octave), pwelch (matlab), welch_psd:
>   Arguments are incompatible; so to do the same thing with each function
>   requires different arguments in each case.  If each function is
>   called with the same arguments, behaviour (except in special cases) is
>   different.  Producing the same behaviour requires care. Where documentation
>   requirees the same behaviour, output is the same to an accuracy of no worse
>   than 10x machine precision (eps*10).
>
> (b) arburg (octave), arburg (matlab), burg_model:
>   For the real-data test case, results are exactly the same --- except that
>   reflection coefficients (k) from octave arburg have the wrong sign.
>   For complex data, the Matlab arburg gives the same (correct) answer as
>   burg_model, but the octave arburg does not give the correct answer.
>
> (c) pburg (octave), pburg (matlab), burg_psd:
>   For 'twosided' spectra Matlab pburg gives the same answers as burg_psd, but
>   for 'onesided' spectra, the first and last values in the spectrum are too
>   small by a factor of 2.
>   The octave pburg does not work for complex data because it calls the octave
>   arburg.  The octave __power.m works correctly for both real and complex
>   data.
>
>
> Section 3. pburg: Full compatibility with Matlab ?
> --------------------------------------------------
>
>>> What would it take to create a fully Matlab compatiable pburg function?
> Arguments are fully compatible, but with optional extra arguments.
> Compatible behaviour is harder to achieve, and I am not sure that this is
> desirable. For example, the Matlab pburg
> a) assumes units of frequency are radians/sec unless sampling
>    frequency, Fs, is specified.  If Fs is specified, the units are Hertz.
>    Octave pburg uses only Hertz.
> b) returns only nfft/2+1 (frequency,spectral-density) coordinates if the
>    spectrum is onesided.  Octave pburg returns nfft values.
> c) plots the spectral density in dB rather than correct physical units;
>    the graph should be a power spectrum of a physical quantity rather
>    than the response of a filter.
> d) has a error in onesided spectra (the first and last values of spectral
>    density are too small by a factor of 2).
> A completely fixed Octave pburg would seem to have more appropriate behaviour.
>
>
> Section 4. Comparison of arguments. Details.
> --------------------------------------------
>
> (a) pwelch (octave), pwelch (matlab), welch_psd:
>   The number of arguments is large and behaviour of the functions is a bit
>   complicated ... for the sake of brevity this comparison/description
>   therefore contains some minor approximations. First we write the first line
>   of each function, adjusting argument names to make them more consistent::
>
>   1) octave [Pxx,Pci,f]=pwelch(x,nfft,Fs,window,overlap,ci,range,units,trend)
>   2) matlab [Pxx,f]    =pwelch(x,window,overlap,nfft,Fs,range)
>   3)        [Pxx,f]=welch_psd(x,window,Fs,overlap,pad,range,units,trend,sloppy)
>
>   +----------+--------------------+--------------------+--------------------+
>   | Argument |  octave pwelch     |   matlab pwelch    |   welch_psd        |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   |  Pxx     | spectral density   | spectral density   | spectral density   |
>   |          | units specified by | ["x" units]/Hertz  | ["x" units]/Hertz  |
>   |          | "units" arg        | or per (rad/s)     |                    |
>   +----------+--------------------+--------------------+--------------------+
>   | Pci      | confidence interval| not available      | not available      |
>   +----------+--------------------+--------------------+--------------------+
>   | f        | frequency (Hertz)  | frequency          | frequency (Hertz)  |
>   |          |                    | (Hertz if given Fs)|                    |
>   |          |                    | (rad/s without Fs) |                    |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   | x        | complex data                                                 |
>   +----------+--------------------+--------------------+--------------------+
>   | window   | windowing vector, or the length of data segment              |
>   |          +--------------------+--------------------+--------------------+
>   |          | default = Hann     | default = Hamming  | default = Hann     |
>   |          | default_length=nfft| default_length=    | default_length=    |
>   |          |                    | length(x)/8        |    sqrt(length(x)) |
>   |          |                    |                    |   rounded up to 2^n|
>   +----------+--------------------+--------------------+--------------------+
>   | Fs       | sampling frequency | sampling frequency | sampling frequency |
>   |          | default=2 Hz       | default= 1 Hz      | default=1 Hz       |
>   +----------+--------------------+--------------------+--------------------+
>   | nfft     | length of FFT      | length of FFT or   | not available; uses|
>   |          | default = 256      | list of frequencies| length(window)+pad |
>   |          |                    | default=256 or 2^n |    as nfft         |
>   |          |                    |   >length(window)  |                    |
>   +----------+--------------------+--------------------+--------------------+
>   | overlap  | number of overlap  | number of overlap  | overlap fraction   |
>   |          | samples, default=  | samples, default=  | default=0.5, i.e.  |
>   |          | length(window)/2   | length(window)/2   | 50%                |
>   +----------+--------------------+--------------------+--------------------+
>   | pad      | not available      | not available      | length of zero     |
>   |          |                    |                    | padding on FFT     |
>   |          |=nfft-length(window)|=nfft-length(window)| default=0          |
>   +----------+--------------------+--------------------+--------------------+
>   | ci       | confidence interval| not available      | not available      |
>   +----------+--------------------+--------------------+--------------------+
>   | range    | 'half'             | 'half', 'onesided' | 'half', 'onesided' |
>   |          | 'whole'            | 'whole', 'twosided'| 'whole', 'twosided'|
>   |          |                    |                    | 'shift'            |
>   +----------+--------------------+--------------------+--------------------+
>   | units    | plot scaling       | not available      | 'plot', 'semilogx' |
>   |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
>   |          |                    |                    | 'loglog', 'squared'|
>   +----------+--------------------+--------------------+--------------------+
>   | trend    | remove trend       | not available      |remove trend, 'long'|
>   |          | 'mean', 'linear'   | -- does not        |'short','no-detrend'|
>   |          | default=none       |    remove mean     | default='long'     |
>   +----------+--------------------+--------------------+--------------------+
>   | sloppy   | not available      | not available      | round nfft to 2^n  |
>   +----------+--------------------+--------------------+--------------------+
>     The above table does not describe the special cases where the octave pwelch
>     and welch_psd are used for calculating cross-power spectra, coherence
>     and/or transfer functions.
>
> (b) arburg (octave), arburg (matlab), burg_model:
>   (1) octave  [a,v,k]=arburg(x,p)
>   (2) matlab  [a,v,k]=arburg(x,p)
>   (3)         [a,v,k]=burg_model(x,p,stop)
>
>   +----------+--------------------+--------------------+--------------------+
>   | Argument |  octave arburg     |   matlab arburg    |   burg_model       |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   | a        | autoregression coefficients                                  |
>   +----------+--------------------+--------------------+--------------------+
>   | v        | mean square residual (noise)                                 |
>   +----------+--------------------+--------------------+--------------------+
>   | k        | reflection coefficients                                      |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   | x        | real data          | complex data       | complex data       |
>   +----------+--------------------+--------------------+--------------------+
>   | p        | number of poles (required arg)                               |
>   +----------+--------------------+--------------------+--------------------+
>   | stop     | not available      | not available      | stopping criterion |
>   |          |                    |                    | 'FPE','AIC'        |
>   |          |                    |                    | default='none'     |
>   +----------+--------------------+--------------------+--------------------+
>
> (c) pburg (octave), pburg (matlab), burg_psd:
>   (1) octave  [Pxx,f]=pburg(x,p,nfft,Fs,range,units)
>   (2) matlab  [Pxx,f]=pburg(x,p,nfft,Fs,range)
>   (3)         [Pxx,f]=burg_psd(x,p,nfft,Fs,range,units,method,stop)
>
>   +----------+--------------------+--------------------+--------------------+
>   | Argument |  octave pburg      |   matlab pburg     |   burg_psd         |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   |  Pxx     | spectral density   | spectral density   | spectral density   |
>   |          | ["x" units]/Hertz  | ["x" units]/Hertz  | ["x" units]/Hertz  |
>   |          |                    | or per (rad/s)     |                    |
>   +----------+--------------------+--------------------+--------------------+
>   | f        | frequency (Hertz)  | frequency          | frequency (Hertz)  |
>   |          |                    | (Hertz if given Fs)|                    |
>   |          |                    | (rad/s without Fs) |                    |
>   +----------+--------------------+--------------------+--------------------+
>   +----------+--------------------+--------------------+--------------------+
>   | x        | real data          | complex data       | complex data       |
>   +----------+--------------------+--------------------+--------------------+
>   | p        | number of poles (required arg)                               |
>   +----------+--------------------+--------------------+--------------------+
>   | nfft     | number of frequency| length of FFT or   | number of frequency|
>   |          | values             | list of frequencies| values or list of  |
>   |          |                    |                    | frequencies        |
>   |          | default = 256      | default=256        | default=256        |
>   +----------+--------------------+--------------------+--------------------+
>   | Fs       | sampling frequency | sampling frequency | sampling frequency |
>   |          | default=2 Hz       | default= 1 Hz      | default=1 Hz       |
>   +----------+--------------------+--------------------+--------------------+
>   | range    | 'half'             | 'half', 'onesided' | 'half', 'onesided' |
>   |          | 'whole'            | 'whole', 'twosided'| 'whole', 'twosided'|
>   |          |                    |                    | 'shift'            |
>   +----------+--------------------+--------------------+--------------------+
>   | units    | plot scaling       | not available      | 'plot', 'semilogx' |
>   |          | 'squared', 'db'    |                    | 'semilogy', 'db'   |
>   |          |                    |                    | 'loglog', 'squared'|
>   +----------+--------------------+--------------------+--------------------+
>   | method   | not available      | not available      | 'FFT', 'poly'      |
>   +----------+--------------------+--------------------+--------------------+
>   | stop     | not available      | not available      | stopping criterion |
>   |          |                    |                    | 'FPE','AIC'        |
>   |          |                    |                    | default='none'     |
>   +----------+--------------------+--------------------+--------------------+
>
>
> Section 5. Comparison of behaviours. Details.
> ---------------------------------------------
> The test platforms are
>    a) Octave 2.1.69 + octave-forge-2006.07.09 on Debian 3.1r1
>       with sourceforge.net/tracker bugfix #1530283   for arburg.m
>       with sourceforge.net/tracker bugfix #1530309   for __power.m
>    b) Matlab 7.0.1.24704(R14)SP1 + sigproc toolbox v6.2.1, on MS-windows XP
>       (this is the most recent Matlab version available to me)
>
> (a) pwelch (octave), pwelch (matlab), welch_psd:
>   The test script is
>
>      format long
>      rand('seed',2038014164);
>      a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
>      white = rand(1,16384);
>      signal = detrend(filter(0.70181,a,white));
>      skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
>      %% REAL-DATA TEST
>     [psd1,f1]=welch_psd(signal,hann(256),25,0.5,0,'whole','plot','no-detrend');
>      if ( strcmp(version,'2.1.69') )
>        disp( 'OCTAVE 2.1.69' )
>        [psd2,f2]=pwelch(signal,256,25,hann(256),128,'whole','squared','none');
>      else
>        disp( 'MATLAB' )
>        [psd2,f2]=pwelch(signal,hann(256),256*0.5,256,25,'twosided');
>        end
>      sum_psd1=sum(psd1)
>      sum_psd2=sum(psd2)
>      maxdiffr=max(abs(psd1./reshape(psd2,1,[])-1))  %% max relative difference
>      %% COMPLEX-DATA TEST
>     [psd3,f3]=welch_psd(skewed,hann(256),25,0.5,0,'whole','plot','no-detrend');
>      if ( strcmp(version,'2.1.69') )
>        disp( 'OCTAVE 2.1.69' )
>        [psd4,f4]=pwelch(skewed,256,25,hann(256),128,'whole','squared','none');
>      else
>        disp( 'MATLAB' )
>        [psd4,f4]=pwelch(skewed,hann(256),256*0.5,256,25,'twosided');
>        end
>      sum_psd3=sum(psd3)
>      sum_psd4=sum(psd4)
>      maxdiffc=max(abs(psd3./psd4'-1))  %% maximum relative difference
>      eps
>
>   +----------+-----------------------+---------------------------+
>   |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
>   |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
>   +----------+-----------------------+---------------------------+
>   |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
>   | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
>   +----------+-----------------------+---------------------------+
>   | real data                                                    |
>   | sum_psd1 | 3.62735601946465      | 3.61040879075053          |
>   | sum_psd2 | 3.62735601946465      | 3.61040879075053          |
>   | maxdiffr | 1.33226762955019e-15  | 1.110223024625157e-015    |
>   +----------+-----------------------+---------------------------+
>   | complex data                                                 |
>   | sum_psd3 | 3.62735601946465      | 3.61040879075054          |
>   | sum_psd4 | 3.62735601946465      | 3.61040879075054          |
>   | maxdiffc | 1.33226762955019e-15  | 1.332267629550188e-015    |
>   +----------+-----------------------+---------------------------+
>     We should have sum_psd1==sum_psd2, maxdiffr==0, sum_psd3==sum_psd4,
>     maxdiffc==0.
>
>
> (b) arburg (octave), arburg (matlab), burg_model:
>   The additional test script (carrying on from previous test) is
>
>      %% REAL-DATA TEST
>      [ar1,vr1,kr1]=burg_model(signal,5);
>      [ar2,vr2,kr2]=arburg(signal,5);
>      maxdifar=max(abs(ar2./ar1-1))  %% maximum relative difference
>      vr1/vr2-1
>      maxdifkr=max(abs(reshape(kr2,1,[])./kr1-1)) %% kr2 has inconsistent shape
>      %% COMPLEX-DATA TEST
>      [ac1,vc1,kc1]=burg_model(skewed,5);
>      [ac2,vc2,kc2]=arburg(skewed,5);
>      maxdifac=max(abs(ac2./ac1-1))  %% maximum relative difference
>      vc2/vc1-1
>      maxdifkc=max(abs(reshape(kc2,1,[])./kc1-1))
>
>   +----------+-----------------------+---------------------------+
>   |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
>   |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
>   +----------+-----------------------+---------------------------+
>   |          | burg_model vs. arburg | burg_model vs. arburg     |
>   | eps      | 2.22044604925031e-16  |  2.220446049250313e-016   |
>   +----------+-----------------------+---------------------------+
>   | real data                                                    |
>   | maxdifar | 8.88178419700125e-16  |  0                        |
>   | vr1/vr2-1| -5.55111512312578e-16 |  0                        |
>   | maxdifkr | 2.00000000000000      |  0                        |
>   +----------+-----------------------+---------------------------+
>   | complex data                                                 |
>   | maxdifac | 2.15561046016939      |  0                        |
>   | vc2/vc1-1| 7.679016 - 0.0026252i |  0                        |
>   | maxdifkc | 1.68536747392439      |  5.590172498962309e-017   |
>   +----------+-----------------------+---------------------------+
>      We should have maxdifar==0, vr1/vr2-1==0, maxdifkr==0, maxdifac==0,
>      vc2/vc1-1==0, maxdifkc==0.
>      ... Octave arburg gives the correct autoregression coefficients and
>      variance (ar1,vr1) for real data but the reflection coefficients have
>      the wrong sign; for complex data, the octave arburg gives the wrong
>      answers.  Matlab pburg and burg_model give exactly the same answer
>      for both real and complex data.
>      Now compare __power with ar_psd using coefficients (ac1,vc1)
>      calculated by burg_model  ...
>
>      [psd3,f3]=ar_psd(ac1,vc1,256,25,'whole','squared');
>      [psd4,f4]=__power(sqrt(vc1),ac1,256,25,'whole','squared');
>      sum_psd3=sum(psd3)
>      %%     sum_psd3 = 3.63449149426036e+00 + 4.89110705062523e-19i
>      sum_psd4=sum(psd4)
>      %%     sum_psd4 = 3.63449149426036
>      maxdiffc=max(abs(psd3./psd4-1))
>      %%     maxdiffc =  3.33066934574812e-16
>
>      ... we conclude that __power works OK for complex filter coefficients
>
> (b) pburg (octave), pburg (matlab), burg_psd:
>   The additional test script (carrying on from previous test) is
>
>      %% REAL-DATA TEST
>      [psd5,f5]=burg_psd(signal,5,256,25,'whole');
>      if ( strcmp(version,'2.1.69') )
>        disp( 'OCTAVE 2.1.69' )
>        [psd6,f6]=pburg(signal,5,256,25,'whole');
>      else
>        disp( 'MATLAB' )
>        [psd6,f6]=pburg(signal,5,256,25,'twosided');
>        end
>      sum_psd5=sum(psd5)
>      sum_psd6=sum(psd6)
>      maxdiffr=max(abs(psd5./reshape(psd6,1,[])-1))  %% max relative difference
>      %% COMPLEX-DATA TEST
>      [psd7,f7]=burg_psd(skewed,5,256,25,'whole');
>      if ( strcmp(version,'2.1.69') )
>        disp( 'OCTAVE 2.1.69' )
>        [psd8,f8]=pburg(skewed,5,256,25,'whole');
>      else
>        disp( 'MATLAB' )
>        [psd8,f8]=pburg(skewed,5,256,25,'twosided');
>        end
>      sum_psd7=sum(psd7)
>      sum_psd8=sum(psd8)
>      maxdiffc=max(abs(psd7./reshape(psd8,1,[])-1))  %% max relative difference
>  
>   +----------+-----------------------+---------------------------+
>   |          | Octave 2.1.69         | Matlab 7.0.1.24704(R14)SP1|
>   |          | Debian3.1r1/Intel-P4  | MS-Windows-XP/Intel-P4    |
>   +----------+-----------------------+---------------------------+
>   |          | welch_psd vs. pwelch  | welch_psd vs. pwelch      |
>   | eps      | 2.22044604925031e-16  | 2.220446049250313e-016    |
>   +----------+-----------------------+---------------------------+
>   | real data                                                    |
>   | sum_psd5 | 3.63449149426037      | 3.62433951286238          |
>   | sum_psd6 | 3.63449149426037      | 3.62433951286238          |
>   | maxdiffr | 1.24344978758018e-14  | 5.551115123125783e-016    |
>   +----------+-----------------------+---------------------------+
>   | complex data                                                 |
>   | sum_psd7 | 3.634491494 +4.9e-19i | 3.62433951286238          |
>   | sum_psd8 | 3.64552153406807      | 3.62433951286238          |
>   | maxdiffc | 6.15692319331572      | 6.661338147750939e-016    |
>   +----------+-----------------------+---------------------------+
>      We should have sum_psd5==sum_psd6, maxdiffr==0, sum_psd7==sum_psd8,
>      maxdiffc==0.
>      Octave pburg doesn't work for complex data because it calls  arburg.
>      Twosided spectra from matlab pburg are OK, but for "onesided" spectrum
>       [psd8,f8]=pburg(skewed,5,256,25,'onesided');
>     the first and last values in the spectrum from matlab are too small by
>     a factor of 2.
>
>
> ------------------------------------------------------------------------
>
> Bugs and Design Problems with Burg and Welch functions
> ------------------------------------------------------
> Peter Lanspeary, 28 Sept 2006
>
> 1) The Burg functions.
>
> If we ignore additional arguments, both octave (arburg, pburg) and my efforts
> (burg_model, purg_psd) have interfaces which are compatible with matlab.
> However, behaviour is different.
>  (a) The octave arburg does not work with complex data, and with real data
>      the returned reflection coefficients have the wrong sign.  I have
>      modified the octave arburg code to fix these problems so that its
>      behaviour is identical to the matlab (R14) arburg.  The modified arburg
>      source is in the third attachment to this email.
>      The "burg_model" function offers the addition of a "stopping criterion"
>      argument which is intended to prevent addition of unecessary poles to the
>      autoregressive model.  
>  (b) The matlab (R14) pburg has a more subtle problem. It is OK for "whole" or
>      "twosided" spectra, but when asked to do a "onesided" spectrum the
>      first and last spectral-density values are only half what they ought
>      to be.  It seems likely that the matlab code uses a periodogram
>      function where it __should__ have calculated
>      2.0*[squared amplitude of FFT].  According to the theory, the Burg
>      spectrum is a continuous function of frequency (Equation 2.28,
>      Kay & Marple, Proc IEEE, 69(11), 1981) and it does not have the
>      end-discontinuities produced by the one-sided periodogram function.
>      I conclude that, for one-sided spectra, the octave __power.m is
>      mathematically correct and the matlab pburg.m is mathematically wrong.
>      This is the most important reason for not exactly reproducing the
>      behaviour of matlab pburg.
>      When calculating a Burg spectrum, the FFT/DFT is simply a convenient
>      numerical method.  There is no requirement in the theory to base the
>      interface on the FFT.  The function "ar_psd" offers the option of
>      achieving the same result by using the function "polyval".  Freqz
>      also uses polyval.
>
> 2) The Welch functions
>
> When I first wrote "welch_psd", I didn't know about "pwelch", so I just did
> what seemed to make most sense.  I wanted it to be easy to use without
> making mistakes, and I wanted the default arguments to give the best possible
> spectra.  The Welch method is a bit difficult and it's easy to make mistakes.
> When I first saw the current matlab interface to pwelch, it seemed to have
> several design flaws:
>   (a) the default number of data segments is 8.  This is far too small if we
>       need a reliable spectrum of noisy data.  Matlab's new
>       "psd(spectrum.welch,...)" doesn't do this; it sets a default window
>       length of 64 samples and pads the FFT to 256 samples, giving 156
>       segments per 10,000 samples, which is much better.
>   (b) segment overlap should be expressed as a proportion or percentage
>       of segment length (like an H/P spectrum analyser) rather than as a
>       number of samples of overlap.  To specify the correct number of overlap
>       samples the user must first know the segment length, which may not
>       be available.  In "psd(spectrum.welch,...)" overlap is expressed
>       as a percentage of segment length.
>   (c) the sampling frequency argument (Fs) should be placed before the
>       overlap and FFT-length arguments because Fs is nearly always used,
>       but the overlap and FFT-length are rarely used.
>       (octave pwelch has Fs in the correct position)
>   (d) If the sampling frequency (Fs) argument is not given, units of
>       returned frequency are radians/sec.  The observed effect is
>       indistinguishable from a default Fs of 2*pi Hz, I had difficulty
>       understanding this feature from the matlab documentation; it is
>       a source of confusion that I would prefer to avoid.
> Mathworks have addressed all of these problems except (d) in their newer
> "psd(spectrum.welch,...)" function.
>
>
> ------------------------------------------------------------------------
>
> ## Copyright (C) 1999 Paul Kienzle
> ##
> ## This program is free software; you can redistribute it and/or modify
> ## it under the terms of the GNU General Public License as published by
> ## the Free Software Foundation; either version 2 of the License, or
> ## (at your option) any later version.
> ##
> ## This program is distributed in the hope that it will be useful,
> ## but WITHOUT ANY WARRANTY; without even the implied warranty of
> ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
> ## GNU General Public License for more details.
> ##
> ## You should have received a copy of the GNU General Public License
> ## along with this program; if not, write to the Free Software
> ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
> ##
> ## Modified 28 Sept 2006 by Peter Lanspeary for complex data.
>
> ## usage:  [a, v, k] = arburg (x, p)
> ##
> ## fits an AR (p)-model using Burg method (a so called maximum entropy model).
> ## x = data vector to estimate
> ## a: AR coefficients
> ## v: variance of white noise
> ## k: reflection coeffients for use in lattice filter
> ##
> ## The power spectrum of the resulting filter can be plotted with
> ## pburg(x, p), or you can plot it directly with power(sqrt(v), a).
> ##
> ## Example
> ##   ## Target system
> ##   pw=[0.2, 0.4, 0.45, 0.95];   # pole angle (nyquist freq. is 1.0)
> ##   pr=[0.98, 0.98, 0.98, 0.96]; # pole distance (0.0<=x<1.0)
> ##   sys_a = real(poly([pr, pr].*exp(1i*pi*[pw, -pw])));
> ##   order = 2*length(pw);
> ##   ## Filter impulse+random gaussian noise to produce signal
> ##   n = 1024;
> ##   s = [1 ; 0.1*randn(n-1,1)];
> ##   x = filter(1,sys_a,s); % AR system output
> ##   ## Determine system from signal
> ##   [a, v] = arburg(x, order);
> ##   ## Plot magnitude response of signal and matched system
> ##   figure(0);
> ##   mag = abs(fft(x))/sqrt(n);
> ##   [h, w] = freqz(sqrt(v), a, [], 2);
> ##   semilogy(2*[0:n/2-1]/n,mag(1:(n/2)),'1;spectrum;');
> ##   hold on;
> ##   semilogy(w,abs(h),sprintf('2;order %d burg;', order));
> ##   hold off;
> ##   ## Plot zero-pole graph of target system and matched system
> ##   figure(1);
> ##   axis("square"); __gnuplot_set__ pointsize 2; grid;
> ##   r = exp(2i*pi*[0:100]/100); plot(real(r), imag(r), "0;;");
> ##   hold on;
> ##   r = roots(sys_a); plot(real(r), imag(r), "1x;system;");
> ##   r = roots(a); plot(real(r), imag(r), "2x;arburg;");
> ##   hold off;
> ##   axis("normal"); __gnuplot_set__ pointsize 1; grid('off');
> ##
> ## See also:
> ## pburg, power, freqz, impz for measuring the characteristics
> ##    of the resulting filter
> ## aryule for alternative spectral estimators
> ##
> ## Note: Orphanidis '85 claims lattice filters are more tolerant of
> ## truncation errors, which is why you might want to use them.  However,
> ## lacking a lattice filter processor, I haven't tested that the lattice
> ## filter coefficients are reasonable.
> ##
> ## Algorithm derived from:
> ##    Sophocles J. Orfanidis (1985).
> ##    Optimum signal processing: An introduction.
> ##    New York: Macmillan.
>
> function [a, v, k] = arburg (x, p)
>
>   if (nargin != 2) usage("[a, v, k] = arburg(x,p)"); end
>
>   k = zeros(1,p);
>   n = length(x);
>   x = reshape(x,1,[]);
>   v = (x*x')/n;
>
>   ## f and b are the forward and backward error sequences
>   f = x(2:n);
>   b = x(1:n-1);
>
>   ## remaining stages i=2 to p
>   for i=1:p
>
>     ## get the i-th reflection coefficient
>     g = (-2*f*b')/(f*f'+b*b');
>     k(i) = g;
>
>     ## generate next filter order
>     if i==1
>       a = [ g ] ;
>     else
>       a = [ g, a+g*conj(a(i-1:-1:1)) ];
>     endif
>
>     ## keep track of the error
>     v = v*(1-g*conj(g));
>
>     ## update the prediction error sequences
>     oldf = f;
>     f = oldf(2:n-i) + g*b(2:n-i);
>     b = b(1:n-i-1) + conj(g)*oldf(1:n-i-1);
>
>   endfor
>   a = [ 1, a(p:-1:1) ] ;
>
> endfunction
>
> %!demo
> %! % use demo('pburg');
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Octave-sources mailing list
> Octave-sources@...
> https://www.cae.wisc.edu/mailman/listinfo/octave-sources

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[Peter Lanspeary]

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
_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[John W. Eaton]

On  1-Nov-2006, Peter Vernon Lanspeary wrote:

| Hi everybody,

The sources list is really for posting source code.  There aren't that
many subscribers.  For discussions, the help list is generally best.
If you want to discuss future enhancements to the core of Octave, use
the maintainers list.  If you want to report a bug, use the bug list.

Thanks,

jwe
_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[dbateman3]

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}))

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources