开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 嵌入式/单片机编程 > DSP编程 > 查看源码
DHMM_GL.cpp
资源名称：VQ-DHMM.rar　[点击查看]
上传用户：avbj512
上传日期：2013-09-18
资源大小：6239k
文件大小：25k
源码类别：

DSP编程

开发平台：

Visual C++

DHMM_GL.cpp：源码内容
  1. // DHMM_GL.cpp:
  2. // Implementation of the DHMM_GL Module.
  3. // That is the transform of previous DHMM Code by GuLiang.
  4. //
  5. // Created 2001/08, By DongMing, MDSR.
  6. //
  7. //////////////////////////////////////////////////////////////////////
  9. #include "stdafx.h"
  10. #include "DHMM_GL.h"
  12. #include <math.h>
  14. #define REESTIMATE_ALL (0x0001)
  15. #define ENDSTATE_ABSORB (0x0002)
  16. #define REESTIMATE_B (0x0004)
  18. //////////////////////////////////////////////////////////////////////
  19. // Private Function Head
  20. void DHMM_VQ(int *DMNU_Sound_attrib, int DMNU_Speaker_num, char *DMNU_Featurefile_name, char *DMNU_Codefile_name,
  21.  int Codeword_num, double K_means_precision, int DMNU_mode,
  22.  int DMAD_n_Code_Word_Num, int DMAD_n_Code_Word_Dim,
  23.  DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Book, DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Word);
  24. void DHMM_train(int *DMNU_Sound_attrib, int Speaker_num, int Train_num, char *DMNU_Featurefile_name,
  25. char *DMNU_Codefile_name, char *DMNU_Patternfile_name, int Status_num, int mode,
  26. WORD_SAMPLE * DMAD_pu_Word_Sample,
  27. DHMM_MODEL * DMAD_pu_DHMM_Model);
  28. double DHMM_Calculate(double ***Alpha, double ***Belta, double *Pi, double **A,
  29.   double **B, int **Out, int Sample_num, int Output_num, int *T,
  30.   int n, int mode, int End_state_mode);
  31. double DHMM_Priority(double *Pi, double **A, double **B, int *Out, int T, int n, int mode);
  33. //////////////////////////////////////////////////////////////////////
  34. // API functions
  35. int DHMM_VQ_Train_Code_Book_GL(DYNA_2DIM_DOUBLE_ARRAY d2dda_Code_Word, int n_Code_Word_Num, int n_Code_Word_Dim,
  36.    DYNA_2DIM_DOUBLE_ARRAY d2dda_Initial_Code_Book, DYNA_2DIM_DOUBLE_ARRAY d2dda_Code_Book, int n_Code_Book_Size)
  37. {
  38. PRO_LOG("tVQ = GL, K-mean precision = %15.4E.n", KMEAN_PRECISION);
  39. DHMM_VQ(NULL, 0, NULL, NULL,
  40. n_Code_Book_Size, KMEAN_PRECISION, 0,
  41. n_Code_Word_Num, n_Code_Word_Dim,
  42. d2dda_Code_Book, d2dda_Code_Word);
  44. return 0;
  45. }
  47. int DHMM_Model_Train_DHMM_Model_GL(WORD_SAMPLE * pu_Word_Sample, int n_Word_Sample_Num,
  48.    DHMM_MODEL * pu_DHMM_Model)
  49. {
  50. PRO_LOG("tModel = GL, max loop = %4d, precision = %15.4E.n", MAX_CALC_NUM, EPSILON_CALC);
  51. DHMM_train(NULL, n_Word_Sample_Num, n_Word_Sample_Num, NULL,
  52. NULL, NULL, pu_DHMM_Model->n_State_Num, ~ENDSTATE_ABSORB,
  53. pu_Word_Sample, pu_DHMM_Model);
  55. return 0;
  56. }
  58. int DHMM_Recog_Forward_Backward_GL(DHMM_MODEL * pu_DHMM_Model,
  59.    WORD_SAMPLE * pu_Word_Sample,
  60.    double * pd_Max_Likelihood)
  61. {
  62. (*pd_Max_Likelihood) = DHMM_Priority(pu_DHMM_Model->pdPi, pu_DHMM_Model->d2dda_A, pu_DHMM_Model->d2dda_B,
  63. pu_Word_Sample->pn_VQed_Feature_Sequence, pu_Word_Sample->n_Feature_Sequence_Len,
  64. pu_DHMM_Model->n_State_Num, ~ENDSTATE_ABSORB);
  66. return 0;
  67. }
  69. //////////////////////////////////////////////////////////////////////
  70. // Original Code
  71. //////////////////////////////////////////////////////////////////////
  73. // DMNU = DongMing doesn't use this when transforming.
  74. // DMAD = DongMing addes this when transforming.
  76. // binary split VQ code-book training (clustering) algorithm
  77. static void DHMM_VQ(int *DMNU_Sound_attrib, int DMNU_Speaker_num, char *DMNU_Featurefile_name, char *DMNU_Codefile_name,
  78. int Codeword_num, double K_means_precision, int DMNU_mode,
  79. int DMAD_n_Code_Word_Num, int DMAD_n_Code_Word_Dim,
  80. DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Book, DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Word)
  81. {
  82. int    i, j, p, m;
  83. int    Dimension; // Dimension: 特征维数
  84. int    Nearer;
  85. int    *Counter; // 记录各码字聚类到的训练向量个数
  86. int    Total_num; // 训练向量集的总帧数
  87. double **Codeword; // VQ码本
  88. double **Sum;
  89. double **s; // 训练向量集
  90. double Min_Dist;
  91. double Dist; // 某一训练向量和某一码字间的距离
  92. double New_D; // 所有训练向量的量化失真之和？
  93. double Old_D;
  95. // Initial Memory Block
  96. Dimension = DMAD_n_Code_Word_Dim;
  97. Total_num = DMAD_n_Code_Word_Num;
  98. Codeword = DMAD_d2dda_Code_Book;
  100. if((Sum = new double *[Codeword_num]) == NULL)
  101. {
  102. DEBUG_PRINTF("Not enough memory for Sum !n");
  103. ASSERT(0);
  104. }
  105. for(i=0; i<Codeword_num; i++)
  106. {
  107. if((Sum[i] = new double [Dimension]) == NULL)
  108. {
  109. DEBUG_PRINTF("Not enough memory for Sum[%d] !n", i);
  110. ASSERT(0);
  111. }
  112. }
  114. if((Counter  = new int [Codeword_num]) == NULL)
  115. {
  116. DEBUG_PRINTF("Not enough memory for Counter !n");
  117. ASSERT(0);
  118. }
  120. s = DMAD_d2dda_Code_Word;
  122. // k-means Clusteration
  124. for(p=0; p<Dimension; p++) Codeword[0][p] = 0; // 第一个码字清零
  125. for(i=0; i<Total_num; i++)
  126. for(p=0; p<Dimension; p++)
  127. Codeword[0][p] += s[i][p];
  128. for(p=0; p<Dimension; p++) Codeword[0][p] /= Total_num; // 第一个码字置为所有训练向量的平均值
  130. for(m=2; m<=Codeword_num; m*=2) // m: 分裂后码本大小
  131. {
  132. DEBUG_PRINTF("tM = %d.n", m);
  134. for(i=m/2-1; i>=0; i--)
  135. for(j=0; j<Dimension; j++) // initialize Codeword[0],Codeword[1],...,Codeword[m-1] by splitting current code words
  136. {
  137. Codeword[i*2+1][j] = Codeword[i][j] * 0.9; // splitting parameter epsilon set to 0.1
  138. Codeword[i*2][j]   = Codeword[i][j] * 1.1;
  139. }
  141. Old_D = 1.0e20;
  142. New_D = 1.0e15;
  144. while(((Old_D - New_D) / New_D) > K_means_precision)
  145. {
  146. for(j=0; j<m; j++)
  147. {
  148. for(p=0; p<Dimension; p++) Sum[j][p] = 0;
  149. Counter[j] = 0;
  150. }
  152. Old_D = New_D;
  153. New_D = 0;
  155. for(i=0; i<Total_num; i++)
  156. {
  157. Dist = 0;
  158. for(p=0; p<Dimension; p++)
  159. Dist += (Codeword[0][p] - s[i][p]) * (Codeword[0][p] - s[i][p]);
  160. Min_Dist = Dist; // 计算该训练向量与第一个码字间的距离
  161. Nearer = 0;
  162. for(j=1; j<m; j++) // 计算该训练向量与第2到第m个码字间的距离，并求最小值
  163. {
  164. Dist = 0;
  165. for(p=0; p<Dimension; p++)
  166. Dist += (Codeword[j][p] - s[i][p]) * (Codeword[j][p] - s[i][p]);
  167. if(Dist < Min_Dist)
  168. {
  169. Min_Dist = Dist;
  170. Nearer = j; // 记录距离最近的码字标号
  171. }
  172. }
  173. for(p=0; p<Dimension; p++) Sum[Nearer][p] += s[i][p]; // 累加聚类到标号为Nearer的码字的所有训练向量
  174. Counter[Nearer] ++; // 记录聚类到标号为Nearer的码字的训练向量的个数
  175. New_D += Min_Dist;
  176. }
  178. for(j=0; j<m; j++)
  179. {
  180. if(Counter[j] > 0)
  181. for(p=0; p<Dimension; p++)
  182. Codeword[j][p] = Sum[j][p] / Counter[j]; // 用聚类到标号为j的码字的所有训练的平均值更新该码字
  183. }
  184. }
  185. }
  186. for(i=0; i<Codeword_num; i++)
  187. {
  188. delete Sum[i];
  189. }
  190. delete Sum;
  191. delete Counter;
  192. }
  194. static void DHMM_train(int *DMNU_Sound_attrib, int Speaker_num, int Train_num, char *DMNU_Featurefile_name,
  195.    char *DMNU_Codefile_name, char *DMNU_Patternfile_name, int Status_num, int mode,
  196.    WORD_SAMPLE * DMAD_pu_Word_Sample,
  197.    DHMM_MODEL * DMAD_pu_DHMM_Model)
  198. // Speaker_num 说话人个数
  199. // Train_num 训练集人数(Train_num=Speaker_num-Test_num)
  200. // Status_num 状态数
  201. // mode 训练模式
  202. {
  203. int     i, j, k, /* n, */ m;
  204. int     Output_num; // VQ码本大小（码字个数）
  205. int     **Output; // for each word item model being trained, the output sequence (observation sequence) of each utterance of this model. （每个元素保存该帧的码字标号）
  206. // 1st dimension: index of the training utterance; 
  207. // 2nd dimension: index of frame
  208. int     *Frame_num; // 各训练句子的帧数
  209. double  *Pi; // 初始概率
  210. double  **A; // A矩阵：Status_num*Status_num
  211. double  **B; // B矩阵：Status_num*Output_num
  212. double  ***Alpha; // Train_num*帧数*状态数的double矩阵
  213. double  ***Belta; // Train_num*帧数*状态数的double矩阵
  214. double  w, u;
  216. // Generate (A, B, Pi)
  217. Output_num = DMAD_pu_DHMM_Model->n_Code_Book_Size;
  219. Pi = DMAD_pu_DHMM_Model->pdPi; // The above block transform to this line. #DongMing#
  221. // allocate memory for A matrix
  222. A = DMAD_pu_DHMM_Model->d2dda_A; // The above block transform to this line. #DongMing#
  224. // allocate memory for B matrix
  225. B = DMAD_pu_DHMM_Model->d2dda_B; // The above block transform to this line. #DongMing#
  227. //事实上，此处为Alpha和Belta分配了过多的内存，没有必要为所有用于训练的说话人分配Alpha和Belta
  228. if((Alpha = new double **[Train_num]) == NULL)
  229. {
  230. DEBUG_PRINTF("Not enough memory for Alpha !n");
  231. ASSERT(0);
  232. }
  234. if((Belta = new double **[Train_num]) == NULL)
  235. {
  236. DEBUG_PRINTF("Not enough memory for Belta !n");
  237. ASSERT(0);
  238. }
  240. if((Output = new int *[Train_num]) == NULL)
  241. {
  242. DEBUG_PRINTF("Not enough memory for Output !n");
  243. ASSERT(0);
  244. }
  246. if((Frame_num = new int [Train_num]) == NULL)
  247. {
  248. DEBUG_PRINTF("Not enough memory for Frame_num !n");
  249. ASSERT(0);
  250. }
  252. for(i=1; i<Status_num; i++) Pi[i] = 0; // 初始状态概率：第一个状态为1，其余为0
  253. Pi[0] = 1;
  255. // 初始化A矩阵、B矩阵
  256. for(i=0; i<Status_num; i++)
  257. {
  258. for(j=0; j<Status_num; j++) A[i][j] = 0;
  259. if(i < Status_num-1)
  260. {
  261. A[i][i]   = 0.5;
  262. A[i][i+1] = 0.5;
  263. }
  264. else
  265. {
  266. A[i][i]   = 1.0;
  267. }
  268. w = 0;
  269. for(j=0; j<Status_num; j++) w = w + A[i][j];
  270. for(j=0; j<Status_num; j++) A[i][j] = A[i][j] / w; // 归一化A矩阵的每一行
  272. w = 0;
  273. for(j=0; j<Output_num; j++) // 随机初始化B矩阵
  274. {
  275. B[i][j] = 1.0 / Output_num * (0.9 + ((rand()%1000) / 5000.0)); // RAND_MAX is defined as 0x7fff in "stdlib.h"
  276. // B[i][j] = 1.0 / Output_num;
  277. w = w + B[i][j];
  278. }
  279. for(j=0; j<Output_num; j++) B[i][j] = B[i][j] / w; // 归一化B矩阵的每一行
  280. }
  282. k = 0; // k: training utterance counter for each word item model being trained
  283. // 装入该词条的训练数据
  284. for(m=0; m<Speaker_num; m++)
  285. {
  286. Frame_num[k] = DMAD_pu_Word_Sample[k].n_Feature_Sequence_Len;
  287. // ????????????? is it necessary to save alpha and beta for each training utterance ??????????????
  288. if((Alpha[k] = new double *[Frame_num[k]]) == NULL)
  289. {
  290. DEBUG_PRINTF("Not enough memory for Alpha[%d] !n", k);
  291. ASSERT(0);
  292. }
  293. for(i=0; i<Frame_num[k]; i++)
  294. if((Alpha[k][i] = new double [Status_num]) == NULL)
  295. {
  296. DEBUG_PRINTF("Not enough memory for Alpha[%d][%d] !n", k, i);
  297. ASSERT(0);
  298. }
  300. if((Belta[k] = new double *[Frame_num[k]]) == NULL)
  301. {
  302. DEBUG_PRINTF("Not enough memory for Belta[%d] !n", k);
  303. ASSERT(0);
  304. }
  305. for(i=0; i<Frame_num[k]; i++)
  306. if((Belta[k][i] = new double [Status_num]) == NULL)
  307. {
  308. DEBUG_PRINTF("Not enough memory for Belta[%d][%d] !n", k, i);
  309. ASSERT(0);
  310. }
  311. if((Output[k] = new int [Frame_num[k]]) == NULL)
  312. {
  313. DEBUG_PRINTF("Not enough memory for Frame_num[%d] !n", k);
  314. ASSERT(0);
  315. }
  316. for(i=0; i<Frame_num[k]; i++)
  317. Output[k][i] = DMAD_pu_Word_Sample[k].pn_VQed_Feature_Sequence[i];
  318. k ++;
  319. }
  321. w = 1;
  322. m = 0;
  323. while(m < MAX_CALC_NUM)
  324. {
  325. // Alpha和Belta作为事先分配好内存的变量，并未在调用DHMM_Calculate后使用，因此应该在DHMM_Calculate函数内部分配内存，
  326. // 而入口参数只传递大小信息（如训练句子数、帧数、状态数等等），甚至可以不必为Alpha和Belta分配如此多的内存。
  327. u = DHMM_Calculate(Alpha, Belta, Pi, A, B, Output, Train_num, Output_num,
  328. Frame_num, Status_num, REESTIMATE_ALL, mode);
  329. //u = DHMM_Calculate(Alpha, Belta, Pi, A, B, Output, Train_num, Output_num,
  330. //Frame_num, Status_num, REESTIMATE_ALL, mode, DMAD_pu_Word_Sample);
  331. // for (i = 0; i < Status_num; i++) for (j = 0; j < Output_num; j++) PRO_LOG("%dt%dt:%15.4E.n", i, j, B[i][j]);
  333. m ++;
  334. if((u > w) && (fabs((u-w)/w) < EPSILON_CALC)) break;
  335. else w = u;
  336. }
  337. // for (i = 0; i < Status_num; i++) for (j = 0; j < Output_num; j++) TRACE("%dt%dt:%15.4E.n", i, j, B[i][j]);
  338. PRO_LOG("tLoop = %4d.n", m);
  340. // 释放内存
  341. for(k=0; k<Train_num; k++)
  342. {
  343. for(i=0; i<Frame_num[k]; i++)
  344. delete Alpha[k][i];
  345. delete Alpha[k];
  347. for(i=0; i<Frame_num[k]; i++)
  348. delete Belta[k][i];
  349. delete Belta[k];
  351. //delete Output[k];
  352. }
  354. delete Alpha;
  355. delete Belta;
  357. //delete Output;
  358. delete Frame_num;
  359. }
  361. static double DHMM_Calculate(double ***Alpha, double ***Belta, double *Pi, double **A,
  362.  double **B, int **Out, int Sample_num, int Output_num, int *T,
  363.  int n, int mode, int End_state_mode)
  364. // 功能：训练某一词条的模型
  366. // input:
  367. // Pi 初始时刻各状态的概率
  368. // Alpha 前向概率
  369. // Belta 后向概率
  370. // A 初始A矩阵
  371. // B 初始B矩阵
  372. // Out 用于训练的各句子的观察序列（输出序列）
  373. // Sample_num 训练集人数
  374. // Output_num VQ码本大小（码字数，即观察值的个数）
  375. // T 用于训练的各句子的帧数
  376. // n 状态数
  377. // mode 训练模式(REESTIMATE_A, REESTIMATE_B, REESTIMATE_ALL)
  378. // End_state_mode 结尾帧状态模式(ENDSTATE_ABSORB = 0x200)
  380. // output:
  381. // A
  382. // B
  384. // caller(r_dhmm.cpp): u = DHMM_Calculate(Alpha, Belta, Pi, A, B, Output, Train_num, Output_num,
  385. //       Frame_num, Status_num, REESTIMATE_ALL, mode);
  386.  
  387. {
  388. int    i, j, k, m;
  389. double w;
  390. double **New_A;
  391. double **New_B;
  392. double **Ct; // to prevent underflow of Alpha and Belta, Ct is the scaling coefficient
  393. double Residue;
  395. if((New_A = new double *[n]) == NULL)
  396. {
  397. DEBUG_PRINTF("Not enough memory for New_A !n");
  398. ASSERT(0);
  399. }
  400. for(i=0; i<n; i++)
  401. {
  402. if((New_A[i] = new double [n]) == NULL)
  403. {
  404. DEBUG_PRINTF("Not enough memory for New_A[%d] !n", i);
  405. ASSERT(0);
  406. }
  407. }
  409. if((New_B = new double *[n]) == NULL)
  410. {
  411. DEBUG_PRINTF("Not enough memory for New_B !n");
  412. ASSERT(0);
  413. }
  414. for(i=0; i<n; i++)
  415. {
  416. if((New_B[i] = new double [Output_num]) == NULL)
  417. {
  418. DEBUG_PRINTF("Not enough memory for New_B[%d] !n", i);
  419. ASSERT(0);
  420. }
  421. }
  423. if((Ct = new double *[Sample_num]) == NULL)
  424. {
  425. DEBUG_PRINTF("Not enough memory for Ct !n");
  426. ASSERT(0);
  427. }
  428. for(i=0; i<Sample_num; i++)
  429. if((Ct[i] = new double [T[i]]) == NULL)
  430. {
  431. DEBUG_PRINTF("Not enough memory for Ct[%d] !n", i);
  432. ASSERT(0);
  433. }
  435. // Calculate Alpha & Belta.
  437. for(m=0; m<Sample_num; m++)
  438. {
  439. if (T[m] > 0)
  440. {
  441. w = 0;
  442. for(i=0; i<n; i++)
  443. {
  444. Alpha[m][0][i] = Pi[i] * B[i][ Out[m][0] ];
  445. w = w + Alpha[m][0][i];
  446. }
  447. Ct[m][0] = 1.0 / w; // scaling coefficient
  448. for(i=0; i<n; i++)
  449. Alpha[m][0][i] = Alpha[m][0][i] * Ct[m][0]; // scaling
  451. // 计算Alpha
  452. for(k=1; k<T[m]; k++) // 第m个训练utterance的各帧
  453. {
  454. Ct[m][k] = 0;
  455. for(i=0; i<n; i++)
  456. {
  457. w = 0;
  458. for(j=0; j<n; j++)
  459. w = w + Alpha[m][k-1][j] * A[j][i];
  460. w = w * B[i][ Out[m][k] ];
  461. Alpha[m][k][i] = w;
  462. Ct[m][k] = Ct[m][k] + w;
  463. }
  464. Ct[m][k] = 1.0 / Ct[m][k];
  465. for(i=0; i<n; i++)
  466. Alpha[m][k][i] = Alpha[m][k][i] * Ct[m][k];
  467. } // end of: for(k=1; k<T[m]; k++)
  469. // 最后一帧状态为吸收态，即只有最后一个状态概率不为0，其余状态概率均为0
  470. if(End_state_mode & ENDSTATE_ABSORB) // End_state_mode是在RECOG.CPP中根据用户输入确定的
  471. {
  472. for(i=0; i<n-1; i++) Alpha[m][ T[m]-1 ][i] = 0;
  473. Belta[m][ T[m]-1 ][n-1] = Ct[m][ T[m]-1 ];
  474. for(i=0; i<n-1; i++) Belta[m][ T[m]-1 ][i] = 0;
  475. }
  476. else
  477. {
  478. for(i=0; i<n; i++) Belta[m][ T[m]-1 ][i] = Ct[m][ T[m]-1 ];
  479. }
  481. // 计算Belta
  482. for(k=T[m]-2; k>=0; k--)
  483. for(i=0; i<n; i++)
  484. {
  485. w = 0;
  486. for(j=0; j<n; j++)
  487. w = w + A[i][j] * B[j][ Out[m][k+1] ] * Belta[m][k+1][j];
  488. Belta[m][k][i] = w * Ct[m][k];
  489. }
  490. }
  491. }
  493. // Re-estimate A[i][j], i = 0 to n-1, j = 0 to n-1.
  495. /*  if((mode == REESTIMATE_A) || (mode == REESTIMATE_ALL))
  496. {
  497. for(i=0; i<n; i++)
  498. for(j=0; j<n; j++)
  499. New_A[i][j] = 0;
  501.   for(m=0; m<Sample_num; m++)
  502.   {
  503.   w = 0;
  504.   for(i=0; i<n; i++)
  505.   w += Alpha[m][T[m]-1][i];
  506.   
  507. for(k=0; k<T[m]-1; k++)
  508. {
  510.   for(i=0; i<n; i++)
  511.   for(j=0; j<n; j++)
  512.   {
  513.   w1 = Alpha[m][k][i] * A[i][j] * B[j][ Out[m][k+1] ] * Belta[m][k+1][j];
  514.   w1 /= w;
  515.   New_A[i][j] += w1;
  516.   }
  517.   }
  518.   }
  519.   
  520. for(i=0; i<n; i++)
  521. {
  522. for(j=0; j<n; j++)
  523. if((j < i) || (j > (i+1)))
  524. New_A[i][j] = 0;
  526.   w = 0;
  527.   for(j=0; j<n; j++) w += New_A[i][j];
  528.   for(j=0; j<n; j++) New_A[i][j] /= w;
  529.   
  530. for(j=0; j<n; j++)
  531. if((j < i) || (j > (i+1)) && (New_A[i][j] < Epsilon_A))
  532. New_A[i][j] = Epsilon_A;
  534.   w = 0;
  535.   for(j=0; j<n; j++) w += New_A[i][j];
  536.   for(j=0; j<n; j++) New_A[i][j] /= w;
  537.   }
  538.   
  539. for(i=0; i<n; i++)
  540. for(j=0; j<n; j++)
  541. A[i][j] = New_A[i][j];
  542.    }*/
  543.    
  544.    // Re-estimate B[i][j], i = 0 to n-1, j = 0 to m-1.
  545.    
  546.    if((mode == REESTIMATE_B) || (mode == REESTIMATE_ALL))
  547.    {
  548.    // clear new B matrix
  549. for(i=0; i<n; i++)
  550. for(j=0; j<Output_num; j++)
  551. New_B[i][j] = 0;
  552.    
  553. for(m=1; m<Sample_num; m++)
  554. {
  555. w = 0;
  556. for(i=0; i<n; i++)
  557. if (T[m] > 0)    // add by wd @ 030331z
  558. w += Alpha[m][T[m]-1][i]; // w: 观察序列概率，即P(O)
  559.    
  560. // reestimate matrix B
  561. for(k=0; k<T[m]; k++)
  562. for(i=0; i<n; i++)
  563. {
  564. New_B[i][Out[m][k]] += Alpha[m][k][i] * Belta[m][k][i] /
  565. Ct[m][k] / w;
  566. }
  567. }
  568.    
  569. for(i=0; i<n; i++)
  570. {
  571. w = 0;
  572. for(j=0; j<Output_num; j++) w += New_B[i][j];
  573. for(j=0; j<Output_num; j++)
  574. {
  575. New_B[i][j] /= w; // 归一化
  576. }
  578. for(j=0; j<Output_num; j++)
  579. if(New_B[i][j] < Epsilon_B) // Epsilon_B=1.0e-6
  580. New_B[i][j] = Epsilon_B;
  581.    
  582. w = 0;
  583. for(j=0; j<Output_num; j++) w += New_B[i][j];
  584. for(j=0; j<Output_num; j++) New_B[i][j] /= w; // 归一化
  586. }
  587.    
  588. for(i=0; i<n; i++) // 更新B矩阵
  589. for(j=0; j<Output_num; j++)
  590. B[i][j] = New_B[i][j];
  592.    }
  593.    
  594.    Residue = 0;
  595.    for(m=0; m<Sample_num; m++)
  596.    if(T[m] > 0) // add by wd @ 030331
  597. Residue += DHMM_Priority(Pi, A, B, Out[m], T[m], n, End_state_mode);
  598.        /*if (T[m] > 0)
  599.    {
  600.    Residue += DHMM_Priority(Pi, A, B, Out[m], T[m], n, End_state_mode);
  601.    }
  602.    else
  603.    Residue += -10; */
  604.    
  605.    for(i=0; i<n; i++) delete New_A[i];
  606.    delete New_A;
  607.    
  608.    for(i=0; i<n; i++) delete New_B[i];
  609.    delete New_B;
  610.    
  611.    for(i=0; i<Sample_num; i++) delete Ct[i];
  612.    delete Ct;
  613.    
  614.    return Residue;
  615. }
  617. static double DHMM_Priority(double *Pi, double **A, double **B, int *Out, int T, int n, int mode)
  619. // input:
  620. //  Pi 初始概率
  621. // A A矩阵
  622. // B B矩阵
  623. // Out 训练集中第m个说话人的观察序列
  624. // T 训练集中第m个说话人的观察序列帧数
  625. // n 状态数
  626. // mode
  628. // return value:
  629. // logP(o1,o2,...,oT|Pi,A,B)
  631. //caller(r_dhmm.cpp): Residue += DHMM_Priority(Pi, A, B, Out[m], T[m], n, End_state_mode);
  633. {
  634. int   i, j, k;
  635. double  w, p;
  636. double  *Alpha1;
  637. double  *Alpha2;
  639. if((Alpha1 = new double [n]) == NULL) { DEBUG_PRINTF("Not enough memory for !n"); ASSERT(0); }
  640. if((Alpha2 = new double [n]) == NULL) { DEBUG_PRINTF("Not enough memory for !n"); ASSERT(0); }
  642. p = 0;
  643. w = 0;
  644. for(i=0; i<n; i++)
  645. {
  646. Alpha1[i] = Pi[i] * B[i][ Out[0] ];
  647. w += Alpha1[i];
  648. }
  649. for(i=0; i<n; i++) Alpha1[i] /= w;
  650. p += log(w);
  652. for(k=1; k<T; k++)
  653. {
  654. for(i=0; i<n; i++)
  655. {
  656. w = 0;
  657. for(j=0; j<n; j++)
  658. w = w + Alpha1[j] * A[j][i];
  659. w = w * B[i][ Out[k] ];
  660. Alpha2[i] = w;
  661. }
  663. if(k < T-1)
  664. {
  665. w = 0;
  666. for(i=0; i<n; i++) w += Alpha2[i];
  667. for(i=0; i<n; i++) Alpha1[i] = Alpha2[i] / w;
  668. p += log(w);
  669. }
  670. }
  672. if(mode & ENDSTATE_ABSORB)
  673. {
  674. p += log(Alpha2[n-1]);
  675. }
  676. else
  677. {
  678. w = 0;
  679. for(i=0; i<n; i++) w += Alpha2[i];
  680. p += log(w);
  681. }
  683. delete Alpha1;
  684. delete Alpha2;
  686. return p;
  687. }
  689. static double DHMM_Calculate(double ***Alpha, double ***Belta, double *Pi, double **A,
  690.  double **B, int Out[140][150], int Sample_num, int Output_num, int *T,
  691.  int n, int mode, int End_state_mode, WORD_SAMPLE * DMAD_pu_Word_Sample)
  692. // 功能：训练某一词条的模型
  694. // input:
  695. // Pi 初始时刻各状态的概率
  696. // Alpha 前向概率
  697. // Belta 后向概率
  698. // A 初始A矩阵
  699. // B 初始B矩阵
  700. // Out 用于训练的各句子的观察序列（输出序列）
  701. // Sample_num 训练集人数
  702. // Output_num VQ码本大小（码字数，即观察值的个数）
  703. // T 用于训练的各句子的帧数
  704. // n 状态数
  705. // mode 训练模式(REESTIMATE_A, REESTIMATE_B, REESTIMATE_ALL)
  706. // End_state_mode 结尾帧状态模式(ENDSTATE_ABSORB = 0x200)
  708. // output:
  709. // A
  710. // B
  712. // caller(r_dhmm.cpp): u = DHMM_Calculate(Alpha, Belta, Pi, A, B, Output, Train_num, Output_num,
  713. //       Frame_num, Status_num, REESTIMATE_ALL, mode);
  714.  
  715. {
  716. int    i, j, k, m;
  717. double w;
  718. double **New_A;
  719. double **New_B;
  720. double **Ct; // to prevent underflow of Alpha and Belta, Ct is the scaling coefficient
  721. double Residue;
  723. if((New_A = new double *[n]) == NULL)
  724. {
  725. DEBUG_PRINTF("Not enough memory for New_A !n");
  726. ASSERT(0);
  727. }
  728. for(i=0; i<n; i++)
  729. {
  730. if((New_A[i] = new double [n]) == NULL)
  731. {
  732. DEBUG_PRINTF("Not enough memory for New_A[%d] !n", i);
  733. ASSERT(0);
  734. }
  735. }
  737. if((New_B = new double *[n]) == NULL)
  738. {
  739. DEBUG_PRINTF("Not enough memory for New_B !n");
  740. ASSERT(0);
  741. }
  742. for(i=0; i<n; i++)
  743. {
  744. if((New_B[i] = new double [Output_num]) == NULL)
  745. {
  746. DEBUG_PRINTF("Not enough memory for New_B[%d] !n", i);
  747. ASSERT(0);
  748. }
  749. }
  751. if((Ct = new double *[Sample_num]) == NULL)
  752. {
  753. DEBUG_PRINTF("Not enough memory for Ct !n");
  754. ASSERT(0);
  755. }
  756. for(i=0; i<Sample_num; i++)
  757. if((Ct[i] = new double [T[i]]) == NULL)
  758. {
  759. DEBUG_PRINTF("Not enough memory for Ct[%d] !n", i);
  760. ASSERT(0);
  761. }
  763. // Calculate Alpha & Belta.
  765. for(m=0; m<Sample_num; m++)
  766. {
  767. w = 0;
  768. for(i=0; i<n; i++)
  769. {
  770. Alpha[m][0][i] = Pi[i] * B[i][ DMAD_pu_Word_Sample[m].pn_VQed_Feature_Sequence[i] ];
  771. w = w + Alpha[m][0][i];
  772. }
  773. Ct[m][0] = 1.0 / w; // scaling coefficient
  774. for(i=0; i<n; i++)
  775. Alpha[m][0][i] = Alpha[m][0][i] * Ct[m][0]; // scaling
  777. // 计算Alpha
  778. for(k=1; k<T[m]; k++) // 第m个训练utterance的各帧
  779. {
  780. Ct[m][k] = 0;
  781. for(i=0; i<n; i++)
  782. {
  783. w = 0;
  784. for(j=0; j<n; j++)
  785. w = w + Alpha[m][k-1][j] * A[j][i];
  786. w = w * B[i][ Out[m][k] ];
  787. Alpha[m][k][i] = w;
  788. Ct[m][k] = Ct[m][k] + w;
  789. }
  790. Ct[m][k] = 1.0 / Ct[m][k];
  791. for(i=0; i<n; i++)
  792. Alpha[m][k][i] = Alpha[m][k][i] * Ct[m][k];
  793. } // end of: for(k=1; k<T[m]; k++)
  795. // 最后一帧状态为吸收态，即只有最后一个状态概率不为0，其余状态概率均为0
  796. if(End_state_mode & ENDSTATE_ABSORB) // End_state_mode是在RECOG.CPP中根据用户输入确定的
  797. {
  798. for(i=0; i<n-1; i++) Alpha[m][ T[m]-1 ][i] = 0;
  799. Belta[m][ T[m]-1 ][n-1] = Ct[m][ T[m]-1 ];
  800. for(i=0; i<n-1; i++) Belta[m][ T[m]-1 ][i] = 0;
  801. }
  802. else
  803. {
  804. for(i=0; i<n; i++) Belta[m][ T[m]-1 ][i] = Ct[m][ T[m]-1 ];
  805. }
  807. // 计算Belta
  808. for(k=T[m]-2; k>=0; k--)
  809. for(i=0; i<n; i++)
  810. {
  811. w = 0;
  812. for(j=0; j<n; j++)
  813. w = w + A[i][j] * B[j][ Out[m][k+1] ] * Belta[m][k+1][j];
  814. Belta[m][k][i] = w * Ct[m][k];
  815. }
  816. }
  818. // Re-estimate A[i][j], i = 0 to n-1, j = 0 to n-1.
  820. /*  if((mode == REESTIMATE_A) || (mode == REESTIMATE_ALL))
  821. {
  822. for(i=0; i<n; i++)
  823. for(j=0; j<n; j++)
  824. New_A[i][j] = 0;
  826.   for(m=0; m<Sample_num; m++)
  827.   {
  828.   w = 0;
  829.   for(i=0; i<n; i++)
  830.   w += Alpha[m][T[m]-1][i];
  831.   
  832. for(k=0; k<T[m]-1; k++)
  833. {
  835.   for(i=0; i<n; i++)
  836.   for(j=0; j<n; j++)
  837.   {
  838.   w1 = Alpha[m][k][i] * A[i][j] * B[j][ Out[m][k+1] ] * Belta[m][k+1][j];
  839.   w1 /= w;
  840.   New_A[i][j] += w1;
  841.   }
  842.   }
  843.   }
  844.   
  845. for(i=0; i<n; i++)
  846. {
  847. for(j=0; j<n; j++)
  848. if((j < i) || (j > (i+1)))
  849. New_A[i][j] = 0;
  851.   w = 0;
  852.   for(j=0; j<n; j++) w += New_A[i][j];
  853.   for(j=0; j<n; j++) New_A[i][j] /= w;
  854.   
  855. for(j=0; j<n; j++)
  856. if((j < i) || (j > (i+1)) && (New_A[i][j] < Epsilon_A))
  857. New_A[i][j] = Epsilon_A;
  859.   w = 0;
  860.   for(j=0; j<n; j++) w += New_A[i][j];
  861.   for(j=0; j<n; j++) New_A[i][j] /= w;
  862.   }
  863.   
  864. for(i=0; i<n; i++)
  865. for(j=0; j<n; j++)
  866. A[i][j] = New_A[i][j];
  867.    }*/
  868.    
  869.    // Re-estimate B[i][j], i = 0 to n-1, j = 0 to m-1.
  870.    
  871.    if((mode == REESTIMATE_B) || (mode == REESTIMATE_ALL))
  872.    {
  873.    // clear new B matrix
  874.    for(i=0; i<n; i++)
  875.    for(j=0; j<Output_num; j++)
  876.    New_B[i][j] = 0;
  877.    
  878.    for(m=1; m<Sample_num; m++)
  879.    {
  880.    w = 0;
  881.    for(i=0; i<n; i++)
  882.    w += Alpha[m][T[m]-1][i]; // w: 观察序列概率，即P(O)
  883.    
  884.    // reestimate matrix B
  885.    for(k=0; k<T[m]; k++)
  886.    for(i=0; i<n; i++)
  887.    {
  888.    New_B[i][Out[m][k]] += Alpha[m][k][i] * Belta[m][k][i] /
  889.    Ct[m][k] / w;
  890.    }
  891.    }
  892.    
  893.    for(i=0; i<n; i++)
  894.    {
  895.    w = 0;
  896.    for(j=0; j<Output_num; j++) w += New_B[i][j];
  897.    for(j=0; j<Output_num; j++)
  898.    {
  899.    New_B[i][j] /= w; // 归一化
  900.    }
  902.    for(j=0; j<Output_num; j++)
  903.    if(New_B[i][j] < Epsilon_B) // Epsilon_B=1.0e-6
  904.    New_B[i][j] = Epsilon_B;
  905.    
  906.    w = 0;
  907.    for(j=0; j<Output_num; j++) w += New_B[i][j];
  908.    for(j=0; j<Output_num; j++) New_B[i][j] /= w; // 归一化
  910.    }
  911.    
  912.    for(i=0; i<n; i++) // 更新B矩阵
  913.    for(j=0; j<Output_num; j++)
  914.    B[i][j] = New_B[i][j];
  915.    }
  916.    
  917.    Residue = 0;
  918.    for(m=1; m<Sample_num; m++)
  919.    Residue += DHMM_Priority(Pi, A, B, Out[m], T[m], n, End_state_mode);
  920.    
  921.    for(i=0; i<n; i++) delete New_A[i];
  922.    delete New_A;
  923.    
  924.    for(i=0; i<n; i++) delete New_B[i];
  925.    delete New_B;
  926.    
  927.    for(i=0; i<Sample_num; i++) delete Ct[i];
  928.    delete Ct;
  929.    
  930.    return Residue;
  931. }