apskdem3.c

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

/*
 * APSKDEMO  ASK/PSK Demodulation phase detection and decision part
 *
 *           Syntax:  [sys, x0] = apskdemo(t,x,u,flag,M,Fc,N,P,R)
 *                      where M is the dimension of N and P.
 *                            N is the distribution of the centric circle.
 *                            P is the initial phase for each circle.
 *                            R is the assigned amplitude information.
 *                      there are three input:
 *                            u[0] is the modulated signal, block input.
 *                            u[1] is the tirgger signal which means the
 *                                  begining of the next input.
 * Wes Wang 5/2/94
 * Copyright (c) 1994-96 The MathWorks, Inc.
 * All Rights Reserved
 * $Revision: 1.1 $  $Date: 1996/04/01 19:01:34 $
 */

/* specify the name of this S-Function. */
#define S_FUNCTION_NAME apskdem3

/* Defines for easy access  the matrices which are passed in */
#define PI          3.14159265358979
#define NUM_ARGS    5
#define NUM_LAYER   ssGetArg(S, 0)
#define FRE_CARR  ssGetArg(S, 1)
#define NUM_IN_L  ssGetArg(S, 2)
#define PHA_IN_L  ssGetArg(S, 3)
#define AMP_IN_L  ssGetArg(S, 4)

/* include simstruc.h for the definition of the SimStruct and macro definitions. */
#include <math.h>
#ifdef MATLAB_MEX_FILE
#include <stdio.h>
#endif
#include "simstruc.h"

#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif

/*
 * mdlInitializeSizes - initialize the sizes array
 */

static void mdlInitializeSizes (S)
    SimStruct *S;
{
    if (ssGetNumArgs(S) == NUM_ARGS) {
      int    i, j, k, mulNum, idphNum, numInL;
        int    CurrentBase, numLayer;
           int    test;
            int    phaseIndex[256];
         double tmp, phases[256];   /* current maximum is 256 point */
           numLayer = mxGetPr(NUM_LAYER)[0];
       
#ifdef MATLAB_MEX_FILE
          if ((mxGetN(NUM_LAYER) != 1) || (mxGetM(NUM_LAYER) != 1)) {
         mexErrMsgTxt("The number of layer must be a scalar");
       }
       if ((mxGetN(FRE_CARR) != 1) || (mxGetM(FRE_CARR) != 1)) {
           mexErrMsgTxt("The carrier frequency must be a scalar");
     }
       if ((mxGetN(NUM_IN_L) * mxGetM(NUM_IN_L) != numLayer)) {
            mexErrMsgTxt("The vector size for number in each layer is not consistant");
         }
       if ((mxGetN(PHA_IN_L) * mxGetM(PHA_IN_L) != numLayer)) {
            mexErrMsgTxt("The vector size for initial phase in each layer is not consistant");
          }               
#endif

        /* block multiple number should be 2, 4, 8, 16, 32, 64, 128, or 256 */
        mulNum = 0;
        for (i = 0; i < numLayer; i++)
            mulNum += mxGetPr(NUM_IN_L)[i];

        /* number of independant phases */
        idphNum = -1;
        CurrentBase = 0;

        /* phaseIndex counted from the very out circle */
        for (i = numLayer-1; i >= 0; i--){
            if (i != numLayer - 1)
                CurrentBase += (int)mxGetPr(NUM_IN_L)[i+1];
            numInL = (int)mxGetPr(NUM_IN_L)[i];
            for (j = 0; j < numInL; j++) {
                tmp = mxGetPr(PHA_IN_L)[i] + ((double)(numInL - j - 1))*PI*2/numInL;
                k = 0; test = 1;

                while ((k <= idphNum) && test) {
                    if (fabs(tmp - phases[k]) <=0.00001) {
                        phaseIndex[CurrentBase + j] = k;
                        test = 0;
                    }
                    k++;
                }
                if (test) {
                    idphNum++;
                    phases[idphNum] = tmp;
                    phaseIndex[CurrentBase + j] = idphNum;
                }
            }
        }
        idphNum++;
        
        ssSetNumContStates(    S, idphNum); /* number of continuous states */
        ssSetNumDiscStates(    S, 0);       /* number of discrete states */
        ssSetNumInputs    (        S, 2);   /* number of inputs */
        ssSetNumOutputs   (        S, 1);   /* number of outputs */
        ssSetDirectFeedThrough(S, 0);       /* direct feedthrough flag */
        ssSetNumSampleTimes(   S, 1);       /* number of sample times */
        ssSetNumInputArgs(     S, NUM_ARGS);/* number of input arguments */
        /* last time of u[1]==1;
         * last calculated output;
         * independent number of state phase shift
         */
        ssSetNumRWork(         S, idphNum + 2);       /* number of real work vector elements */
        /* (1) state number
         * (2) real work vector
         * (3) point vector/ same as the multi number
         */
        ssSetNumIWork(         S, mulNum + 2);      /* number of integer work vector elements */
        ssSetNumPWork(         S, 0);      /* number of pointer work vector elements */
        } else {
#ifdef MATLAB_MEX_FILE
       char err_msg[256];
      sprintf(err_msg, "Wrong number of input arguments passed to S-function MEX-file.\n"
             "%d input arguments were passed in when expecting %d input arguments.\n", ssGetNumArgs(S) + 4, NUM_ARGS + 4);
           mexErrMsgTxt(err_msg);
#endif
    }
}

