/*****************************************************************************
*                                                                            *
*   HILB.SPL     Copyright (C) 2000 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Simple Hilbert Transform                                    *
*                                                                            *
*   Revisions:    8 Feb 2000  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_HILB

    HILB

    Purpose: Calculates a simple Hilbert transform of a real series

    Syntax:  HILB(s)

              s - the input series (real)


    Returns: A complex series or array

    Example:
             W1: Gnorm(1000, 0.1)
             W2: Hilb(w1)
             W3: Real(w2)
             W4: Imag(w2)

             The real part of HILB is the same as the input series.
             The imaginary part of HILB is the Hilbert Transform.

    Remarks:
             HILB uses the FFT and IFFT to calculate the Hilbert transform.

    See Also:
             Demodfm

#endif



/* simple Hilbert transformer */
hilb(s, len)
{
        local nr, nc, h, y;

        if (argc < 2)
        {
                if (argc < 1) error("hilb - input series required");

                len = -1;
        }

        (nr, nc) = size(s);

        if (len < 0)
        {
                len = nr;
        }

        y = fft(real(s), len);

        /* Hilbert weighting */
        if (len != 1)
        {
                h = {1, 2 * ones(fix((len - 1) / 2), 1), ones(1 - (len % 2), 1), zeros(fix((len - 1) / 2), 1)};

                y = y * h[.., ones(1, nc)];
        }
        
        y = ifft(y);

        /* repair deltax */
        setdeltax(y, deltax(s));

        return(y);
}