View Raw SPL
#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);
}