/*****************************************************************************
*                                                                            *
*   SOSFILT.SPL  Copyright (C) 2015 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Filters a series with 2nd order section form coefficients   *
*                                                                            *
*   Revisions:    3 Aug 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_SOSFILT

    SOSFILT  

    Purpose: 
             Filters a series with coefficients in second order section form.
                                                                        
    Format:  
             SOSFILT(sos, series, zi)

             (y, zf) = SOSFILT(sos, series, zi)

                 sos - A Nx6 array where the first 3 columns are the numerator
                       terms and the last 3 columns are the denominator terms.
                       Each row represents a 2nd order stage.

              series - A series, the input data to filter.

                  zi - Optional. A series, the initial conditions.

    Returns: 
             A series, the filtered result.

             (y, zf) = SOSFILT(sos, series, zi) returns the filtered series
             and final conditions.

    Example:
             W1: butterworth(1, 100.0, 10)
             W2: gsin(100, 1/100, 4.0) + gsin(100, 1/100, 40.0)
             W3: cas2sos(w1)
             W4: cascade(w2, w1)
             W5: sosfilt(w3, w2)
             W6: w4-w5

             W1 contains a 10 Hz lowpass Butterworth filter. W2
             contains the sum of a 4 Hz and 40 Hz sinewave. W3
             converts the cascade form coefficients of W1 into a 6x6
             array of second order sections. W4 removes the 40 Hz
             component by applying the cascade IIR filter coefficients. 
             W5 performs the same filtering with the SOS coefficients
             and W6 displays the difference.

    Example:
             W1: cas2sos(butterworth(1, 100.0, 10))
             W2: gsin(200, 1/100, 4.0) + gsin(200, 1/100, 40.0)
             W3: extract(w2, 1, 100)
             W4: extract(w2, 101, -1)
             W5: (y, zf) = sosfilt(w1, w3);y
             W6: sosfilt(w1, w4, zf)
             W7: concat(w5, w6)
             W8: sosfilt(w1, w2)
             W9: w7 - w8

             W1 contains a 10 Hz lowpass Butterworth filter in SOS form.

             W2 contains 200 samples of the sum of a 4 Hz and 40 Hz sinewave. 

             W3 and W4 split the input series into two sections.

             W5 and W6 filter each section where the final conditions
             of the first section are used as the initial conditions of
             the second section.

             W7 combines the filtered sections and W8 filters the
             original data in one step.

             W9 shows the data filtered in sections is identical to the
             unsectioned filtered data.

    Remarks:
             SOSFILT filters a series with filter coefficients in second order
             section form. The filter coefficients represent the following
             Z transform: 
  

                               -1       -2                -1       -2
                    b10 + b11 z  + b12 z       b20 + b21 z  + b22 z
         H(z) = G * ________________________ *  ________________________ ...
                               -1       -2                -1       -2
                      1 + a11 z  + a12 z         1 + a21 z  + a22 z


             where the SOS form is an Nx6 array:

             {{b10, b11, b12, 1, a11, a12},
              {b20, b21, b22, 1, a21, a22}, 

                            ...

              {bN0, bN1, bN2, 1, aN1, aN2}}

             Each row of the SOS coefficients represents a second order stage.

             Cascade coefficients are similar, but represented as a single
             column series with the coefficients in the following order:

             G, b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...


             See CASCADE to filter a series with coefficients in cascade form.

             See SOS2CAS to convert SOS form coefficients to cascade form.

    See Also:
             Cas2dir
             Cas2sos
             Cascade
             Filteq
             Sos2cas
             Sos2zp
#endif


/* filter data with an Nx6 SOS filter */
sosfilt(sos, series, zi)
{
        local c, y, zf;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("sosfilt - input SOS filter and series required");

                        sos = {};
                }

                zi = {};
        }

        if (numcols(sos) != 6)
        {
                error("sosfilt - Nx6 input SOS array required");
        }

        /* convert to cascade */
        c = sos2cas(sos);

        setdeltax(c,  deltax(series));
        setdeltax(zi, deltax(series));

        /* handles multiple columns of series */
        if (outargc > 1)
        {
                (y, zf) = cascade(series, c, zi);

                return(y, zf);
        }
        else
        {
                y = cascade(series, c, zi);

                return(y);
        }
}