src/source_fit.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 fitting functions
00022   Copyright (C) 2005 Vivien Chappelier, Herve Jegou
00023 */
00024 
00025 
00026 #include <math.h>
00027 #include <it/source_fit.h>
00028 #include <it/distance.h>
00029 #include <it/math.h>
00030 #include <it/io.h>
00031 
00032 /*--------------------------------------------------------------------*/
00033 static double GenGaussian (double alpha, double beta, double x)
00034 {
00035   return exp (-pow (fabs (x) / alpha, beta));
00036 }
00037 
00038 
00039 /*--------------------------------------------------------------------*/
00040 /* Return the discrete pdf law associated with the Generalized 
00041    Gaussian of parameters gamma, beta and C                           */
00042 vec source_pdf_GG (vec symbols, double alpha, double beta)
00043 {
00044   idx_t i;
00045   vec  pdf = vec_new_zeros (vec_length (symbols));
00046   for (i = 0; i < vec_length (symbols); i++)
00047     pdf[i] = GenGaussian (alpha, beta, symbols[i]);
00048 
00049   vec_normalize (pdf, 1);
00050   return pdf;
00051 }
00052 
00053 
00054 
00055 /*--------------------------------------------------------------------*/
00056 void source_estim_GG_from_histo (vec pdf, vec symbols, double *palpha,
00057          double *pbeta)
00058 {
00059   double step_alpha, step_beta;
00060   double dK, last_dK;
00061   int  i;
00062 
00063   /* derivatives of the Kullback distance on alpha and beta */
00064   double derdK_alpha, derdK_beta;
00065   vec  pdf_GG, pdf_GGbisalpha, pdf_GGbisbeta;
00066 
00067   /* Kind of a gradient descent with the derivative of kullback distance.   */
00068   step_alpha = *palpha / 10;
00069   step_beta = *pbeta / 10;
00070   last_dK = 666;
00071 
00072   /* Initial Kullback distance */
00073   pdf_GG = source_pdf_GG (symbols, *palpha, *pbeta);
00074   dK = vec_distance_kullback_leibler (pdf, pdf_GG);
00075   vec_delete (pdf_GG);
00076 
00077   /* "Geometric" decreasing of the convergence step */
00078   for (i = 0; i < 5; i++) {
00079     do {
00080       pdf_GGbisalpha = source_pdf_GG (symbols, *palpha + IT_EPSILON, *pbeta);
00081       derdK_alpha =
00082   (vec_distance_kullback_leibler (pdf, pdf_GGbisalpha) -
00083    dK) / IT_EPSILON;
00084       vec_delete (pdf_GGbisalpha);
00085 
00086       pdf_GGbisbeta = source_pdf_GG (symbols, *palpha, *pbeta + IT_EPSILON);
00087       derdK_beta =
00088   (vec_distance_kullback_leibler (pdf, pdf_GGbisbeta) -
00089    dK) / IT_EPSILON;
00090       vec_delete (pdf_GGbisbeta);
00091 
00092       *palpha = *palpha - derdK_alpha * step_alpha;
00093       *pbeta = *pbeta - derdK_beta * step_beta;
00094 
00095       last_dK = dK;
00096 
00097       pdf_GG = source_pdf_GG (symbols, *palpha, *pbeta);
00098       dK = vec_distance_kullback_leibler (pdf, pdf_GG);
00099       vec_delete (pdf_GG);
00100       //              printf( "alpha = %f, beta = %f, dK = %f\n", *palpha, *pbeta, dK );
00101     }
00102     while (dK + IT_EPSILON * 1e3 < last_dK);
00103 
00104     step_alpha *= 0.5;
00105     step_beta *= 0.5;
00106   }
00107 }

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