/*****************************************************************************
*                                                                            *
*   FMDEMOD.SPL  Copyright (C) 2000-2013 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    FM Demodulation using the Hilbert Transform                 *
*                                                                            *
*   Revisions:    8 Feb 2000  RRR  Creation                                  *
*                22 May 2013  RRR  optional fmin/fmax modulation range       *
*                                                                            *
*****************************************************************************/

#if @HELP_FMDEMOD

    FMDEMOD

    Purpose: Demodulates an FM series using the Hilbert Transform

    Syntax:  FMDEMOD(s, fc, fdev)

                 s - the input series

                fc - Optional, a real. The center frequency of the frequency
                     modulated result. Defaults to -1, estimate the center
                     frequency using the FFT.

              fdev - Optional, a real. The frqeuency deviation. Defaults to 1.

    Returns: A series or array

    Example:
             W1: gtriwave(1000,.001, 4)
             W2: W1 * 100 + 20
             W3: cos(integ(2*pi*w2))
             W4: demodfm(w3)

             W2 represents the scaled information signal and W3 is the
             resulting frequency modulated signal. The amplitude of W2
             determines the instantaneous frequency of W3. The instantaneous
             frequency of W3 ranges from:

                             min(w3) == 20 Hz

                                     to

                             max(w3) == 120 Hz


             W4 is the demodulated result.

    Example:
             W1: gtriwave(1000, .001, 4) * 100 + 20
             W2: fmmod(w1, 200, 2)
             W3: fmdemodfm(w2, 200, 2)
             
             Same as above except the minimum and maximum modulation
             frequencies are specified in W4 to scale the result.

    Remarks:
             FMDEMOD uses HILB to calculate the Hilbert Transform.

             If FC is < 0, the center frequency is estimated using
             SINFIT3.

             See FMMOD to frequency modulate a series.

             See DEMODAM to demodulate an amplitude modulated series.

    See Also:
             Demodam
             Hilb
             Modfm

#endif


/* demodulate an FM series using the Hilbert transform */
ITERATE fmdemod(s, fc, fdev)
{
        local h, fm, fit, coef;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("fmdemod - input series required");
                        
                        /* estimate frequency */
                        fc = -1;
                }

                fdev = 1;
        }

        if (fc < 0)
        {
                /* estimate center frequency */
                (fit, coef) = sinfit3(s);

                fc = coef[2];
        }

        /* hilbert transform */
        h  = hilb(s) * exp(-i * 2 * pi * fc * xvals(s));
        fm = deriv(unwrap(phase(h))) / (2 * pi * fdev);

        /* restore the offset */
        setxoffset(fm, xoffset(s));

        return(fm);
}