/*****************************************************************************
*                                                                            *
*   FIRSAMP.SPL   Copyright (C) 2000, 2017 DSP Development Corporation       *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     FIR filter design using frequency sampling                 *
*                                                                            *
*   Revisions:    14 Mar 2000  RRR  Creation                                 *
*                 16 Nov 2017  RRR  interpolation factor                     *
*                                                                            *
*****************************************************************************/


#if @HELP_FIRSAMP

    FIRSAMP

    Purpose: Designs an arbitrary FIR filter using frequency sampling

    Syntax:  FIRSAMP(f, m, fac)

                f - An XY series specifying the desired frequencies
                    (in Hertz) and magnitudes of the filter OR an
                    explicit series specifying the frequencies only

                m - An optional series, explicit desired magnitudes

              fac - Optional. An integer, the interpolation factor used
                    to design the filter. Defaults to 2. 

    Returns: A series

    Example:
             W1: {0, 1, 10, 20, 30}
             W2: {1, 1, 2,  0, 0}
             W3: xy(w1, w2)
             W4: firsamp(w3)
             W5: mag(fft(w4, bestpow2(length(w4))))
             W6: extract(w5, 1, int(length(w5)/2));overp(w3, lred)

             W4 contains a 121 point linear phase FIR filter. The filter
             has unity gain from at 0 and 1 Hz and a gain of 2 at
             10 Hz. W5 shows the magnitude response of the filter. The
             original specification is overplotted onto the actual
             response in W6.

    Example:
             W1: {0, 1, 10, 20, 30}
             W2: {1, 1, 2,  0, 0}
             W3: xy(w1, w2)
             W4: firsamp(w3, 10)
             W5: mag(fft(w4, bestpow2(length(w4))))
             W6: extract(w5, 1, int(length(w5)/2));overp(w3, lred)

             Same as above except the interpolation factor is increased
             to 10. W4 contains a 599 point linear phase FIR filter.

    Example:
             W7: firsamp(w1, w2, 10)

             Same as above except the frequencies and magnitudes
             are specified explicitly.

    Remarks:
             FIRSAMP sorts the input frequencies in ascending order.
             If a frequency of 0 Hz is not specified, a 0 Hz term
             equal to the magnitude of the frequency nearest 0
             is added to the list.

             FIRSAMP produces non-causal filters with linear phase.

             The input filter frequency response specification is
             interpolated to a regular interval and the IFFT is applied
             to derive the filter coefficients.  Increasing FAC
             produces a filter with more coefficients but the frequency
             response will match the input specification more closely.

             To reduce ringing at sharp transitions, multiply the
             resulting filter coefficients with a windowing function.
             For example:

             W1: xy({0, 0.3, 0.4, 0.6, 0.7, 0.9}, {0, 1, 0, 0, 0.5, 0.5})
             W2: firsamp(w1, 5)
             W3: hamming(w2)
             W4: filtmag(w2, {1}, 32*1024)
             W5: filtmag(w3, {1}, 32*1024)

             Note: FILTMAG is available in DADiSP/Filters.

    See Also:
             Fft
             Fftshift
             Xyinterp
#endif


/* FIR design via frequency sampling */
firsamp(f, m, fac)
{
        local s, h, dx;

        /* inputs */
        (f, m, fac) = firsamp_parse_args(f, m, fac);

        /* sort values in ascending order */
        h = grade(f, 1);
        f = f[h];
        m = m[h];

        /* get XY plot so we can interpolate it */
        m = xy(f, m);

        /* add DC term if required */
        if (f[1] != 0)
        {
                m = {xy({0}, {m[1]}), m};
                printf("firsamp: Setting DC Value to %g", m[1]);
        }

        /* get automatic interpolation interval - one point series  */
        dx = xyinterp(m, -1, 1);

        /* interpolate to equally spaced samples */
        m = xyinterp(m, castreal(dx / fac));

        /* form full magnitude response */
        m = {m, (rev(extract(m, 2, -1)))};

        /* inverse transform to get impulse response */
        h = real(ifft(m));

        /* shift to non causal - automatically sets linear phase */
        h = fftshift(h);

        return(h);
}


/* input args */
firsamp_parse_args(f, m, fac)
{
        if (argc < 3)
        {
                if (argc < 2)
                {
                        m = yvals(f);
                        f = xvals(f);

                        fac = -1;
                }
                else
                {
                        if (isscalar(m))
                        {
                                fac = m;
                                m   = yvals(f);
                                f   = xvals(f);
                        }
                        else
                        {
                                fac = -1;
                        }
                }
        }

        if (fac <= 0)
        {
                fac = 2;
        }

        return(f, m, fac);
}