surandom.c

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

/*
 * SURANDOM   SIMULINK uniform random number generator.
 *
 *      Syntax:  [sys, x0] = SURANDOM(t, x, u, flag, seed, up_bd, lo_bd)
 *      This function output random numbers. The numbers are uniformly
 *      distributed between the low boundary and the high boundary.
 *      The output is a vector when the seed is a vector.
 *
 * Wes Wang  August 25, 1994
 * Copyright (c) 1994-96 by The MathWorks, Inc.
 * All Rights Reserved
 * $Revision: 1.1 $  $Date: 1996/04/01 19:06:01 $
 */

#define S_FUNCTION_NAME surandom


#ifdef MATLAB_MEX_FILE
#include <stdio.h>    /* needed for declaration of sprintf */
#include "mex.h"      /* needed for declaration of mexErrMsgTxt */
#endif

/*
 * need to include simstruc.h for the definition of the SimStruct and
 * its associated macro definitions.
 */

#include "simstruc.h"

/*
 * Defines for easy access of the input parameters
 */

#define NUM_ARGS       3
#define INISEED        ssGetArg(S,0)
#define UPBOUND        ssGetArg(S,1)
#define LOBOUND        ssGetArg(S,2)

#define LOW_SEED   1          /* minimum seed value */
#define HIGH_SEED  2147483646 /* maximum seed value */
#define START_SEED 1144108930 /* default seed value */

/*
 * mdlInitializeSizes - called to initialize the sizes array stored in
 *                      the SimStruct.  The sizes array defines the
 *                      characteristics (number of inputs, outputs,
 *                      states, etc.) of the S-Function.
 */

static double urand(seed)
	unsigned int *seed;
{
	unsigned int hi, lo;
#if (UINT_MAX == 0xffffffff)
	int test;
#else
	long test;
#endif

#define IA 16807	/* magic multiplier = 7^5 */
#define IM 2147483647	/* modulus = 2^31 - 1 */
#define IQ 127773	/*IM div IA */
#define IR 2836	/* IM modulo IA */
#define S  4.656612875245797e-10	/* reciprocal of 2^31-1 */

	hi = *seed / IQ;
	lo = *seed % IQ;
	test = IA * lo - IR * hi;

	*seed = ((test < 0) ? (unsigned int)(test + IM) : (unsigned int)test);

	return ((double) (*seed * S));

#undef IA
#undef IM
#undef IQ
#undef IR
#undef S
}

