/*****************************************************************************
*                                                                            *
*   SFREQ.SPL    Copyright (C) 2009 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Evaluates the frequency response of an s transform          *
*                                                                            *
*   Revisions:   17 Aug 2009  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_SFREQ

    SFREQ

    Purpose: Evaluates the frequency response of a Laplace transform in Hertz

    Syntax:  SFREQ(b, a, N)

             (h, f) = SFREQ(b, a, N)

                b - a series, numerator (i.e. zero) coefficients

                a - a series, denominator (i.e. pole) coefficients,

                N - an optional integer, number of output samples,
                    defaults to 200.


    Alternative Syntax:

             SFREQ(b, a, f)

             (h, w) = SFREQ(b, a, w)

                b - a series, numerator (i.e. zero) coefficients

                a - a series, denominator (i.e. pole) coefficients,

                f - a series, the frequencies in Hertz to evaluate
                    the system function.


    Returns: Displays the magnitude and phase response in two windows.

             h = SFREQ(b, a) returns the complex frequency response
             as one XY series.

             (h, f) = FREQA(b, a) returns the complex frequency response
             and frequency values in Hertz as two separate series.

    Example:

             h = sfreq({1}, {1, -0.5, 0.2})

             h contains 200 samples of the frequency response of the
             s transform:

                           1
             H(s) = --------------
                     2               
                    s - 0.5s + 0.2


             The normalized frequency values range from 0 to 10 radians/s.

    Example:

             sfreq({1}, {0.5, -0.2}, 1024)

             Returns 1024 samples of the magnitude and phase response
             of the system in two separate windows.

    Example:

             b = {0.2, 0.3, 1.0}
             a = {1.0, 0.4, 1.0}
             w = logspace(-1, 1)

             sfreq(b, a, w)            

             Displays 100 samples of the magnitude and phase response of 
             the system:


                        2
                    0.2s  + 0.3s + 1
             H(s) = ----------------
                       2
                      s  + 0.4s + 1


             The frequency values range from 0.1 rad/sec to 10 rad/sec.


    Remarks:
             SFREQ displays the magnitude and phase response of the
             continuous system specified by the Laplace transform:

                     B(s)   b[1]*s^n + b[2]*s^(n-1) + ... + b[n+1]
             H(s) =  ---- = --------------------------------------
                     A(s)   a[1]*s^m + a[2]*s^(m-1) + ... + a[m+1]


             If no output arguments are supplied, the magnitude and phase
             response are displayed in two Windows.

             The frequency values, f, are in Hertz where:

                 f = w / (2 * pi)

             See FREQS to display a continuous complex frequency response
             with angular frequency values in radians/sec.

             SFREQ(b, a) or SFREQ(b, a, N) chooses frequencies to best
             capture the magnitude characteristics of the system.

             SFREQ(b, a, w) where w is a series, computes the frequency 
             response at each frequency sample of w.

    See Also:
             Filteq
             Freqs
             Freqz
             Invfreqs
             Invfreqz
             Mag
             Phase
             Residues
             Residuez
             Zfreq
#endif


/* evaluate continuous frequency response in Hertz */
sfreq(b, a, f)
{
        local w;

        w = (argc < 3) ? 200 : f;

        if (isarray(w))
        {
                /* convert to radian frequency */
                w *= 2 * pi;
        }

        (h, w) = freqs(b, a, w);

        /* convert to Hertz */
        f = w / (2 * pi);
        setvunits(f, "Hertz");

        if (outargc == 0)
        {
                freqsplot(h, f);
        }
        else if (outargc > 1) 
        {
                return(h, f);
        }
        else 
        {
                return(xy(f, h));
        }
}