apskdemo.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 demodulated amplitude.
 *                            u[2] 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:46 $
 */

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

/* 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, 3);   /* 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[2]==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(x0, S)
     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];

        char err_msg[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[2] != 0) && (tmpTime > 0)) {
        /* major computation */
        /* u[1] is the estimated amplitude */
        /* x[0] to x[staNum - 1] is the correlation result */

        int    *mulNum    = ssGetIWork(S) + 1;
        double *iniPhase   = ssGetRWork(S) + 2;
        double mini, tmp, tmpTime;
        int    i, decis, numBase, numInL, numLayer, decLay;
          char   err_msg[256];

        numLayer = mxGetPr(NUM_LAYER)[0];
        /* first find the best amplitude match, find the minimum distance */
        mini = 1000000; decLay = 0;
        for (i = 0; i < numLayer; i++) {
            tmp = fabs(u[1] - mxGetPr(AMP_IN_L)[i]);
            if (mini >= tmp) {
                decLay = i;
                mini = tmp;
            }
        }
/*
#ifdef MATLAB_MEX_FILE
        sprintf(err_msg, "decLay %d, mini %6.3f, u[1] %6.3f numLayer %d\n",
            decLay, mini, u[1], numLayer);
        mexPrintf(err_msg);
#endif
 */
        /*then find the phase match inside the amplitude range, find the max correlation*/
        /* base number, counting from the very end */
        numBase = 0;
        i = numLayer - 1;
        while( i > decLay) {
            numBase += mxGetPr(NUM_IN_L)[i];
            i--;
        }
        decis = *mulNum - numBase;
        numInL = mxGetPr(NUM_IN_L)[decLay];
        if (numInL > 1) {
            double  maxi;
            int     *phaseInd  = ssGetIWork(S) + 2 + numBase;
            int    phaseIndex;

            maxi = 0.0;
            for (i = 0; i < numInL; i++) {
                phaseIndex = *phaseInd;
                *phaseInd++;
                tmp = x[phaseIndex];
/*
#ifdef MATLAB_MEX_FILE
                sprintf(err_msg, "phaseIndex %d, numBase %d, x[?] %6.3f, tmp %6.3f, maxi %6.3f, iniPhase %6.3f\n",
                phaseIndex, numBase, x[phaseIndex], tmp, maxi, iniPhase[phaseIndex]);
                mexPrintf(err_msg);
#endif
 */
                if (tmp > maxi) {
                    decis = *mulNum - numBase - i;
                    maxi = tmp;
                }
/*
#ifdef MATLAB_MEX_FILE
                sprintf(err_msg, "decision %d, maxi %6.3f\n",
                decis, maxi);
                mexPrintf(err_msg);
#endif
 */
            }
        }
        *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