static void mdlInitializeSizes(S)
  SimStruct *S;
{
    /*
     * Set-up size information.
     */ 
    
    if (ssGetNumArgs(S) == NUM_ARGS) {
	int numOutput, iniSeed;

	numOutput = mxGetN(INISEED) * mxGetM(INISEED);
	if (numOutput < 1) {
#ifdef MATLAB_MEX_FILE
	    char err_msg[256];
	    sprintf(err_msg, "Input variable is empty.");
	    mexErrMsgTxt(err_msg);
#endif    
        }
	if ((mxGetN(UPBOUND) * mxGetM(UPBOUND)) != (mxGetN(LOBOUND) * mxGetM(LOBOUND))) {
#ifdef MATLAB_MEX_FILE
	    char err_msg[256];
	    sprintf(err_msg, "The dimension for the upbound and the low bound should be the same.");
	    mexErrMsgTxt(err_msg);
#endif    
	}

	if (((mxGetN(UPBOUND) * mxGetM(UPBOUND)) != 1) && (mxGetN(UPBOUND) * mxGetM(UPBOUND)) != numOutput) {
#ifdef MATLAB_MEX_FILE
	    char err_msg[256];
	    sprintf(err_msg, "Dimension for boundaries must match the dimension for seeds.");
	    mexErrMsgTxt(err_msg);
#endif    
	}
    ssSetNumContStates(    S, 0);
	ssSetNumDiscStates(    S, 0);
	ssSetNumInputs(        S, 0);
	ssSetNumOutputs(       S, numOutput);
	ssSetDirectFeedThrough(S, 0);
	ssSetNumInputArgs(     S, NUM_ARGS);
	ssSetNumSampleTimes(   S, 1);
	ssSetNumRWork(         S, numOutput);
	ssSetNumIWork(         S, numOutput);
	ssSetNumPWork(         S, 0);
    } 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 - initializes the array of sample times stored in
 *                            the SimStruct associated with this S-Function.
 */

static void mdlInitializeSampleTimes(S)
  SimStruct *S;
{
    /*
     * Note, blocks that are continuous in nature should have a single
     * sample time of 0.0.
     */

    ssSetSampleTimeEvent(S, 0, 0.0);
    ssSetOffsetTimeEvent(S, 0, 0.0);
}

/*
 * mdlInitializeConditions - initializes the states for the S-Function
 */

static void mdlInitializeConditions(x0, S)
  double *x0;
  SimStruct *S;
{
    int     numOutput, i;
    double  *RWork         = ssGetRWork(S);       /* Real Work Vector */
    int     *IWork         = ssGetIWork(S);       /* Integer Work Vector */
    double  *Pseed         = mxGetPr(INISEED);
    double  *Lo_Bd         = mxGetPr(LOBOUND);
    double  *Up_Bd         = mxGetPr(UPBOUND);

    numOutput = mxGetN(INISEED) * mxGetM(INISEED);

    for (i = 0; i < numOutput; i++) {
      IWork[i] = RWork[i];
      if ((IWork[i] < LOW_SEED) || (IWork[i] > HIGH_SEED))
	  IWork[i] = START_SEED + i;
      RWork[i] = urand((unsigned int *) &IWork[i]);
    }

    numOutput = mxGetN(UPBOUND) * mxGetM(UPBOUND);
    for (i = 0; i < numOutput; i++) {
      if ((Up_Bd[i] - Lo_Bd[i]) < 0) {
#ifdef MATLAB_MEX_FILE
	mexErrMsgTxt("Uniform random distribution low bound cannot be larger than up bound.");
#endif
      }
    }
}

/*
 * mdlOutputs - computes the outputs of the S-Function
 */

static void mdlOutputs(y, x, u, S, tid)
  double *y, *x, *u;
  SimStruct *S;
  int tid;
{
    int     numOutput, i, mFlag;
    double  *RWork         = ssGetRWork(S);       /* Real Work Vector */
    int     *IWork         = ssGetIWork(S);       /* Integer Work Vector */
    double  *Lo_Bd         = mxGetPr(LOBOUND);
    double  *Up_Bd         = mxGetPr(UPBOUND);

    numOutput = mxGetN(INISEED) * mxGetM(INISEED);

    if (mxGetN(UPBOUND) * mxGetM(UPBOUND) <= 1) {
      mFlag = 1;
    } else {
      mFlag = 0;
    }

    for (i = 0; i < numOutput; i++) {
      y[i] = RWork[i];
      if (mFlag) {
	     y[i] = y[i] * (Up_Bd[0] - Lo_Bd[0]) + Lo_Bd[0];
      } else {
	     y[i] = y[i] * (Up_Bd[i] - Lo_Bd[i]) + Lo_Bd[i];
      }      
      RWork[i] = urand((unsigned int *) &IWork[i]);
    }
}

/*
 * mdlUpdate - computes the discrete states of the S-Function
 */

static void mdlUpdate(x, u, S, tid)
  double  *x, *u;
  SimStruct *S;
  int tid;
{
}

/*
 * mdlDerivatives - computes the derivatives of the S-Function
 */

static void mdlDerivatives(dx, x, u, S, tid)
  double *dx, *x, *u;
  SimStruct *S;
  int tid;
{
}

/*
 * mdlTerminate - called at termination of model execution.
 */

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