I was amazed when use randn command at Matlab. randn command will generate random data every we call that command. After I search at google, I found how to make this happen. I get this code at seismic unix source code. This code will generate random noise or white noise with Gaussian method. Code for main.c is :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
#include <stdio.h> #include <stdlib.h> #include "frannor.h" int main(int argc, char **argv) { double *noise=NULL; int i; int ndata; unsigned int seed; if(argc!=2) { fprintf(stderr, "usage : random 5"); exit(0); } ndata = atoi(argv[1]); if (-1 == (seed = (unsigned int) time((time_t *) NULL))) { fprintf(stderr, "time() failed to set seed"); exit(0); } srannor(seed); //seed random number generator noise = (double*) calloc (ndata, sizeof(double)); for(i=0; i<ndata; i++) //create random data { /* Compute noise vector elements in [-1, 1] */ /* GAUSS METHOD. frannor gives elements in N(0,1)--ie. pos & negs */ noise[i] = (double) frannor(); } for(i=0; i<ndata; i++) printf("%3.2f ", noise[i]); printf("\n"); free(noise); //free allocated memory return(1); } |
This is code for frannor.h :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
#ifndef FRANNOR_H_ #define FRANNOR_H_ #include <math.h> #define ABS(x) ((x) < 0 ? -(x) : (x)) #define AA 12.37586 #define Bsu 0.4878992 #define Csu 12.67706 #define C1 0.9689279 #define C2 1.301198 #define PC 0.01958303 #define XN 2.776994 #define OXN 0.3601016 #define NBITS 24 /* macro to generate a random number uni uniform on [0,1) */ #define UNI(uni) \ uni = u[i]-u[j]; if (uni<0.0) uni += 1.0; u[i] = uni; \ if (--i<0) i = 16; if (--j<0) j = 16 float frannor(void); void srannor (int seed); #endif /* FRANNOR_H_ */ |
This is code for frannor.c :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 |
/* Copyright (c) Colorado School of Mines, 2010.*/ /* All rights reserved. */ /*********************** self documentation **********************/ /***************************************************************************** FRANNOR - functions to generate a pseudo-random float normally distributed with N(0,1); i.e., with zero mean and unit variance. frannor return a normally distributed random float srannor seed random number generator for normal distribution ****************************************************************************** Function Prototypes: float frannor (void); void srannor (int seed); ****************************************************************************** frannor: Input: (none) Returned: normally distributed random float srannor: Input: seed different seeds yield different sequences of random numbers. ****************************************************************************** Notes: Adapted from subroutine rnor in Kahaner, Moler and Nash (1988) which in turn was based on an algorithm by Marsaglia and Tsang (1984). ****************************************************************************** References: "Numerical Methods and Software", D. Kahaner, C. Moler, S. Nash, Prentice Hall, 1988. Marsaglia G. and Tsang, W. W., 1984, A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions: SIAM J. Sci. Stat. Comput., v. 5, no. 2, p. 349-359. ***************************************************************************** Author: Dave Hale, Colorado School of Mines, 01/21/89 *****************************************************************************/ /**************** end self doc ********************************/ #include "frannor.h" static float v[]={ 0.3409450, 0.4573146, 0.5397793, 0.6062427, 0.6631691, 0.7136975, 0.7596125, 0.8020356, 0.8417227, 0.8792102, 0.9148948, 0.9490791, 0.9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917, 1.1332738, 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635, 1.2982650, 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929, 1.4569915, 1.4834526, 1.5100121, 1.5367061, 1.5635712, 1.5906454, 1.6179680, 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422, 1.7896223, 1.8200099, 1.8510770, 1.8829044, 1.9155830, 1.9492166, 1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713, 2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117, 2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943 }; /* internal state variables for uniform random number generator */ static int i=16,j=4; static float u[]={ 0.8668672834288, 0.3697986366357, 0.8008968294805, 0.4173889774680, 0.8254561579836, 0.9640965269077, 0.4508667414265, 0.6451309529668, 0.1645456024730, 0.2787901807898, 0.06761531340295, 0.9663226330820, 0.01963343943798, 0.02947398211399, 0.1636231515294, 0.3976343250467, 0.2631008574685 }; float frannor(void) /************************************************************************** return a normally distributed random float *************************************************************************** Returned: normally distributed random float **************************************************************************/ { int k; float uni,vni,rnor,x,y,s,bmbx,xnmx; /* uni is uniform on [0,1) */ UNI(uni); /* vni is uniform on [-1,1) */ vni = uni+uni-1.0; /* k is in range [0,63] */ k = ((int)(u[i]*128))%64; /* fast part */ rnor = vni*v[k+1]; if (ABS(rnor)<=v[k]) return rnor; /* slow part */ x = (ABS(rnor)-v[k])/(v[k+1]-v[k]); UNI(y); s = x+y; if (s<=C2) { if (s<=C1) return rnor; bmbx = Bsu-Bsu*x; if (y<=Csu-AA*exp(-0.5*bmbx*bmbx)) { if (exp(-0.5*v[k+1]*v[k+1])+y*PC/v[k+1] <= exp(-0.5*rnor*rnor)) return rnor; do { UNI(y); x = OXN*log(y); UNI(y); } while (-2.0*log(y)<=x*x); xnmx = XN-x; return (rnor>=0.0 ? ABS(xnmx) : -ABS(xnmx)); } } bmbx = Bsu-Bsu*x; return (rnor>=0.0 ? ABS(bmbx) : -ABS(bmbx)); } void srannor (int seed) /***************************************************************************** seed random number generator ****************************************************************************** Input: seed different seeds yield different sequences of random numbers. *****************************************************************************/ { int ii,jj,ia,ib,ic,id; float s,t; i = 16; j = 4; ia=ABS(seed)%32707; ib=1111; ic=1947; for (ii=0; ii<17; ii++) { s = 0.0; t = 0.5; for (jj=0; jj<64; jj++) { id = ic-ia; if (id<0) { id += 32707; s += t; } ia = ib; ib = ic; ic = id; t *= 0.5; } u[ii] = s; } } |
To compile this source code, use command :
gcc main.c frannor.c -lm -o random
This is output from program :
|
1 2 3 4 |
toto@toto-laptop:/home/toto/$ ./random 6 0.08 1.78 -1.35 -1.03 1.62 0.07 toto@toto-laptop:/home/toto/$ ./random 6 1.12 1.32 0.73 -0.34 0.25 -0.58 |
Source : http://www.cwp.mines.edu/cwpcodes/index.html
My name is Toto Sugito. This is my notes when I try something. Maybe, this is NOT THE BEST SOLUTION for your problems. If you have other method or idea that better with my post, please share in this blog. Thank for visiting my site.