/*****************************************************************************
*                                                                            *
*   FINTEG.SPL   Copyright (C) 2007-2012 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    frequency domain integration                                *
*                                                                            *
*   Revisions:   23 May 2007  RRR  Creation                                  *
*                10 Nov 2010  RRR  integration order                         *
*                 9 Mar 2012  RRR  optional high pass filter                 *
*                                                                            *
*****************************************************************************/

#if @HELP_FINTEG

    FINTEG

    Purpose: Performs frequency domain integration.

    Syntax:  FINTEG(s, n, fc)

              s  - the input series

              n  - Optional. An integer, the integration order. Defaults
                   to 1, first order integration.

              fc - Optional. A real, the high pass cutoff frequency. Defaults
                   to -1, no high pass filtering.

    Returns: A series or array.

    Example:
             W1: gnorm(1000, 1)
             W2: integ(W1)
             W3: finteg(W1)

             W2 contains a standard time domain integration of W1.
             W3 contains the integration as performed in the frequency
             domain.

    Example:
             W1: gnorm(10000, 1/1000);setvunits("G")
             W2: finteg(W1, 2, 5)

             W1 contains data sampled at 1000 Hz. W2 high pass filters the
             data with a cutoff frequency of 5 Hz and double integrates the
             result.

    Example:
             n = 1000;
             x = extract(linspace(0, 2*pi, n), 1, n-1);

             W1: -sin(x + 2*cos(3*x)) * (1-6*sin(3*x))
             W2: cos(x + 2*cos(3*x))
             W3: finteg(W1)
             W4: abs(W2 - W3)

             W1 contains 999 samples of the periodic function 

             -sin(x + 2*cos(3*x)) * (1-6*sin(3*x))

             W2 contains the analytic integral of W1 

             cos(x + 2*cos(3*x))

             W3 computes the frequency domain integration and W4 displays the
             difference between the analytic and calculated. In this case,
             the difference is negligible.
            
    Remarks:
             FINTEG performs integration in the frequency domain
             by calculating:

                  F(f)/(2*pi*f*i)

            where:     F(f) is the FFT of the original series
                       f    is the frequency range of the fft
                       i    = sqrt(-1)

            F(1) is set to 0 to demean the result. The integration
            calculation is inverse transformed back into the time
            domain.

            This method performs a de-meaned integration, i.e the DC
            offset of both the input and output series are removed.

            The result is converted back into the time domain.

            To view the result in the frequency domain, try:

                spectrum(finteg(s))

            finteg(s, n) computes the nth order integration with:

                  F(f) / (2*pi*f*i)^n

            where n is a positive integer.

            If fc >= 0, the input data is high passed filtered. The filter
            is implemented by zeroing out the frequency domain components
            of the input <= fc. See ZHPFILT for more details.

    See Also:
             Deriv
             Fderiv
             FFT
             Finteg
             Integ
             Spectrum
             Zhpfilt
#endif


/* frequency domain integration */
ITERATE finteg(s, n, fc)
{
        local len, df, f, y, w, u, j;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("finteg - input series required");
                
                        n = 1;
                }

                fc = -1;
        }

        /* length and delta F */
        len = length(s);
        df  = rate(s) / len;

        /* generate frequency values */
        if (len % 2 == 0)
        {
                /* even */
                f = df * concat(0..len/2, (-len/2 + 1)..-1);
        }
        else
        {
                /* odd */
                f = df * concat(0..(len-1)/2, -(len-1)/2..-1);
        }

        /* set frequency spacing */
        f.deltax = df;

        /* frequency scaling */
        w = (2 * pi * f * i);
        
        if (n != 1) w = w ^ n;

        /* frequency domain integration */
        y = fft(s) / w;

        /* high pass FFT filter */
        if (fc >= 0)
        {
                /* multiply by frequency domain high pass mask */
                y *= zhpfilt(s, fc, 1);
        }

        /* demean - also removes 1/0 pathology */
        y[1] = 0;

        /* inverse transform */
        y = ifft(y);

        /* convert to real */
        if (isreal(s))
        {
                y = real(y);
        }

        /* units */
        u = s[1..2];
        
        loop (j = 1..n)
        {
                u = integ(u);
        }
        
        y.vunits = u.vunits;

        /* offset */
        y.xoffset = s.xoffset;

        return(y);
}