/*****************************************************************************
*                                                                            *
*   ADDPHASE.SPL Copyright (C) 2015 DSP Development Corporation              *
*                           All Rights Reserved                              *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Adds a constant phase value to a series                     *
*                                                                            *
*   Revisions:   28 Aug 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_ADDPHASE

    ADDPHASE

    Purpose: Adds a constant phase value to a series.

    Syntax:  ADDPHASE(s, p)

                 s - A series or array.
 
                 p - A real, the phase value in radians.

    Returns: A series or array

    Example:
             W1: gsin(1000, 1/1000, 3)
             W2: gcos(1000, 1/1000, 3)
             W3: addphase(w1, pi/2)
             W4: w2-w3

             W1 contains a 3 Hz sine wave with a sample rate of 1000 Hz.

             W2 contains a 3 Hz cosine wave with a sample rate of 1000 Hz.

             W3 adds pi/2 radians of phase shift to the sine wave.

             W4 displays the difference between the cosine and the phase
             shifted sine.

    Example:
             W1: 1..10
             W2: -w1
             W3: addphase(w1, pi)
             W4: w2-w3

             W1 contains a linear ramp.

             W2 inverts the ramp.

             W3 adds pi radians of phase shift to the sine wave.

             W4 displays the difference between the inverted ramp and
             the phase shifted ramp.

    Example:
             W1: gsqrwave(1000, 1/1000, 3)
             W2: gtriwave(1000, 1/1000, 3)
             W3: addphase(w1, pi/3)
             W4: addphase(w2, pi/3)

             W1 contains a 3 Hz square wave with a sample rate of 1000 Hz.

             W2 contains a 3 Hz triangle wave with a sample rate of 1000 Hz.

             W3 adds pi/3 radians of phase shift to the square wave
             resulting in significant distortion.

             W4 adds pi/3 radians of phase shift to the triangle wave
             resulting in significant distortion.

    Example:
             W1: gcos(1000, 1/1000, 3)
             W2: gcos(1000, 1/1000, 9)
             W3: w1 + w2
             W4: addphase(w3, pi/4)
             W5: addphase(w1, pi/4) + addphase(w2, pi/4)
             W6: w4 - w5

             W1 and W2 each contain a cosine wave.

             W3 sums the cosines.

             W4 adds pi/4 radians of phase shift to the sum.

             W5 adds pi/4 radians of phase shift to each component and
             then sums.

             W6 shows the results are identical to machine precision.

    Remarks:
             ADDPHASE adds a phase constant to the PHASE of the FFT
             of the input series. The result is transformed back to
             the time domain.

             If the series is purely real, the phase input value is
             converted to an anti-symmetric constant series and added
             to the original phase. The original DC phase value is
             preserved or set to pi if the phase input is an odd
             multiple of pi.  If the series is real and even length,
             the original Nyquist phase value is also preserved or set
             to -pi if the phase input is an odd multiple of phase.
             This handling of DC and Nyquist phase inverts the input
             series for a phase input of +- N*pi for N odd.

    See Also:
             Fft
             Magnitude
             Phase
             Phasediff
#endif


/* add p radians of phase */
ITERATE addphase(s, p)
{
        local f, n, np, p_dc, is_real;

        if (argc < 2)
        {
                if (argc < 1) s = 0;

                p = 0;
        }

        if (not(isarray(s))) error("addphase - input series required");

        /* FFT */
        f = fft(s);

        /* get length to determine even/odd processing */
        n = length(s);

        /* real flag */
        is_real = isreal(s);

        if (is_real)
        {
                /* phase at DC, pi if p odd multiple of pi */
                p_dc = ((abs(p) / pi) % pi == 1.0) ? pi : 0;

                if (n % 2 == 0)
                {
                        /* build antisymmetric even length phase values */
                        np = (n - 2) / 2;

                        /* constant phase series */
                        p  = ones(np, 1) * p;

                        /* handle DC and Nyquist */
                        p = {p_dc, p, -p_dc, -p};
                }
                else
                {
                        /* build antisymmetric odd length phase values */
                        np = (n - 1) / 2;

                        /* constant phase series */
                        p  = ones(np, 1) * p;

                        /* handle DC */
                        p = {p_dc, p, -p};
                }
        }

        /* add phase constant */
        p = phase(f) + p;

        /* recombine with magnitude */
        f = mag(f) * exp(i * p);

        /* untransform */
        s = ifft(f);

        if (is_real)
        {
                s = real(s);
        }

        return(s);
}