[FFmpeg-devel] [PATCH] Box-Muller gaussian generator

[Michael Niedermayer]

On Fri, Dec 11, 2009 at 08:38:35PM +0100, Stefano Sabatini wrote:
> On date Friday 2009-12-11 13:23:54 +0100, Michael Niedermayer encoded:
> > On Fri, Dec 11, 2009 at 01:11:29PM +0100, Stefano Sabatini wrote:
> > > On date Friday 2009-12-11 01:49:21 +0100, Michael Niedermayer encoded:
> > > > On Fri, Dec 11, 2009 at 01:40:08AM +0100, Michael Niedermayer wrote:
> > > > > On Sun, Dec 06, 2009 at 11:13:37PM +0100, Stefano Sabatini wrote:
> > > > > > On date Sunday 2009-12-06 19:40:20 +0100, Vitor Sessak encoded:
> > > > > > > Stefano Sabatini wrote:
> > > > > > >> Hi,
> > > > > > >>
> > > > > > >> I'm going to use this for a noise filter, I believe this could be
> > > > > > >> useful in other part of libav* though.
> > > > > > >
> > > > > > >> +double av_bmg_get(AVBMG *bmg)
> > > > > > >> +{
> > > > > > >> +    double x1, x2, w, y;
> > > > > > >> +    AVLFG *lfg = &bmg->lfg;
> > > > > > >> +
> > > > > > >> +    do {
> > > > > > >> +        x1 = 2.0 * ranf(lfg) - 1.0;
> > > > > > >> +        x2 = 2.0 * ranf(lfg) - 1.0;
> > > > > > >> +        w = x1*x1 + x2*x2;
> > > > > > >> +    } while (w >= 1.0);
> > > > > > >> +
> > > > > > >> +    w = sqrt((-2.0 * log(w)) / w);
> > > > > > >> +    y = x2 * w;
> > > > > > >> +
> > > > > > >> +    return bmg->mean + y * bmg->std_dev;
> > > > > > >> +}
> > > > > > >
> > > > > > > This generates two numbers to use only one. Isn't it faster to store the  
> > > > > > > second number for the next call or have a function to fill a buffer of  
> > > > > > > random numbers?
> > > > > > 
> > > > > > Implemented the first idea.
> > > > > 
> > > > > void av_bmg_get(AVLFG *lfg, double out[2]);
> > > > >     double x1, x2, w;
> > > > > 
> > > > >     do {
> > > > >         x1 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
> > > > >         x2 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
> > > > >         w = x1*x1 + x2*x2;
> > > > >     } while (w >= 1.0);
> > > > >
> > > > >     w = sqrt((-2.0 * log(w)) / w);
> > > 
> > > If av_lfg_get() returns UINT_MAX/2 two times in a row then we have:
> > > x1 = x2 = w = 0
> > > while this isn't very likely (I don't even know if av_lfg() can return
> > > two times in a row the same number) we'll have a division-by-zero.
> > 
> > > Is this the bug you mentioned in your first mail?
> > 
> > yes
> > 
> > 
> > > Also how it is fixed in this second variant?
> > 
> > UINT_MAX/2.0 == 0xFFFFFFFF/2.0 is not an integer, av_lfg_get() returns
> > integers
> 
> Updated patch.
[...]


> +#undef printf
> +
[...]
>      }
>      av_log(NULL, AV_LOG_ERROR, "final value:%X\n", x);
[...]
> +            printf("%f\n%f\n",

please dont do this


[...]
> Index: libavfilter-soc/ffmpeg/libavutil/lfg.h
> ===================================================================
> --- libavfilter-soc.orig/ffmpeg/libavutil/lfg.h	2009-12-11 20:28:37.000000000 +0100
> +++ libavfilter-soc/ffmpeg/libavutil/lfg.h	2009-12-11 20:35:26.000000000 +0100
> @@ -51,4 +51,12 @@
>      return c->state[c->index++ & 63] = 2*a*b+a+b;
>  }
>  
> +/**
> + * Returns the next number generated by a Box-Muller generator using
> + * the random number issued by lfg.
> + *
> + * @param out[2] array where are placed the two generated numbers

I think you should mention something about gaussian distributions

[...]
-- 
Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB

Into a blind darkness they enter who follow after the Ignorance,
they as if into a greater darkness enter who devote themselves
to the Knowledge alone. -- Isha Upanishad
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
URL: <http://lists.mplayerhq.hu/pipermail/ffmpeg-devel/attachments/20091211/ce6749ac/attachment.pgp>