diff options
| -rw-r--r-- | macros/pburg.sci | 215 |
1 files changed, 130 insertions, 85 deletions
diff --git a/macros/pburg.sci b/macros/pburg.sci index 78198ed..f677157 100644 --- a/macros/pburg.sci +++ b/macros/pburg.sci @@ -1,93 +1,138 @@ -function [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion) +function [psd,f_out] = pburg(x, poles, varargin) //Calculate Burg maximum-entropy power spectral density. -//Calling Sequence -//[psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion) -//All but the first two arguments are optional and may be empty. -//Parameters -// x: [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] -// Fs: [real scalar] sampling frequency (Hertz) [default=1] -// range: 'half', 'onesided' : frequency range of the spectrum is from 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', 'centerdc' : 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" are real, the default range is 'half', otherwise default range is 'whole'. -// method: 'fft': use FFT to calculate power spectral density. 'poly': calculate spectral density 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)". This argument is ignored and a spectrum is not plotted if the caller requires a returned value. -// criterion: [optional string arg] model-selection criterion. Limits the number of poles so that spurious poles are not added when the whitened data has no more information in it (see Kay & Marple, 1981). Recognized values are 'AKICc' -- approximate corrected Kullback information criterion (recommended), 'KIC' -- Kullback information criterion 'AICc' -- corrected Akaike information criterion 'AIC' -- Akaike information criterion 'FPE' -- final prediction error" criterion The default is to NOT use a model-selection criterion. -//Description -//This function is being called from Octave + +//Calling Sequence: +//[psd,f_out] = pburg(x, poles, freq, Fs, range, method, plot_type, criterion) + +//Parameters: +//All but the first two parameters are optional and may be empty. +// +//x- [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] +//Fs:[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 from 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', 'centerdc'- same as 'whole' but with the first half of the spectrum swapped with second half to put the +//zero-frequency value in the middle. If "freq" is vector, 'shift' is ignored. If model coefficients "ar_coeffs" are real, the default +//range is 'half', otherwise default range is 'whole'. +// +//method: +//'fft'- use FFT to calculate power spectral density. +//'poly'- calculate spectral density 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)". This argument is ignored and a +//spectrum is not plotted if the caller requires a returned value. +// +//criterion: [optional string arg] model-selection criterion. Limits the number of poles so that spurious poles are not +//added when the whitened data has no more information in it (see Kay & Marple, 1981). Recognized values are- +//'AKICc' -- approximate corrected Kullback information criterion (recommended), +//'KIC' -- Kullback information criterion +//'AICc' -- corrected Akaike information criterion +//'AIC' -- Akaike information criterion +//'FPE' -- final prediction error" criterion +// +//The default is to NOT use a model-selection criterion +// +// RETURNED VALUES: +//If return values are not required by the caller, the spectrum is plotted and nothing is returned. +//psd: [real vector] power-spectral density estimate. +//f_out: [real vector] frequency values. + +//Description: //This function is a wrapper for arburg and ar_psd. -//The functions "arburg" and "ar_psd" do all the work. -//See "help arburg" and "help ar_psd" for further details. -//Examples -//a = [1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746]; -//[psd,f_out] = pburg(a,2); -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 2 | rhs > 8) -error("Wrong number of input arguments.") -end +//Examples: +//a = [1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746]; +//[psd, f_out] = pburg(a, 2); -select(rhs) - - case 2 then - if(lhs==1) - psd = callOctave("pburg",x,poless) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles) - else - error("Wrong number of output argments.") - end + funcprot(0); + if (nargin < 2) + error('pburg: need at least 2 args.'); + end + nvarargin = length(varargin); + criterion = []; + for iarg = 1:nvarargin + arrgh = varargin(iarg); + if (type(arrgh) == 10 && ( ~strcmp(arrgh,'AKICc') ||... + ~strcmp(arrgh,'KIC') || ~strcmp(arrgh,'AICc') ||... + ~strcmp(arrgh,'AIC') || ~strcmp(arrgh,'FPE') ) ) + criterion = arrgh; + if (nvarargin > 1) + varargin(iarg) = []; + else + varargin = list(); + end + end + end + [ar_coeffs,residual] = arburg(x,poles,criterion); + if (nargout == 0) + ar_psd(ar_coeffs,residual,varargin(:)); + elseif (nargout == 1) + psd = ar_psd(ar_coeffs,residual,varargin(:)); + elseif (nargout >= 2) + [psd,f_out] = ar_psd(ar_coeffs,residual,varargin(:)); + end - case 3 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq) - else - error("Wrong number of output argments.") - end - case 4 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq,Fs) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq,Fs) - else - error("Wrong number of output argments.") - end - case 5 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq,Fs,range) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range) - else - error("Wrong number of output argments.") - end - case 6 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq,Fs,range,method) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method) - else - error("Wrong number of output argments.") - end - case 7 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type) - else - error("Wrong number of output argments.") - end - case 8 then - if(lhs==1) - psd = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type,criterion) - elseif(lhs==2) - [psd,f_out] = callOctave("pburg",x,poles,freq,Fs,range,method,plot_type,criterion) - else - error("Wrong number of output argments.") - end - end endfunction +//tests: + +//fs = 1000; +//t = 0:1/fs:1-1/fs; +//x = cos(2*%pi*100*t); +//order = 4; +//[pxx, f] = pburg(x, order, [], fs); +//figure; +//plot(f, 10*log10(pxx)); +//title('PSD Estimate using Burg Method - Sinusoidal Signal'); +//xlabel('Frequency (Hz)'); +//ylabel('Power/Frequency (dB/Hz)'); + +//fs = 1000; +//t = 0:1/fs:1-1/fs; +//x = cos(2*%pi*100*t); +//orders = [2, 4, 8, 16]; +//figure; +//for i = 1:length(orders) +// order = orders(i); +// [pxx, f] = pburg(x, order, [], fs); +// subplot(length(orders), 1, i); +// plot(f, 10*log10(pxx)); +// title(['PSD Estimate using Burg Method - Order ' string(order)]); +// xlabel('Frequency (Hz)'); +// ylabel('Power/Frequency (dB/Hz)'); +//end + +//fs = 1000; +//t = 0:1/fs:0.1-1/fs; +//x = cos(2*%pi*100*t); +//order = 4; +//[pxx, f] = pburg(x, order, [], fs); +//figure; +//plot(f, 10*log10(pxx)); +//title('PSD Estimate using Burg Method - Short Data Segment'); +//xlabel('Frequency (Hz)'); +//ylabel('Power/Frequency (dB/Hz)'); +//fs = 1000; +//t = 0:1/fs:1-1/fs; +//x = cos(2*%pi*100*t) + cos(2*%pi*200*t); +//order = 4; +//[pxx, f] = pburg(x, order, [], fs); +//figure; +//plot(f, 10*log10(pxx)); +//title('PSD Estimate using Burg Method - Multicomponent Signal'); +//xlabel('Frequency (Hz)'); +//ylabel('Power/Frequency (dB/Hz)'); |