/*
 * mdlInitializeSampleTimes - initialize the sample times array
 *
 * This function is used to specify the sample time(s) for your S-function.
 * If your S-function is continuous, you must specify a sample time of 0.0.
 * Sample times must be registered in ascending order.
 */

static void mdlInitializeSampleTimes(S)
    SimStruct *S;
{
    ssSetSampleTimeEvent(S, 0, 0.0);
    ssSetOffsetTimeEvent(S, 0, 0.0);
}

/*
 * mdlInitializeConditions - initialize the states
 * Initialize the states, Integers and real-numbers
 */

static void mdlInitializeConditions(double *x0, SimStruct *S)
{
        double *resetTime   = ssGetRWork(S);
        double *lastOut   = ssGetRWork(S) + 1;
        double *iniPhases = ssGetRWork(S) + 2;
        int    *staNum    = ssGetIWork(S);
        int    *mulNum    = ssGetIWork(S) + 1;
        int    *phaseInd  = ssGetIWork(S) + 2;

            int    i, j, k, mlNum, idphNum, test, CurrentBase, numLayer,numInL;
             double tmp, phases[256];   /* current maximum is 256 point */
           int    phaseIndex[256];

        /* multiple number */
            numLayer = mxGetPr(NUM_LAYER)[0];
        mlNum = 0;
        for (i = 0; i < numLayer; i++) {
            mlNum += (int)mxGetPr(NUM_IN_L)[i];
        }

        /* number of independant phases */
        idphNum = -1;
        CurrentBase = 0;
        /* phaseIndex counted from the very out circle */
        for (i = numLayer-1; i >= 0; i--){
            if (i != numLayer - 1)
                CurrentBase += (int)mxGetPr(NUM_IN_L)[i+1];
            numInL = (int)mxGetPr(NUM_IN_L)[i];
            for (j = 0; j < numInL; j++) {
                tmp = mxGetPr(PHA_IN_L)[i] + ((double)(numInL - j - 1))*PI*2/numInL;
                k = 0; test = 1;
                while ((k <= idphNum) && test) {
                    if (fabs(tmp - phases[k]) <=0.00001) {
                        phaseIndex[CurrentBase + j] = k;
                        test = 0;
                    }
                    k++;
                }
                if (test) {
                    idphNum++;
                    phases[idphNum] = tmp;
                    phaseIndex[CurrentBase + j] = idphNum;
                }
            }
        }
        idphNum++;
        
        *resetTime = 0;
        *lastOut = 0;
        for (i = 0; i < idphNum; i++) {
               *x0++ = 0.;
             *iniPhases++ = phases[i];
        }
        *staNum = idphNum;
        *mulNum = mlNum;
        for (i = 0; i < mlNum; i++)
            *phaseInd++ = phaseIndex[i];
}

/*
 * mdlOutputs - compute the outputs
 *
 * In this function, you compute the outputs of your S-function
 * block.  The outputs are placed in the y variable.
 */

