/*****************************************************************************
*                                                                            *
*   FDERIV.SPL   Copyright (C) 2010 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    frequency domain derivative                                 *
*                                                                            *
*   Revisions:   10 Nov 2010  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_FDERIV

    FDERIV

    Purpose: Calculates the frequency domain derivative.

    Syntax:  FDERIV(s, n)

              s  - A series, the input.

              n  - Optional, an integer. The derivative order. Defaults
                   to 1, first derivative.

    Returns: A series or array.

    Example:
             W1: demean(integ(gnorm(1000, 1)))
             W2: finteg(W1)
             W3: deriv(W2)
             W4: fderiv(W2)

             W1 contains zero mean synthesized data that is frequency
             domain integrated W2.  W3 contains a standard time domain
             derivative of W1.  W3 contains the derivative as performed
             in the frequency domain to recover W1.

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

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

             W1 contains 999 samples of the periodic function 

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

             W2 contains the analytic derivative of W1 

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

             W3 computes the frequency domain derivative and W4 displays the
             difference between the analytic and calculated. In this case,
             the difference is negligable.
            
    Remarks:
             FDERIV computes the derivative 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)

            The result is converted back into the time domain.

            To view the result in the frequency domain, try:

                spectrum(fderiv(s))

            fderiv(s, n) computes the nth derivative with:

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

            where n is a positive integer.

    See Also:
             Deriv
             FFT
             Finteg
             Integ
             Spectrum
#endif


/* frequency domain derivative */
ITERATE fderiv(s, n)
{
        local len, f, y, w, u, j;

        if (argc < 2)
        {
                if (argc < 1) error("fderiv - input series required");
                
                n = 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 derivative */
        y = fft(s) * w;

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

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

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

        y.vunits = u.vunits;

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

        return(y);
}