gfdeconv.c

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

/*============================================================================
 *      Syntax:     [q, r] = gfdeconv(b, a, p)
 *GFDECONV GF(P) polynomial deconvolution or GF(P^M) elements dividing.
 *       [Q, R] = GFDECONV(B, A) computes the quotient Q and remainder R of
 *       deconvolution B by A in GF(2).
 *
 *       [Q, R] = BFDECONV(B, A, P) computes the quotient and remainder R of
 *       deconvolution B by A in GF(P) when P is a scalar prime number.
 *       When P is a matrix containing the tuple of all elements in GF(p^M),
 *       this function takes A and B as indices (power number of the
 *       exponential form) of GF(p^M) elements. The output Q is
 *       alpha^C = alpha^A / alpha^B in GF(p^M). The computation is
 *       element-by-element computation. You can generate the tuple of all
 *       elements in GF(Q^M) by P = GFTUPLE([-1:Q^M-2]', M, Q). 
 *   
 *       In polynomial computation, A, B, and C are in ascending order, i.e.,
 *       A = [a_0, a_1, a_2,..., a_(n-1), a_n] represents
 *       A(X) = a_0 + a_1 X + a_2 X^2 +...+ a_(n-1) X^(n-1) + a_n X^n
 *       a_i must be a element in GF(P).
 *
 *       In power representation form, [-Inf, 0, 1, 2, ...] represents 
 *       [0, 1, alpha, alpha^2, ...] in GF(p^m).
 *
 *      See also GFADD, GFMUL, GFTUPLE, GFPRETTY
 *
 *============================================================================
 *     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:13 $
 *==========================================================================
 */
#include <math.h>
#include "mex.h"
#include "gflib.c"
void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
    int     mb, nb, ma, na, np, mp, len_b, len_a, len_p, len_q, len_r, i;
    int     *pbb, *paa, *pp, *pqq, *prr, *Iwork1, *Iwork2;
    double  *pb, *pa, *p, *pq, *pr;

    if ( nrhs < 2 ){
        mexErrMsgTxt("Not enough input for GFDECONV!");
    }else if ( nrhs == 2 ){
        len_p = 1;
    }else if ( nrhs > 2 ){
        p = mxGetPr(prhs[2]);
        np= mxGetM(prhs[2]);
        mp= mxGetN(prhs[2]);
        len_p = np*mp;
    }
    /* get input arguments */
    pb = mxGetPr(prhs[0]);
    pa = mxGetPr(prhs[1]);
    mb = mxGetM(prhs[0]);
    nb = mxGetN(prhs[0]);
    ma = mxGetM(prhs[1]);
    na = mxGetN(prhs[1]);
    len_b = mb*nb;
    len_a = ma*na;

    /* variable type conversion for calling functions in gflib.c */
    pbb = (int *)mxCalloc(len_b, sizeof(int));
    paa = (int *)mxCalloc(len_a, sizeof(int));
    pp  = (int *)mxCalloc(len_p, 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];
    if ( nrhs == 2 ){
        pp[0] = 2;
    }else{
        for (i=0; i < len_p; i++)
            pp[i] = (int) p[i];
    }

	/* truncate input */
    Iwork1 = (int *)mxCalloc(len_b+len_a, sizeof(int));
    gftrunc(pbb, &len_b, len_p,Iwork1);
    gftrunc(paa, &len_a, len_p,Iwork1+len_b); 
	
    if(len_a > len_b)
        len_q = 1;
    else
        len_q = len_b - len_a + 1;
    len_r = len_b;
    pqq = (int *)mxCalloc(len_q, sizeof(int));
    for(i=0; i<len_q; i++)
        pqq[i] = 0;
    prr = (int *)mxCalloc(len_r, sizeof(int));
    for(i=0; i<len_r; i++)
        prr[i] = 0;
    /* case of polynomial calculation */

    if (len_p <= 1){
        /* input check up */    
        for (i=0; i < len_b; i++){
            if (pb[i] < 0 || pb[i] != floor(pb[i]) || pb[i] >= pp[0])
                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] >= pp[0] )
                mexErrMsgTxt("The polynomial coeficients must be in GF(P)");
        }
        /* call gfdeconv() in gflib.c */
        Iwork2 = (int *)mxCalloc(14*len_b-2*len_a+3, sizeof(int));
        gfdeconv(pbb, len_b, paa, len_a, *pp, pqq, len_q, prr, &len_r, Iwork2);
	} else {
        /* call gfpdeconv() in gflib.c */
        Iwork2 = (int *)mxCalloc(len_a-1+(2+mp)*np, sizeof(int));
        gfpdeconv(pbb, len_b, paa, len_a, pp, np, mp, pqq, prr, &len_r, Iwork2);
	}
    pq = mxGetPr(plhs[0]=mxCreateFull(1, len_q, 0));
    for(i=0; i < len_q; i++){
        if(pqq[i] == -Inf)
            pq[i] = -mexGetInf();
        else
            pq[i] = (double)pqq[i];
	}
    pr = mxGetPr(plhs[1]=mxCreateFull(1, len_r, 0));
    for(i=0; i < len_r; i++){
        if(prr[i] < 0)
            pr[i] = -mexGetInf();
        else
            pr[i] = (double)prr[i];
    }
	return;
}
/* end of gfdeconv.c */