Burg and Welch power spectra
Peter Vernon Lanspeary <[email protected]> Thu, 28 Sep 2006 19:00:27 +0930
| Newsgroups | gmane.comp.gnu.octave.sources |
|---|---|
| Message-ID | <[email protected]> |
--hHWLQfXTYDoKhP50
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
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
--hHWLQfXTYDoKhP50
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="comparison_Welch+Burg_fns.txt"
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.
--hHWLQfXTYDoKhP50
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename=tmp
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.
--hHWLQfXTYDoKhP50
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="arburg.m"
## 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');
--hHWLQfXTYDoKhP50
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources
--hHWLQfXTYDoKhP50--