Reply on Burg and Welch power spectra
Peter Vernon Lanspeary <[email protected]> Tue, 26 Sep 2006 00:58:35 +0930
| Newsgroups | gmane.comp.gnu.octave.sources |
|---|---|
| Message-ID | <[email protected]> |
--MGYHOYXEY6WxJCY8 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline > 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 [email protected] >Motorola Labs - Paris +33 1 69 35 48 04 (Ph) >Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob) >91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (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 --MGYHOYXEY6WxJCY8 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 (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. --MGYHOYXEY6WxJCY8 Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="ar_psd.m" %% 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 --MGYHOYXEY6WxJCY8 Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="burg_model.m" %% 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 --MGYHOYXEY6WxJCY8 Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="burg_psd.m" %% 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 --MGYHOYXEY6WxJCY8 Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="welch_psd.m" %% 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 elseif ( strcmp(arg, 'coher' ) ) if ( ~do_coher ) n_results = n_results+1; do_coher = n_results; end elseif ( strcmp(arg, 'ypower' ) ) if ( ~do_ypower ) n_results = n_results+1; do_ypower = n_results; end else error( strcat( 'welch_psd: string arg error: illegal value: ', arg ) ); return; end end %% %% default action is to calculate power spectrum if ( ~n_results ) n_results = 1; do_power = 1; end need_Pxx = do_power || do_trans || do_coher; need_Pxy = do_cross || do_trans || do_coher; need_Pyy = do_coher || do_ypower; %% %% end of decoding and checking optional arguments. Phew! %% %% calculate/adjust segment length, window function log_two = log(2); nearly_one = 0.99999999999; if ( is_win <= 1 ) if ( is_win ) seg_len = window; else seg_len = 2 ^ ceil( log( sqrt(x_len) ) * nearly_one / log_two ); end xx = seg_len - 1; window = 0.5 - 0.5 * cos( 2 * pi / xx * [0 : xx] ); win_meansq = (window * window.') / seg_len; else seg_len = length(window); win_meansq = (window * window.') / seg_len; end if ( x_len < seg_len ) error( 'welch_psd: error: "x" is shorter than segment/window' ); return; end overlap = fix(seg_len * overlap); %% %% calculate FFT length fft_len = seg_len + padding; if ( is_sloppy ) fft_len = 2 ^ ceil( log( fft_len ) * nearly_one / log_two ); end %% %% Remove mean from the data if ( rm_mean == 2 ) n_fft = floor(x_len/seg_len)*seg_len; if ( need_Pxx || need_Pxy ) avg_x = sum( x(1:n_fft) ) / n_fft; end if ( need_y || need_Pxy) avg_y = sum( y(1:n_fft) ) / n_fft; end elseif ( rm_mean == 0 ) avg_x = 0; avg_y = 0; end %% %% Calculate periodograms xx = zeros(1,fft_len); yy = xx; Pxx = xx; Pxy = xx; Pyy = xx; n_fft = 0; for start_seg = [1:seg_len-overlap:x_len-seg_len+1] if ( need_Pxx || need_Pxy ) if ( rm_mean == 1 ) avg_x = sum( x(start_seg:start_seg+seg_len-1) ) / seg_len; end xx(1:seg_len) = window .* ( x(start_seg:start_seg+seg_len-1)-avg_x ); fft_x = fft(xx); end if ( need_Pxy || need_Pyy ) if ( rm_mean == 1 ) avg_y = sum( y(start_seg:start_seg+seg_len-1) ) / seg_len; end yy(1:seg_len) = window .* ( y(start_seg:start_seg+seg_len-1) - avg_y ); fft_y = fft(yy); end if ( need_Pxx ) Pxx = Pxx + fft_x .* conj(fft_x); end if ( need_Pxy ) Pxy = Pxy + fft_x .* conj(fft_y); end if ( need_Pyy ) Pyy = Pyy + fft_y .* conj(fft_y); end n_fft = n_fft +1; end %% %% For one-sided spectra, contributions from negative frequencies are added %% to the positive side of the spectrum -- but not at zero or Nyquist %% (half sampling) frequencies. This keeps power equal in time and spectral %% domains. %% if ( range == 0 ) if ( ~ rem(fft_len,2) ) %% one-sided, fft_len is even psd_len = fft_len/2+1; if ( need_Pxx ) Pxx = Pxx(1:psd_len) + [0 Pxx(fft_len:-1:psd_len+1) 0]; end if ( need_Pxy ) Pxy = Pxy(1:psd_len) + conj([0 Pxy(fft_len:-1:psd_len+1) 0]); end if ( need_Pyy ) Pyy = Pyy(1:psd_len) + [0 Pyy(fft_len:-1:psd_len+1) 0]; end else %% one-sided, fft_len is odd psd_len = (fft_len+1)/2; if ( need_Pxx ) Pxx = Pxx(1:psd_len) + [0 Pxx(fft_len:-1:psd_len+1)]; end if ( need_Pxy ) Pxy = Pxy(1:psd_len) + conj([0 Pxy(fft_len:-1:psd_len+1)]); end if ( need_Pyy ) Pyy = Pyy(1:psd_len) + [0 Pyy(fft_len:-1:psd_len+1)]; end end else %% two-sided (and shifted) psd_len = fft_len; end %% %% Put all results in matrix, one row per spectrum spectra = zeros(n_results,psd_len); plot_title = zeros(n_results,17); scale = n_fft * seg_len * sample_f * win_meansq; if ( do_power ) spectra(do_power,:) = Pxx / scale; plot_title(do_power,:) = str2mat('power spectrum x '); end if ( do_cross ) spectra(do_cross,:) = Pxy / scale; plot_title(do_cross,:) = str2mat('cross spectrum '); end if ( do_trans ) spectra(do_trans,:) = Pxy ./ Pxx; plot_title(do_trans,:) = str2mat('transfer function'); end if ( do_coher ) spectra(do_coher,:) = Pxy .* conj(Pxy) ./ Pxx ./ Pyy; plot_title(do_coher,:) = str2mat('coherence '); end if ( do_ypower ) spectra(do_ypower,:) = Pyy / scale; plot_title(do_ypower,:) = str2mat('power spectrum y '); end freq = [0:psd_len-1] * ( sample_f / fft_len ); %% %% range='shift': Shift zero-frequency to the middle if ( range == 2 ) len2 = fix((fft_len+1)/2); spectra = [ spectra(:,len2+1:fft_len), spectra(:,1:len2)]; freq = [ freq(len2+1:fft_len)-sample_f, freq(1:len2)]; end %% %% Plot the spectrum if there are no return variables. %% Only amplitudes are plotted. Power and coherence have no phase. %% Phase of cross spectrum and transfer function are not plotted. if ( nargout >= 2 ) varargout{1} = spectra; varargout{2} = freq; elseif ( nargout == 1 ) varargout{1} = spectra; else for ii = 1: n_results if ( n_results > 1 ) figure(); end if ( plot_type == 1 ) plot(freq,abs(spectra(ii,:))); elseif ( plot_type == 2 ) semilogx(freq,abs(spectra(ii,:))); elseif ( plot_type == 3 ) semilogy(freq,abs(spectra(ii,:))); elseif ( plot_type == 4 ) loglog(freq,abs(spectra(ii,:))); elseif ( plot_type == 5 ) %% db plot(freq,10*log10(abs(spectra(ii,:)))); end title( char(plot_title(ii,:)) ); end end end --MGYHOYXEY6WxJCY8 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 --MGYHOYXEY6WxJCY8--