Re: Burg and Welch power spectra
David Bateman <[email protected]> Thu, 28 Sep 2006 23:57:56 +0200
| Newsgroups | gmane.comp.gnu.octave.sources |
|---|---|
| Organization | Motorola CRM |
| Message-ID | <[email protected]> |
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
> [email protected]
> https://www.cae.wisc.edu/mailman/listinfo/octave-sources