static void mdlOutputs(y, x, u, S, tid)
    double *y, *x, *u;
    SimStruct *S;
    int tid;
{
    double *resetTime = ssGetRWork(S);
    double *lastOut   = ssGetRWork(S) + 1;
    double tmpTime;
    
    /* Time cumulated so far */
    tmpTime = ssGetT(S) - *resetTime;

    /* In the case of the output to be renewed */
    if ((u[1] != 0) && (tmpTime > 0)) {
        /* major computation */
        /* x[0] to x[staNum - 1] is the correlation result */

        int    *staNum    = ssGetIWork(S) ;        
        int    *mulNum    = ssGetIWork(S) + 1;
        int    *phaseInd  = ssGetIWork(S) + 2;
        double *iniPhase   = ssGetRWork(S) + 2;
        double mini, tmp, ampInL, submini, submaxi;
        int    i, j, decis, numBase, numInL, numLayer, decLay, subdecis, subsup;
        int    phaseIndex;              
        int    yesprint = 0;

        numLayer = mxGetPr(NUM_LAYER)[0];
        /* first find the best amplitude + phase match */
        numBase = 0;
        numInL = 0;
        for (i= numLayer-1; i>=0; i--) {
            numBase += numInL;
            numInL = (int)mxGetPr(NUM_IN_L)[i];
            ampInL = mxGetPr(AMP_IN_L)[i];
            submini = 10000000000.0;
            subdecis = numBase;
            submaxi = 0;
            subsup = numBase + 1;
            for (j=0; j<numInL; j++) {
                phaseIndex = *phaseInd;
                *phaseInd++;
#ifdef MATLAB_MEX_FILE
                if (yesprint) {
                        char   err_msg[256];
                    sprintf(err_msg, "phaseInd %d, staNum %d, i %d, j %d, numInL %d\n",
                     phaseIndex, *staNum, i, j, numInL);
                    mexPrintf(err_msg);
                }
#endif
                /* find the minimum */
                if (ampInL > 0)
                    tmp = fabs(x[phaseIndex] - ampInL * tmpTime / 2) / ampInL;
                else
                    tmp = fabs(x[phaseIndex] - ampInL * tmpTime / 2);                
                if (submini >= tmp) {
                    subdecis = *mulNum - numBase - j;
                    submini = tmp;
                }
                /* find the maximum */
                if (submaxi <= x[phaseIndex]) {
                    subsup = *mulNum - numBase - j;
                    submaxi = x[phaseIndex];
                }
#ifdef MATLAB_MEX_FILE
                if (yesprint) {
                        char   err_msg[256];                
                    sprintf(err_msg, "decision %d, mini %8.5f, tmp %8.5f, x[phaseInd] %8.5f, phaseInd %d, order %d, phase %6.3f\n",
                        decis, mini, tmp, x[phaseIndex], phaseIndex, numBase+j,iniPhase[phaseIndex]);
                    mexPrintf(err_msg);
                    sprintf(err_msg, "subdecis %d, submini %8.5f, subsup %d, submaxi %8.5f\n",
                        subdecis, submini, subsup, submaxi);
                    mexPrintf(err_msg);                    
                }
#endif
            }
            if (i == numLayer - 1) {
                mini = submini;
                decis = subdecis;
            } else if ((mini >= submini) && (subdecis == subsup)) {
                mini = submini;
                decis = subdecis;
            }
        }
        *lastOut = decis - 1;
        *resetTime = ssGetT(S);        
    }
    y[0] = *lastOut;
}


/*
 * mdlUpdate - perform action at major integration time step
 *
 * This function is called once for every major integration time step.
 * Discrete states are typically updated here, but this function is useful
 * for performing any tasks that should only take place once per integration
 * step.
 */

static void mdlUpdate(x, u, S, tid)
    double *x, *u;
    SimStruct *S;
    int tid;
{
    double *resetTime   = ssGetRWork(S);
    double CTime     = ssGetT(S);
    
    /* if input 3 is a nonzero value, reset the state of the integrator */
    if (*resetTime >= CTime) {
        int i;
        int    *staNum    = ssGetIWork(S);
        for (i = 0; i < *staNum; i++)
        x[i] = 0.0;
    }
}

/*
 * mdlDerivatives - compute the derivatives
 *
 * In this function, you compute the S-function block's derivatives.
 * The derivatives are placed in the dx variable.
 */

static void mdlDerivatives(dx, x, u, S, tid)
    double *dx, *x, *u;
    SimStruct *S;
    int tid;
{
    int    i;
    int    *staNum    = ssGetIWork(S);
    double *resetTime   = ssGetRWork(S);
    double *iniPhases = ssGetRWork(S) + 2;
    double CTime     = ssGetT(S);
    double fct;
    
    fct = (mxGetPr(FRE_CARR)[0]) * (CTime - *resetTime) * 2 * PI;
    for (i = 0; i < *staNum; i++) {
        dx[i] = sin(fct + iniPhases[i]) * u[0];
      }
}

/*
 * mdlTerminate - called when the simulation is terminated.
 *
 * In this function, you should perform any actions that are necessary
 * at the termination of a simulation.  For example, if memory was allocated
 * in mdlInitializeConditions, this is the place to free it.
 */

static void mdlTerminate(S)
    SimStruct *S;
{
}

#ifdef  MATLAB_MEX_FILE    /* Is this file being compiled as a MEX-file? */
#include "simulink.c"      /* MEX-file interface mechanism */
#else
#include "cg_sfun.h"       /* Code generation registration function */
#endif