src/source.c

Go to the documentation of this file.
00001 /*
00002    libit - Library for basic source and channel coding functions
00003    Copyright (C) 2005-2005 Vivien Chappelier, Herve Jegou
00004 
00005    This library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Library General Public
00007    License as published by the Free Software Foundation; either
00008    version 2 of the License, or (at your option) any later version.
00009 
00010    This library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Library General Public License for more details.
00014 
00015    You should have received a copy of the GNU Library General Public
00016    License along with this library; if not, write to the Free
00017    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00018 */
00019 
00020 /*
00021   Source generators
00022   Copyright (C) 2005 Vivien Chappelier, Herve Jegou
00023 */
00024 
00025 
00026 #include <it/random.h>
00027 #include <it/source.h>
00028 #include <it/source_func.h>
00029 #include <it/io.h>
00030 #include <it/math.h>
00031 
00032 /* generate a vector of random binary values with the given */
00033 /* probability of the 0 symbol */
00034 bvec source_binary (idx_t size, double p0)
00035 {
00036   bvec v = bvec_new (size);
00037   idx_t i;
00038 
00039   for (i = 0; i < size; i++)
00040     v[i] = (it_rand () >= p0);
00041 
00042   return (v);
00043 }
00044 
00045 
00046 /* generate a vector of values uniformly distributed in [a,b[ */
00047 vec source_uniform (idx_t size, double a, double b)
00048 {
00049   vec  v = vec_new (size);
00050   idx_t i;
00051 
00052   for (i = 0; i < size; i++)
00053     v[i] = a + (b - a) * it_rand ();
00054 
00055   return (v);
00056 }
00057 
00058 /* generate a vector of values uniformly distributed in [a,b[ */
00059 ivec source_uniform_int (idx_t size, int a, int b)
00060 {
00061   ivec  v = ivec_new (size);
00062   idx_t i;
00063 
00064   for (i = 0; i < size; i++)
00065     v[i] = floor ((double)a + (b - a) * it_rand ());
00066 
00067   return (v);
00068 }
00069 
00070 /* generate a vector of independent values drawn from a */
00071 /* gaussian distribution of given mean and standard deviation */
00072 vec source_gaussian (idx_t size, double mean, double std)
00073 {
00074   vec  v = vec_new (size);
00075   idx_t i;
00076 
00077   for (i = 0; i < size; i++)
00078     v[i] = mean + std * it_randn ();
00079 
00080   return (v);
00081 }
00082 
00083 
00084 /* generate a stationnary random vector from a probability
00085    density function using the acceptance-rejection method.
00086    the pdf is assumed to be zero outside [a, b].
00087 */
00088 vec source_pdf (idx_t size, double a, double b, it_function_t pdf,
00089     it_args_t args)
00090 {
00091   idx_t i;
00092   vec  v = vec_new (size);
00093 
00094   for (i = 0; i < size; i++)
00095     v[i] = it_randpdf (a, b, pdf, args);
00096 
00097   return (v);
00098 }
00099 
00100 
00101 /*--------------------------------------------------------------------*/
00102 ivec source_memoryless (idx_t size, vec pdf)
00103 {
00104   idx_t t;
00105   double x;
00106   int  s, e, h;
00107   vec  pdfcs;
00108   ivec S = ivec_new (size);
00109 
00110   it_assert (is_valid_pdf (pdf, IT_EPSILON), "Invalid probability law");
00111   pdfcs = vec_cum_sum (pdf);
00112   vec_ins (pdfcs, 0, 0.0);
00113 
00114   for (t = 0; t < size; t++) {
00115     /* generate a uniformly distributed value  */
00116     /* and search the integral of the pdf for  */
00117     /* the abciss corresponding to that random */
00118     /* value.                                  */
00119     x = it_rand ();
00120     s = 0;
00121     e = vec_length (pdfcs) - 1;
00122 
00123     while (e > s + 1) {
00124       /* compare with the median of the range */
00125       /* and narrow down the range. */
00126       h = (s + e) / 2;
00127 
00128       if (x < pdfcs[h])
00129   e = h;
00130       else
00131   s = h;
00132     }
00133 
00134     if (x < pdfcs[e])
00135       S[t] = s;
00136     else
00137       S[t] = e;
00138   }
00139   vec_delete (pdfcs);
00140   return S;
00141 }

Hosted by
Copyright (C) 2005-2006 Hervé Jégou
Vivien Chappelier
Francois Cayre
libit logo courtesy of Jonathan Delhumeau