/*****************************************************************************
*                                                                            *
*   NOTCHCOMB.SPL Copyright (C) 2008 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Designs a notch comb filter in direct form                 *
*                                                                            *
*   Revisions:    14 Nov 2008  RRR  Creation                                 *
*                                                                            *
*****************************************************************************/

#if @HELP_NOTCHCOMB

    NOTCHCOMB

    Purpose: Designs a notch comb filter in direct form.                                                         

    Syntax:  (b, a) = NOTCHCOMB(k, bw, Ab, fs)

                k - An integer, the notch filter order. The order is
                    equal to Fs/Fo where Fs is the sample frequency and
                    Fo is the fundamental notch filter frequency.
                    Defaults to 10.

               bw - A real, the bandwidth of each notch. The bandwidth
                    is equal to (Fo/(Fs/2))/Q where Q is the Q factor of
                    the filter. Defaults to K * 2 / 35.

               Ab - An optional real, the dB level for bandwidth BW.
                    Defaults to 10*log10(.5) or -3 dB.

               fs - Optional real, the sample rate of the filter. Defaults
                    to 1.0.
    Returns:
             The numerator and denominator coefficients of the notch filter
             in direct form.

    Example:
             Fs = 1800; Fo = 60; Q = 10;
             (b, a) = notchcomb(Fs/Fo, (Fo/(Fs/2))/10, -3, Fs);
             W1: mag(zfreq(b, a, 16*1024))

             A notch comb filter is designed with a sample rate of 1800 Hz.
             The notch frequency is 60 Hz and the Q factor is 10 for a 
             notch bandwidth of 0.00667. The coefficients are returned in 
             direct form. W1 displays the magnitude response of the filter.

             The filter removes multiples of 60 Hz line noise from the input
             data.

    Remarks:
             NOTCHCOMB returns the filter coefficients in direct form.

             The filter order K specifies the number of notch frequencies
             over the range 0 to Fs. The specific notch frequncies are
             located at F(n) = (n * Fs) / K.

             The transfer function if the notch comb filter is:

                                 1 - z^-k
                   H(z) = gain * ----------
                                 1 - a*z^-k

             where the gain is determined to normalize the filter output
             to unity gain.

             See NOTCHFILT to design and filter data with a notch comb
             filter.

  See Also:
             FILTEQ
             NOTCHFILT
             PEAKCOMB
             PEAKFILT

  References:
             [1] Sophocles J. Orfanidis
                 Introduction To Signal Processing
                 Prentice-Hall, 1996

#endif


/* design a notch comb filter in direct form */
notchiir(fs, f0, Q, Ab)
{
        local w0, BW, Gb, beta, gain, b, a;

        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1)
                                {
                                        fs = 1.0;
                                }

                                f0 = fs / 4;
                        }

                        Q = 10;
                }

                /* -3.01 dB */
                Ab = 10 * log10(0.5);
        }

        /* sanity checks */
        if (fs <= 0)
        {
                error(sprintf("%s - fs (%g) must be positive", __FUNC__, fs));
        }

        if (f0 >= fs / 2)
        {
                error(sprintf("%s - f0 (%g) must be < fs/2 (%g)", __FUNC__, f0, fs/2));
        }

        if (Q <= 0)
        {
                error(sprintf("%s - Q (%g) must be positive", __FUNC__, Q));
        }

        /* normalized fundamental */
        w0 = 2 * f0 * pi / fs;

        /* bandwidth */
        BW = w0 / Q;

        /* db attenuation at BW frequencies */
        Ab = abs(Ab);

        Gb   = 10^(-Ab / 20);
        beta = (sqrt(1 - Gb^2) / Gb) * tan(BW / 2);
        gain = 1 / (1 + beta);

        /* direct form coefficients */
        b = gain * {1,  -2 * cos(w0), 1};
        a = {1,  -2 * gain * cos(w0), (2 * gain - 1)};

        setrate(b, fs);
        setrate(a, fs);

        if (outargc > 1)
        {
                /* direct form */
                return(b, a);
        }
        else
        {
                /* cascade form */
                c = dir2cas(b, a);

                return(c);
        }
}