gffilter.c

数字通信第四版原书的例程

/*============================================================================
 *      Syntax:     y = gffilter(b, a, x, p)
 * GFFILTER GF(P) polynomial filter computation.
 *   Y = FILTER(B, A, X) filters the data in vector X with the
 *   filter described by vectors A and B to create the filtered
 *   data Y in GF(2).
 *
 *    Y = FILTER(B, A, X, P) filters the data in GF(P).
 *
 *   The filter uses the standard difference equation:
 *   y(n) = b_nb*x(n) + b_(nb-1)*x(n-1) + ... + b_0*x(n-nb)
 *              - a_(na-1)*y(n-1) - ... - a_0*y(n-na)
 *
 *   Different from all of the other GF functions in this toolbox, this
 *   this function uses descending ordered polynomial, i.e.
 *     A = [a_n, a_(n-1), ..., a_2, a_1, a_0] represents
 *         A(X) = a_n X^n + ... + a_2 X^2 + a_1 X + a_0
 *    a_i must be a non-negative integer less than P.
 *
 *       See also GFADD, GFMUL, GFDIV, GFTUPLE
 *
 * PS: One may directly use the built-in function filter.
 *     However, the computation may not be accurate if
 *     the computation order is high.
 *============================================================================
 *     Original designed by Wes Wang,
 *     Jun Wu,     The Mathworks, Inc.
 *     Dec-12, 1995
 *
 *     Copyright (c) 1995-96 by The MAthWorks, Inc.
 *     All Rights Reserved
 *     $Revision: 1.1 $  $Date: 1996/04/01 18:14:19 $
 *==========================================================================
 */
#include <math.h>
#include "mex.h"
#include "gflib.c"
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
    int     mb, nb, ma, na, mx, nx, len_b, len_a, len_x, p, i;
    int     *pxx, *pyy, *pbb, *paa, *Iwork;
    double  *pb, *pa, *px, *py;

    if ( nrhs < 3 ){
        mexErrMsgTxt("Not enough input for GFFILTER!");
    } else if ( nrhs == 3 ){
        p = 2;
   	} else if ( nrhs == 4 ){
        p = (int)mxGetScalar(prhs[3]);
    }
    pb = mxGetPr(prhs[0]);
    pa = mxGetPr(prhs[1]);
    px = mxGetPr(prhs[2]);    
    mb = mxGetM(prhs[0]);
    nb = mxGetN(prhs[0]);
    ma = mxGetM(prhs[1]);
    na = mxGetN(prhs[1]);
    mx = mxGetM(prhs[2]);
    nx = mxGetN(prhs[2]);

    /* variable assignment. */
    len_b = mb*nb;
    len_a = ma*na;   	
    len_x = mx*nx;

    /* input check up */    
    for (i=0; i < len_b; i++){
        if (pb[i] < 0 || pb[i] != floor(pb[i]) || pb[i] >= p)
            mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
	}
    for (i=0; i < len_a; i++){
        if (pa[i] < 0 || pa[i] != floor(pa[i]) || pa[i] >= p )
            mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
    }
    for (i=0; i < len_x; i++){
        if (px[i] < 0 || px[i] != floor(px[i]) || px[i] >= p)
            mexErrMsgTxt("The input in Galois space must be in GF(P)");
	}

    /* variable type conversion for calling functions in gflib.c */
    pbb = (int *)mxCalloc(len_b, sizeof(int));
    paa = (int *)mxCalloc(len_a, sizeof(int));
    pxx = (int *)mxCalloc(len_x, sizeof(int));
    for (i=0; i < len_b; i++)
        pbb[i] = (int) pb[i];        
    for (i=0; i < len_a; i++)
        paa[i] = (int) pa[i];
    for (i=0; i < len_x; i++)
        pxx[i] = (int) px[i];

	/* truncate input */
    Iwork = (int *)mxCalloc(2*(len_b+len_a+len_x-1), sizeof(int));
    gftrunc(pbb, &len_b, 1,Iwork); 
    gftrunc(paa, &len_a, 1,Iwork+len_b); 

    /* assign (int *)pyy for output of gffilter() */
    pyy = (int *)mxCalloc(len_x, sizeof(int));
    /* call function gffilter() in gflib.c */
    gffilter(pbb, len_b, paa, len_a, pxx, len_x, p, pyy, Iwork+len_b+len_a);

    py = mxGetPr(plhs[0] = mxCreateFull(1, len_x, 0));
    for (i=0; i < len_x; i++)
        py[i] = (double)pyy[i];
    return;
}
/*--end of gffilter.c--*/