srandbit.c

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

/*========================================================================
 *     Syntax: [sys, x0, str, ts] = srandbit(t, x, u, flag, n, prob, seed)
 *
 * SRANDBIT.c --- file name
 * COMMU TOOLBOX : Outputs a zero vector with random bit ones as Source.
 * SRANDDBIT outputs a zero vector with random bit ones.
 *           This block outputs outputs a zero vector with ones in the vector.
 *       	The probability of ones in the vector is controlled by the input
 *       	vector prob.
 *
 * Version 1.0
 *
 * Original designed by Wes Wang,
 * Jun Wu,  The Mathworks, Inc.
 * Nov. 14, 1995

 * Copyright (c) 1995-96 by The MAthWorks, Inc.
 * All Rights Reserved
 * $Revision: 1.1 $  $Date: 1996/04/01 19:05:20 $
 *========================================================================
 *
 * This function has no input and one output. The output is a zero vector
 * with ones in the vector. The probability of ones in the vector is
 * controlled by the input vector prob[].
 *
 * n    is the element number of output vector. 
 * prob is the probability of ones in the output vector. The sum of all
 *           elements in prob cannot be greater than one.
 * seed is the seed of the random number generator.
 *
 *========================================================================*/
#define S_FUNCTION_NAME srandbit

#ifdef MATLAB_MEX_FILE
#include <stdio.h>      
#endif
#include <stdlib.h>     

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

#define NUM_ARGS    3           /* three additional input arguments */
#define N       ssGetArg(S,0)   /* the length of  output vector */
#define PROB    ssGetArg(S,1)   /* the probability of ones in output vector */
#define SEED    ssGetArg(S,2)   /* the seed of output */

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

/*
 * uniform mxRandom number generator
 */
static double Urand(seed)  		/* pointer to a running seed */
  unsigned int *seed;
{
#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	*/

	    unsigned int hi, lo;

    	hi = *seed / IQ;
	    lo = *seed % IQ;
	    test = IA * lo - IR * hi;	/* never overflows	*/
	
	    *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;
{
    /* print warning messages. */  
    int     i;
    int     NumOutput = mxGetPr(N)[0];
    int     len_prob= mxGetM(PROB) * mxGetN(PROB);
    double  sum=0.0;

    for ( i = 0;   i < len_prob; i++)
        sum = sum + mxGetPr(PROB)[i];


    ssSetNumContStates(     S, 0);          /* number of continuous states */
    ssSetNumDiscStates(     S, 0);          /* number of discrete states */
    ssSetNumOutputs(        S, NumOutput);  /* number of outputs */
    ssSetNumInputs(         S, 0);          /* number of inputs */
    ssSetDirectFeedThrough( S, 0);          /* direct feedthrough flag */
    ssSetNumSampleTimes(    S, 1);          /* number of sample times */
    ssSetNumInputArgs(      S, NUM_ARGS);   /* number of input arguments */
    ssSetNumRWork(          S, 2*NumOutput+2);
                    /*  1st: dividing [0 -- len_prob];
                     *  2nd: div_zeros[len_prob+1 -- 2*len_prob];
                     *  3rd: Rand[2*len_prob+1 -- last one];
                     */
    ssSetNumIWork(          S, 3*NumOutput+1);
                    /*  1st: indx[0 -- NumOutput-1];
                     *  2nd: indx_zeros[NumOutput -- 2*NumOutput-1];
                     *  3rd: ind_loc[2*NumOutput -- 3*NumOutput-1];
                     *  4th: Seed[3*NumOutput -- last one];
                     */
    ssSetNumPWork(          S, 0);
}

/*
 * mdlInitializeConditions - initialize the states
 * Initialize the states, Integers and real-numbers
 */
static void mdlInitializeConditions(x0, S)
  double *x0;
  SimStruct *S;
{
    int     i;    
    int     len_prob = mxGetM(PROB) * mxGetN(PROB);
    int     NumOutput = mxGetPr(N)[0];

    int     *indx = ssGetIWork(S);
    int     *ind_zeros = ssGetIWork(S)+NumOutput;
    int     *ind_loc = ssGetIWork(S)+2*NumOutput;
    int     *Seed = ssGetIWork(S)+3*NumOutput;

    double  *dividing = ssGetRWork(S);
    double  *div_zeros= ssGetRWork(S)+NumOutput+1;
    double  *Rand = ssGetRWork(S)+2*NumOutput+1;

    double  *Pseed = mxGetPr(SEED);

    dividing[len_prob] = 0.0;
    for (i=0; i < len_prob; i++){
        dividing[i] = 0.0;
        div_zeros[i]= 0.0;
        indx[i] = 0;
    }
    for (i=0; i < NumOutput; i++){
        ind_zeros[i] = 0;
        ind_loc[i] = 0;
	}

    if ( Pseed[0] < LOW_SEED || Pseed[0] > HIGH_SEED)
        Seed[0] = START_SEED;
    else
        Seed[0] = (int)Pseed[0];
}

/*
 * 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);
}

/*
 * 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;
{
    int     NumOutput = mxGetPr(N)[0];
    int     len_prob = mxGetM(PROB) * mxGetN(PROB);
    double  *prob    = mxGetPr(PROB);       /* inputed probability by user */    

    int     *indx = ssGetIWork(S);
    int     *ind_zeros = ssGetIWork(S)+NumOutput;
    int     *ind_loc = ssGetIWork(S)+2*NumOutput;
    int     *Seed = ssGetIWork(S)+3*NumOutput;

    double  *dividing = ssGetRWork(S);
    double  *div_zeros= ssGetRWork(S)+NumOutput+1;
    double  *Rand = ssGetRWork(S)+2*NumOutput+1;

    int     i, j;
    int     len_indx=0, len_ind_zeros=0, len_ind_loc=0;
    
    for ( i=len_prob-1; i > -1; i-- )
        dividing[i] = dividing[i+1] + prob[i];

    for (i=0; i<NumOutput; i++)
        y[i] = 0;

    len_indx = 0;       
    Rand[0] = Urand((unsigned int *) &Seed[0]);
    for ( i=0; i < len_prob; i++ ){
         if (dividing[i] >= Rand[0]){        
            indx[len_indx] = i+1;
            len_indx ++;
         }
    }
    if ( len_indx != 0 ){ 
        for (i=0; i < len_indx; i++){         
            len_ind_zeros = 0;
            for (j=0; j < NumOutput; j++){
                if (y[j] <= 0.5){
	                ind_zeros[len_ind_zeros] = j;
	                len_ind_zeros++;
                }
            }
            if ( len_ind_zeros >= 1 ){
                for (j=0; j < len_ind_zeros; j++)
                    div_zeros[j] =  (double)j/(double)len_ind_zeros;

                Rand[0] = Urand((unsigned int *) &Seed[0]);
                len_ind_loc = 0;
                for (j=0; j < len_ind_zeros; j++){
                    if (div_zeros[j] <= Rand[0])
                        len_ind_loc++;
                }
                j = ind_zeros[len_ind_loc];
    	        y[j] = 1;
			}
		}
	}
}
/*
 * 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;
{
}

/*
 * 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;
{
}

/*
 * 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