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

#include 

#if @HELP_FREQS

    FREQS

    Purpose: Evaluates the frequency response of a Laplace transform

    Syntax:  FREQS(b, a, N)

             (h, w) = FREQS(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:

             FREQS(b, a, w)

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

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

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

                w - a series, the angular frequencies in radians/s to 
                    evaluate the system function.


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

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

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

    Example:
             h = freqs({1}, {1, -0.5, 0.2})

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

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


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

    Example:
             freqs({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)

             freqs(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 angular frequency values range from 0.1 rad/sec to 10 rad/sec.

    Remarks:
             FREQS 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 angular frequency values, w, are in radians/s where:

                 w = 2 * pi * f

             See SFREQ to display a continuous complex frequency response
             with frequency values in Hertz.

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

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

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


/* evaluate continuous frequency response in rad/s */
freqs(b, a, w)
{
        local s, h;

        if (argc < 3)
        {
                if (argc < 2) error("freqs - b and a input coeffiecients required");
                
                w = 200;
        }

        if (length(w) == 1)
        {
                /* find important frequencies */
                w = freqsvals(b, a, w);
        }

        /* complex frequency */        
        s = i * w;

        /* evaluate system */
        h = polyval(b, s) / polyval(a, s);

        /* units */
        setvunits(w, "Radians/s");
        setvunits(h, getvunits(b));

        /* default to line style */
        setplotstyle(w, 0);
        setplotstyle(h, 0);

        if (outargc == 0)
        {
                /* plot response */
                freqsplot(h, w);
        }
        else if (outargc == 1)
        {
                /* response as XY */
                return(xy(w, h));
        }
        else 
        {
                /* separate complex response a normalized frequencies */
                return(h, w);
        }
}