00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00041
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
00064 double derdK_alpha, derdK_beta;
00065 vec pdf_GG, pdf_GGbisalpha, pdf_GGbisbeta;
00066
00067
00068 step_alpha = *palpha / 10;
00069 step_beta = *pbeta / 10;
00070 last_dK = 666;
00071
00072
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
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
00101 }
00102 while (dK + IT_EPSILON * 1e3 < last_dK);
00103
00104 step_alpha *= 0.5;
00105 step_beta *= 0.5;
00106 }
00107 }