#define BANDFREQ_POWER 1
#define BANDFREQ_MEAN 2
#define BANDFREQ_MEDIAN 3
/* process args for bandpower, meanfreq, medfreq - returns Pxx and possibly frange */
freqband_process_args(func, s, f, frange, psdflag)
{
local pxx, fmin, fmax, dx, ferr, idx1, idx2, len, freqs;
/* get args */
(s, f, frange, psdflag) = freqband_parse_args(func, s, f, frange, psdflag);
if (isempty(s))
{
return({}, {});
}
/* function as numeric type */
type = freqband_process_type_from_func(func);
/* process args */
if (psdflag || length(f) > 2)
{
pxx = refseries(s);
if (psdflag == 1)
{
/* convert psd to power spectrum */
pxx = pxx * deltax(pxx);
}
if (not(isarray(f)))
{
if (isscalar(f) && f > 0)
{
/* delta F of pxx */
f = {0, 1.0 / f};
}
else
{
error(sprintf("%s - PSD frequency values must be a series", func));
}
}
fmin = min(f);
fmax = max(f);
if (not(all(f == xvals(pxx))))
{
/* update power spectrum to use provided frequencies */
dx = abs(f[2] - f[1]);
/* make sure we are using a copy */
pxx = pxx * 1.0;
setdeltax(pxx, dx);
(fmin, fmax) = xminmax(pxx);
}
}
else
{
if ((length(frange) > 0) && (type == BANDFREQ_POWER))
{
/* bandpower - use hamming window with mean squared correction */
pxx = powspec(hamming(s, 3));
}
else
{
pxx = powspec(s);
}
/* update delta F */
setdeltax(pxx, deltax(pxx) * f / rate(s));
/* min and max frequencies */
(fmin, fmax) = xminmax(pxx);
}
if (length(frange) > 0)
{
if (length(frange) == 2)
{
ferr = (frange[1] < fmin || frange[1] > (fmax + eps(fmax)) ||
frange[2] < fmin || frange[2] > (fmax + eps(fmax)));
}
else
{
ferr = 1;
}
if (ferr)
{
error(sprintf("%s - FRANGE must be a two element frequency series between %g and %g", func, fmin, fmax));
}
if (type == BANDFREQ_MEDIAN)
{
/* medfreq frrguency range */
return(pxx, frange);
}
else
{
/* Pxx frequencies for bandpower and meanfreq */
freqs = xvals(pxx);
/* convert frequencies to indices */
idx1 = find(freqs >= frange[1], 1);
idx2 = find(frange[2] <= freqs, 1);
idx1 = idx1[1];
idx2 = idx2[1];
len = idx2 - idx1 + 1;
/* extract range of interest */
pxx = extract(pxx, idx1, len);
return(pxx);
}
}
return(pxx, frange);
}
/* type integer from function string */
freqband_process_type_from_func(func)
{
local type;
func = tolower(func);
if (func == "bandpower")
{
type = BANDFREQ_POWER;
}
else if (func == "meanfreq")
{
type = BANDFREQ_MEAN;
}
else
{
type = BANDFREQ_MEDIAN;
}
return(type);
}
/* parse arguments */
freqband_parse_args(func, x, f, frange, psdflag)
{
if (argc < 5)
{
if (argc < 4)
{
if (argc < 3)
{
if (argc < 2 || not(isarray(x)))
{
error(sprintf("%s - input series required", func));
}
f = {};
}
frange = {};
}
psdflag = "";
}
if (length(x) == 0)
{
return({});
}
if (isempty(f))
{
f = rate(x);
}
if (isstring(frange))
{
psdflag = frange;
frange = {};
}
if (isstring(f))
{
psdflag = f;
f = rate(x);
}
if (not(isstring(psdflag)))
{
error(sprintf("%s - PSDFLAG must be a string", func));
}
else if (strlen(psdflag) > 0)
{
if (tolower(psdflag) == "psd")
{
psdflag = 1;
}
else if (tolower(psdflag) == "powspec")
{
psdflag = 2;
}
else
{
error(sprintf("%s - unknown flag '%s', PSDFLAG must 'psd' or 'powspec' if specified", func, psdflag));
}
}
else
{
psdflag = 0;
}
if (psdflag)
{
if (argc < 4)
{
if (argc < 3)
{
f = xvals(x);
}
else
{
if (length(f) == 2)
{
frange = f;
f = xvals(x);
}
else
{
frange = {};
}
}
}
else if (isarray(frange) && isarray(f))
{
if (length(frange) == 0 && length(f) == 2)
{
frange = f;
f = rate(x);
}
}
}
else
{
if (argc < 4)
{
if (argc < 3)
{
if (argc < 2)
{
f = rate(x);
}
if (isarray(f))
{
frange = f;
f = rate(x);
}
else
{
frange = {};
}
}
else
{
if (length(f) == 2)
{
frange = f;
f = rate(x);
}
else
{
frange = {};
}
}
}
}
if (not(isarray(frange)))
{
error(sprintf("%s - FRANGE must be a two element series of frequency values", func));
}
return(x, f, frange, psdflag);
}