/*****************************************************************************
*                                                                            *
*   ZBSFILT.SPL Copyright (C) 2013 DSP Development Corporation               *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Band stop filtering by FFT zeroing                          *
*                                                                            *
*   Revisions:    2 Sep 2013  RRR  Creation, from zlpfilt.spl                *
*                                                                            *
*****************************************************************************/


#if @HELP_ZBSFILT

    ZBSFILT

    Purpose: Performs band stop filtering by zeroing the FFT values.

    Syntax:  ZBSFILT(s, fs1, fs2, mask)

                 s - A series, the input series.

               fs1 - A real, the first band stop cutoff frequency. Defaults
                     to rate(s) / 10.

               fs2 - A real, the second band stop cutoff frequency. Defaults
                     to rate(s) / 2.

              mask - Optional. An integer, return the frequency domain mask
                     only. Defaults to 0, return processed input.
 
    Returns:
             A series, the band stop filtered output.

    Example:
             W1: gsin(1000,0.001,10)+gsin(1000,0.001,20)+gsin(1000,0.001,30)+gsin(1000,0.001,40)
             W2: zbsfilt(w1, 15, 35)

             W1 contains the sum of 10 Hz, 20 Hz, 30 Hz and 40 Hz sinewaves.

             W2 filters out the 20 Hz and 30 Hz components leaving a series
             consisting of 10 Hz and 40 Hz sine components.

    Remarks:
             ZBSFILT performs band stop filtering by zeroing out the
             stop band frequency components of the FFT of the input.

             FFT zeroing is a crude but fast filtering operation.
             However, the pass band and transition edge of the
             equivalent filter will exhibit significant ripple. For a
             flatter response, use traditional filter design methods
             such as those offered by DADiSP/Filters.

             If fs1 < 0 or fs2 < 0, no filtering is performed.

             If mask == 1, the frequency domain binary mask used to zero
             out the stopband frequencies is returned. 

             See ZBPFILT for a band pass filter.

  See Also:
             Bandstop
             DADiSP/Filters
             Zbpfilt
             Zhpfilt
             Zlpfilt
#endif


/* FFT zeroing band stop filter */
ITERATE zbsfilt(s, fs1, fs2, mask)
{
        local len, bw, f, q;

        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1) error("zbsfilt - input series required");

                                fs1 = rate(s) / 10;
                        }

                        fs2 = rate(s) / 2;
                }

                /* default to process input */
                mask = 0;
        }

        len = length(s);

        /* bin frequency width */
        bw  = rate(s) / len;

        /* build band pass frequency mask */
        f = ones(len, 1, bw);

        if (fs1 < 0 || fs2 < 0)
        {
                /* no filtering, mask is inverted here */
                return(mask ? not(f) : s);
        }

        /* index of cutoff frequencies */
        q1 = int(fs1 / bw) + 1;
        q2 = int(fs2 / bw) + 1;

        if (q1 > len / 2)
        {
                q1 = len;
        }

        if (q2 > len / 2)
        {
                q2 = len;
        }

        /* DC range passband */
        f[1..q2] = 0;

        /* DC range stopband */
        f[1..q1] = 1;

        /* nyquist range pass frequency */
        if (q2 >= 2 && q2 < len)
        {
                f[(len - q2 + 2)..len] = 0;
        }

        /* nyquist range stop frequency */
        if (q1 >= 2 && q1 < len)
        {
                f[(len - q1 + 2)..len] = 1;
        }

        if (mask)
        {
                /* low bandpass mask */
                return(f);
        }
        else
        {
                /* frequency domain processing */
                f = ifft(fft(s) * f);

                if (isreal(s))
                {
                        f = real(f);
                }

                setxoffset(f, xoffset(s));

                return(f);
        }
}