src/libsphinxbase/feat/feat.c

   1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
   2 /* ====================================================================
   3  * Copyright (c) 1999-2004 Carnegie Mellon University.  All rights
   4  * reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions
   8  * are met:
   9  *
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  *
  13  * 2. Redistributions in binary form must reproduce the above copyright
  14  *    notice, this list of conditions and the following disclaimer in
  15  *    the documentation and/or other materials provided with the
  16  *    distribution.
  17  *
  18  * This work was supported in part by funding from the Defense Advanced
  19  * Research Projects Agency and the National Science Foundation of the
  20  * United States of America, and the CMU Sphinx Speech Consortium.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
  23  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
  24  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  25  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
  26  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  27  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  28  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  29  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  30  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  31  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  32  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  33  *
  34  * ====================================================================
  35  *
  36  */
  37 /*
  38  * feat.c -- Feature vector description and cepstra->feature computation.
  39  *
  40  * **********************************************
  41  * CMU ARPA Speech Project
  42  *
  43  * Copyright (c) 1996 Carnegie Mellon University.
  44  * ALL RIGHTS RESERVED.
  45  * **********************************************
  46  *
  47  * HISTORY
  48  * $Log$
  49  * Revision 1.22  2006/02/23  03:59:40  arthchan2003
  50  * Merged from branch SPHINX3_5_2_RCI_IRII_BRANCH: a, Free buffers correctly. b, Fixed dox-doc.
  51  *
  52  * Revision 1.21.4.3  2005/10/17 04:45:57  arthchan2003
  53  * Free stuffs in cmn and feat corectly.
  54  *
  55  * Revision 1.21.4.2  2005/09/26 02:19:57  arthchan2003
  56  * Add message to show the directory which the feature is searched for.
  57  *
  58  * Revision 1.21.4.1  2005/07/03 22:55:50  arthchan2003
  59  * More correct deallocation in feat.c. The cmn deallocation is still not correct at this point.
  60  *
  61  * Revision 1.21  2005/06/22 03:29:35  arthchan2003
  62  * Makefile.am s  for all subdirectory of libs3decoder/
  63  *
  64  * Revision 1.4  2005/04/21 23:50:26  archan
  65  * Some more refactoring on the how reporting of structures inside kbcore_t is done, it is now 50% nice. Also added class-based LM test case into test-decode.sh.in.  At this moment, everything in search mode 5 is already done.  It is time to test the idea whether the search can really be used.
  66  *
  67  * Revision 1.3  2005/03/30 01:22:46  archan
  68  * Fixed mistakes in last updates. Add
  69  *
  70  *
  71  * 20.Apr.2001  RAH (rhoughton@mediasite.com, ricky.houghton@cs.cmu.edu)
  72  *              Adding feat_free() to free allocated memory
  73  *
  74  * 02-Jan-2001  Rita Singh (rsingh@cs.cmu.edu) at Carnegie Mellon University
  75  *              Modified feat_s2mfc2feat_block() to handle empty buffers at
  76  *              the end of an utterance
  77  *
  78  * 30-Dec-2000  Rita Singh (rsingh@cs.cmu.edu) at Carnegie Mellon University
  79  *              Added feat_s2mfc2feat_block() to allow feature computation
  80  *              from sequences of blocks of cepstral vectors
  81  *
  82  * 12-Jun-98    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
  83  *              Major changes to accommodate arbitrary feature input types.  Added
  84  *              feat_read(), moved various cep2feat functions from other files into
  85  *              this one.  Also, made this module object-oriented with the feat_t type.
  86  *              Changed definition of s2mfc_read to let the caller manage MFC buffers.
  87  *
  88  * 03-Oct-96    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
  89  *              Added unistd.h include.
  90  *
  91  * 02-Oct-96    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
  92  *              Added check for sf argument to s2mfc_read being within file size.
  93  *
  94  * 18-Sep-96    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
  95  *              Added sf, ef parameters to s2mfc_read().
  96  *
  97  * 10-Jan-96    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
  98  *              Added feat_cepsize().
  99  *              Added different feature-handling (s2_4x, s3_1x39 at this point).
 100  *              Moved feature-dependent functions to feature-dependent files.
 101  *
 102  * 09-Jan-96    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 103  *              Moved constant declarations from feat.h into here.
 104  *
 105  * 04-Nov-95    M K Ravishankar (rkm@cs.cmu.edu) at Carnegie Mellon University
 106  *              Created.
 107  */
 108
 109
 110 /*
 111  * This module encapsulates different feature streams used by the Sphinx group.  New
 112  * stream types can be added by augmenting feat_init() and providing an accompanying
 113  * compute_feat function.  It also provides a "generic" feature vector definition for
 114  * handling "arbitrary" speech input feature types (see the last section in feat_init()).
 115  * In this case the speech input data should already be feature vectors; no computation,
 116  * such as MFC->feature conversion, is available or needed.
 117  */
 118
 119 #include <assert.h>
 120 #include <string.h>
 121 #ifdef HAVE_CONFIG_H
 122 #include <config.h>
 123 #endif
 124
 125 #ifdef _MSC_VER
 126 #pragma warning (disable: 4244 4996)
 127 #endif
 128
 129 #include "sphinxbase/fe.h"
 130 #include "sphinxbase/feat.h"
 131 #include "sphinxbase/bio.h"
 132 #include "sphinxbase/pio.h"
 133 #include "sphinxbase/cmn.h"
 134 #include "sphinxbase/agc.h"
 135 #include "sphinxbase/err.h"
 136 #include "sphinxbase/ckd_alloc.h"
 137 #include "sphinxbase/prim_type.h"
 138 #include "sphinxbase/glist.h"
 139
 140 #define FEAT_VERSION    "1.0"
 141 #define FEAT_DCEP_WIN           2
 142
 143 #ifdef DUMP_FEATURES
 144 static void
 145 cep_dump_dbg(feat_t *fcb, mfcc_t **mfc, int32 nfr, const char *text)
 146 {
 147     int32 i, j;
 148
 149     E_INFO("%s\n", text);
 150     for (i = 0; i < nfr; i++) {
 151         for (j = 0; j < fcb->cepsize; j++) {
 152             fprintf(stderr, "%f ", MFCC2FLOAT(mfc[i][j]));
 153         }
 154         fprintf(stderr, "\n");
 155     }
 156 }
 157 static void
 158 feat_print_dbg(feat_t *fcb, mfcc_t ***feat, int32 nfr, const char *text)
 159 {
 160     E_INFO("%s\n", text);
 161     feat_print(fcb, feat, nfr, stderr);
 162 }
 163 #else /* !DUMP_FEATURES */
 164 #define cep_dump_dbg(fcb,mfc,nfr,text)
 165 #define feat_print_dbg(fcb,mfc,nfr,text)
 166 #endif
 167
 168 int32 **
 169 parse_subvecs(char const *str)
 170 {
 171     char const *strp;
 172     int32 n, n2, l;
 173     glist_t dimlist;            /* List of dimensions in one subvector */
 174     glist_t veclist;            /* List of dimlists (subvectors) */
 175     int32 **subvec;
 176     gnode_t *gn, *gn2;
 177
 178     veclist = NULL;
 179
 180     strp = str;
 181     for (;;) {
 182         dimlist = NULL;
 183
 184         for (;;) {
 185             if (sscanf(strp, "%d%n", &n, &l) != 1)
 186                 E_FATAL("'%s': Couldn't read int32 @pos %d\n", str,
 187                         strp - str);
 188             strp += l;
 189
 190             if (*strp == '-') {
 191                 strp++;
 192
 193                 if (sscanf(strp, "%d%n", &n2, &l) != 1)
 194                     E_FATAL("'%s': Couldn't read int32 @pos %d\n", str,
 195                             strp - str);
 196                 strp += l;
 197             }
 198             else
 199                 n2 = n;
 200
 201             if ((n < 0) || (n > n2))
 202                 E_FATAL("'%s': Bad subrange spec ending @pos %d\n", str,
 203                         strp - str);
 204
 205             for (; n <= n2; n++) {
 206                 gnode_t *gn;
 207                 for (gn = dimlist; gn; gn = gnode_next(gn))
 208                     if (gnode_int32(gn) == n)
 209                         break;
 210                 if (gn != NULL)
 211                     E_FATAL("'%s': Duplicate dimension ending @pos %d\n",
 212                             str, strp - str);
 213
 214                 dimlist = glist_add_int32(dimlist, n);
 215             }
 216
 217             if ((*strp == '\0') || (*strp == '/'))
 218                 break;
 219
 220             if (*strp != ',')
 221                 E_FATAL("'%s': Bad delimiter @pos %d\n", str, strp - str);
 222
 223             strp++;
 224         }
 225
 226         veclist = glist_add_ptr(veclist, (void *) dimlist);
 227
 228         if (*strp == '\0')
 229             break;
 230
 231         assert(*strp == '/');
 232         strp++;
 233     }
 234
 235     /* Convert the glists to arrays; remember the glists are in reverse order of the input! */
 236     n = glist_count(veclist);   /* #Subvectors */
 237     subvec = (int32 **) ckd_calloc(n + 1, sizeof(int32 *));     /* +1 for sentinel */
 238     subvec[n] = NULL;           /* sentinel */
 239
 240     for (--n, gn = veclist; (n >= 0) && gn; gn = gnode_next(gn), --n) {
 241         gn2 = (glist_t) gnode_ptr(gn);
 242
 243         n2 = glist_count(gn2);  /* Length of this subvector */
 244         if (n2 <= 0)
 245             E_FATAL("'%s': 0-length subvector\n", str);
 246
 247         subvec[n] = (int32 *) ckd_calloc(n2 + 1, sizeof(int32));        /* +1 for sentinel */
 248         subvec[n][n2] = -1;     /* sentinel */
 249
 250         for (--n2; (n2 >= 0) && gn2; gn2 = gnode_next(gn2), --n2)
 251             subvec[n][n2] = gnode_int32(gn2);
 252         assert((n2 < 0) && (!gn2));
 253     }
 254     assert((n < 0) && (!gn));
 255
 256     /* Free the glists */
 257     for (gn = veclist; gn; gn = gnode_next(gn)) {
 258         gn2 = (glist_t) gnode_ptr(gn);
 259         glist_free(gn2);
 260     }
 261     glist_free(veclist);
 262
 263     return subvec;
 264 }
 265
 266 void
 267 subvecs_free(int32 **subvecs)
 268 {
 269     int32 **sv;
 270
 271     for (sv = subvecs; sv && *sv; ++sv)
 272         ckd_free(*sv);
 273     ckd_free(subvecs);
 274 }
 275
 276 int
 277 feat_set_subvecs(feat_t *fcb, int32 **subvecs)
 278 {
 279     int32 **sv;
 280     int32 n_sv, n_dim, i;
 281
 282     if (subvecs == NULL) {
 283         subvecs_free(fcb->subvecs);
 284         ckd_free(fcb->sv_buf);
 285         ckd_free(fcb->sv_len);
 286         fcb->n_sv = 0;
 287         fcb->subvecs = NULL;
 288         fcb->sv_len = NULL;
 289         fcb->sv_buf = NULL;
 290         fcb->sv_dim = 0;
 291         return 0;
 292     }
 293
 294     if (fcb->n_stream != 1) {
 295         E_ERROR("Subvector specifications require single-stream features!");
 296         return -1;
 297     }
 298
 299     n_sv = 0;
 300     n_dim = 0;
 301     for (sv = subvecs; sv && *sv; ++sv) {
 302         int32 *d;
 303
 304         for (d = *sv; d && *d != -1; ++d) {
 305             ++n_dim;
 306         }
 307         ++n_sv;
 308     }
 309     if (n_dim > feat_dimension(fcb)) {
 310         E_ERROR("Total dimensionality of subvector specification %d "
 311                 "> feature dimensionality %d\n", n_dim, feat_dimension(fcb));
 312         return -1;
 313     }
 314
 315     fcb->n_sv = n_sv;
 316     fcb->subvecs = subvecs;
 317     fcb->sv_len = ckd_calloc(n_sv, sizeof(*fcb->sv_len));
 318     fcb->sv_buf = ckd_calloc(n_dim, sizeof(*fcb->sv_buf));
 319     fcb->sv_dim = n_dim;
 320     for (i = 0; i < n_sv; ++i) {
 321         int32 *d;
 322         for (d = subvecs[i]; d && *d != -1; ++d) {
 323             ++fcb->sv_len[i];
 324         }
 325     }
 326
 327     return 0;
 328 }
 329
 330 /**
 331  * Project feature components to subvectors (if any).
 332  */
 333 static void
 334 feat_subvec_project(feat_t *fcb, mfcc_t ***inout_feat, uint32 nfr)
 335 {
 336     uint32 i;
 337
 338     if (fcb->subvecs == NULL)
 339         return;
 340     for (i = 0; i < nfr; ++i) {
 341         mfcc_t *out;
 342         int32 j;
 343
 344         out = fcb->sv_buf;
 345         for (j = 0; j < fcb->n_sv; ++j) {
 346             int32 *d;
 347             for (d = fcb->subvecs[j]; d && *d != -1; ++d) {
 348                 *out++ = inout_feat[i][0][*d];
 349             }
 350         }
 351         memcpy(inout_feat[i][0], fcb->sv_buf, fcb->sv_dim * sizeof(*fcb->sv_buf));
 352     }
 353 }
 354
 355 mfcc_t ***
 356 feat_array_alloc(feat_t * fcb, int32 nfr)
 357 {
 358     int32 i, j, k;
 359     mfcc_t *data, *d, ***feat;
 360
 361     assert(fcb);
 362     assert(nfr > 0);
 363     assert(feat_dimension(fcb) > 0);
 364
 365     /* Make sure to use the dimensionality of the features *before*
 366        LDA and subvector projection. */
 367     k = 0;
 368     for (i = 0; i < fcb->n_stream; ++i)
 369         k += fcb->stream_len[i];
 370     assert(k >= feat_dimension(fcb));
 371     assert(k >= fcb->sv_dim);
 372
 373     feat =
 374         (mfcc_t ***) ckd_calloc_2d(nfr, feat_dimension1(fcb), sizeof(mfcc_t *));
 375     data = (mfcc_t *) ckd_calloc(nfr * k, sizeof(mfcc_t));
 376
 377     for (i = 0; i < nfr; i++) {
 378         d = data + i * k;
 379         for (j = 0; j < feat_dimension1(fcb); j++) {
 380             feat[i][j] = d;
 381             d += feat_dimension2(fcb, j);
 382         }
 383     }
 384
 385     return feat;
 386 }
 387
 388 void
 389 feat_array_free(mfcc_t ***feat)
 390 {
 391     ckd_free(feat[0][0]);
 392     ckd_free_2d((void **)feat);
 393 }
 394
 395 static void
 396 feat_s2_4x_cep2feat(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 397 {
 398     mfcc_t *f;
 399     mfcc_t *w, *_w;
 400     mfcc_t *w1, *w_1, *_w1, *_w_1;
 401     mfcc_t d1, d2;
 402     int32 i, j;
 403
 404     assert(fcb);
 405     assert(feat_cepsize(fcb) == 13);
 406     assert(feat_n_stream(fcb) == 4);
 407     assert(feat_stream_len(fcb, 0) == 12);
 408     assert(feat_stream_len(fcb, 1) == 24);
 409     assert(feat_stream_len(fcb, 2) == 3);
 410     assert(feat_stream_len(fcb, 3) == 12);
 411     assert(feat_window_size(fcb) == 4);
 412
 413     /* CEP; skip C0 */
 414     memcpy(feat[0], mfc[0] + 1, (feat_cepsize(fcb) - 1) * sizeof(mfcc_t));
 415
 416     /*
 417      * DCEP(SHORT): mfc[2] - mfc[-2]
 418      * DCEP(LONG):  mfc[4] - mfc[-4]
 419      */
 420     w = mfc[2] + 1;             /* +1 to skip C0 */
 421     _w = mfc[-2] + 1;
 422
 423     f = feat[1];
 424     for (i = 0; i < feat_cepsize(fcb) - 1; i++) /* Short-term */
 425         f[i] = w[i] - _w[i];
 426
 427     w = mfc[4] + 1;             /* +1 to skip C0 */
 428     _w = mfc[-4] + 1;
 429
 430     for (j = 0; j < feat_cepsize(fcb) - 1; i++, j++)    /* Long-term */
 431         f[i] = w[j] - _w[j];
 432
 433     /* D2CEP: (mfc[3] - mfc[-1]) - (mfc[1] - mfc[-3]) */
 434     w1 = mfc[3] + 1;            /* Final +1 to skip C0 */
 435     _w1 = mfc[-1] + 1;
 436     w_1 = mfc[1] + 1;
 437     _w_1 = mfc[-3] + 1;
 438
 439     f = feat[3];
 440     for (i = 0; i < feat_cepsize(fcb) - 1; i++) {
 441         d1 = w1[i] - _w1[i];
 442         d2 = w_1[i] - _w_1[i];
 443
 444         f[i] = d1 - d2;
 445     }
 446
 447     /* POW: C0, DC0, D2C0; differences computed as above for rest of cep */
 448     f = feat[2];
 449     f[0] = mfc[0][0];
 450     f[1] = mfc[2][0] - mfc[-2][0];
 451
 452     d1 = mfc[3][0] - mfc[-1][0];
 453     d2 = mfc[1][0] - mfc[-3][0];
 454     f[2] = d1 - d2;
 455 }
 456
 457
 458 static void
 459 feat_s3_1x39_cep2feat(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 460 {
 461     mfcc_t *f;
 462     mfcc_t *w, *_w;
 463     mfcc_t *w1, *w_1, *_w1, *_w_1;
 464     mfcc_t d1, d2;
 465     int32 i;
 466
 467     assert(fcb);
 468     assert(feat_cepsize(fcb) == 13);
 469     assert(feat_n_stream(fcb) == 1);
 470     assert(feat_stream_len(fcb, 0) == 39);
 471     assert(feat_window_size(fcb) == 3);
 472
 473     /* CEP; skip C0 */
 474     memcpy(feat[0], mfc[0] + 1, (feat_cepsize(fcb) - 1) * sizeof(mfcc_t));
 475     /*
 476      * DCEP: mfc[2] - mfc[-2];
 477      */
 478     f = feat[0] + feat_cepsize(fcb) - 1;
 479     w = mfc[2] + 1;             /* +1 to skip C0 */
 480     _w = mfc[-2] + 1;
 481
 482     for (i = 0; i < feat_cepsize(fcb) - 1; i++)
 483         f[i] = w[i] - _w[i];
 484
 485     /* POW: C0, DC0, D2C0 */
 486     f += feat_cepsize(fcb) - 1;
 487
 488     f[0] = mfc[0][0];
 489     f[1] = mfc[2][0] - mfc[-2][0];
 490
 491     d1 = mfc[3][0] - mfc[-1][0];
 492     d2 = mfc[1][0] - mfc[-3][0];
 493     f[2] = d1 - d2;
 494
 495     /* D2CEP: (mfc[3] - mfc[-1]) - (mfc[1] - mfc[-3]) */
 496     f += 3;
 497
 498     w1 = mfc[3] + 1;            /* Final +1 to skip C0 */
 499     _w1 = mfc[-1] + 1;
 500     w_1 = mfc[1] + 1;
 501     _w_1 = mfc[-3] + 1;
 502
 503     for (i = 0; i < feat_cepsize(fcb) - 1; i++) {
 504         d1 = w1[i] - _w1[i];
 505         d2 = w_1[i] - _w_1[i];
 506
 507         f[i] = d1 - d2;
 508     }
 509 }
 510
 511
 512 static void
 513 feat_s3_cep(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 514 {
 515     assert(fcb);
 516     assert(feat_n_stream(fcb) == 1);
 517     assert(feat_window_size(fcb) == 0);
 518
 519     /* CEP */
 520     memcpy(feat[0], mfc[0], feat_cepsize(fcb) * sizeof(mfcc_t));
 521 }
 522
 523 static void
 524 feat_s3_cep_dcep(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 525 {
 526     mfcc_t *f;
 527     mfcc_t *w, *_w;
 528     int32 i;
 529
 530     assert(fcb);
 531     assert(feat_n_stream(fcb) == 1);
 532     assert(feat_stream_len(fcb, 0) == feat_cepsize(fcb) * 2);
 533     assert(feat_window_size(fcb) == 2);
 534
 535     /* CEP */
 536     memcpy(feat[0], mfc[0], feat_cepsize(fcb) * sizeof(mfcc_t));
 537
 538     /*
 539      * DCEP: mfc[2] - mfc[-2];
 540      */
 541     f = feat[0] + feat_cepsize(fcb);
 542     w = mfc[2];
 543     _w = mfc[-2];
 544
 545     for (i = 0; i < feat_cepsize(fcb); i++)
 546         f[i] = w[i] - _w[i];
 547 }
 548
 549 static void
 550 feat_1s_c_d_dd_cep2feat(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 551 {
 552     mfcc_t *f;
 553     mfcc_t *w, *_w;
 554     mfcc_t *w1, *w_1, *_w1, *_w_1;
 555     mfcc_t d1, d2;
 556     int32 i;
 557
 558     assert(fcb);
 559     assert(feat_n_stream(fcb) == 1);
 560     assert(feat_stream_len(fcb, 0) == feat_cepsize(fcb) * 3);
 561     assert(feat_window_size(fcb) == FEAT_DCEP_WIN + 1);
 562
 563     /* CEP */
 564     memcpy(feat[0], mfc[0], feat_cepsize(fcb) * sizeof(mfcc_t));
 565
 566     /*
 567      * DCEP: mfc[w] - mfc[-w], where w = FEAT_DCEP_WIN;
 568      */
 569     f = feat[0] + feat_cepsize(fcb);
 570     w = mfc[FEAT_DCEP_WIN];
 571     _w = mfc[-FEAT_DCEP_WIN];
 572
 573     for (i = 0; i < feat_cepsize(fcb); i++)
 574         f[i] = w[i] - _w[i];
 575
 576     /*
 577      * D2CEP: (mfc[w+1] - mfc[-w+1]) - (mfc[w-1] - mfc[-w-1]),
 578      * where w = FEAT_DCEP_WIN
 579      */
 580     f += feat_cepsize(fcb);
 581
 582     w1 = mfc[FEAT_DCEP_WIN + 1];
 583     _w1 = mfc[-FEAT_DCEP_WIN + 1];
 584     w_1 = mfc[FEAT_DCEP_WIN - 1];
 585     _w_1 = mfc[-FEAT_DCEP_WIN - 1];
 586
 587     for (i = 0; i < feat_cepsize(fcb); i++) {
 588         d1 = w1[i] - _w1[i];
 589         d2 = w_1[i] - _w_1[i];
 590
 591         f[i] = d1 - d2;
 592     }
 593 }
 594
 595 static void
 596 feat_1s_c_d_ld_dd_cep2feat(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 597 {
 598     mfcc_t *f;
 599     mfcc_t *w, *_w;
 600     mfcc_t *w1, *w_1, *_w1, *_w_1;
 601     mfcc_t d1, d2;
 602     int32 i;
 603
 604     assert(fcb);
 605     assert(feat_n_stream(fcb) == 1);
 606     assert(feat_stream_len(fcb, 0) == feat_cepsize(fcb) * 4);
 607     assert(feat_window_size(fcb) == FEAT_DCEP_WIN * 2);
 608
 609     /* CEP */
 610     memcpy(feat[0], mfc[0], feat_cepsize(fcb) * sizeof(mfcc_t));
 611
 612     /*
 613      * DCEP: mfc[w] - mfc[-w], where w = FEAT_DCEP_WIN;
 614      */
 615     f = feat[0] + feat_cepsize(fcb);
 616     w = mfc[FEAT_DCEP_WIN];
 617     _w = mfc[-FEAT_DCEP_WIN];
 618
 619     for (i = 0; i < feat_cepsize(fcb); i++)
 620         f[i] = w[i] - _w[i];
 621
 622     /*
 623      * LDCEP: mfc[w] - mfc[-w], where w = FEAT_DCEP_WIN * 2;
 624      */
 625     f += feat_cepsize(fcb);
 626     w = mfc[FEAT_DCEP_WIN * 2];
 627     _w = mfc[-FEAT_DCEP_WIN * 2];
 628
 629     for (i = 0; i < feat_cepsize(fcb); i++)
 630         f[i] = w[i] - _w[i];
 631
 632     /*
 633      * D2CEP: (mfc[w+1] - mfc[-w+1]) - (mfc[w-1] - mfc[-w-1]),
 634      * where w = FEAT_DCEP_WIN
 635      */
 636     f += feat_cepsize(fcb);
 637
 638     w1 = mfc[FEAT_DCEP_WIN + 1];
 639     _w1 = mfc[-FEAT_DCEP_WIN + 1];
 640     w_1 = mfc[FEAT_DCEP_WIN - 1];
 641     _w_1 = mfc[-FEAT_DCEP_WIN - 1];
 642
 643     for (i = 0; i < feat_cepsize(fcb); i++) {
 644         d1 = w1[i] - _w1[i];
 645         d2 = w_1[i] - _w_1[i];
 646
 647         f[i] = d1 - d2;
 648     }
 649 }
 650
 651 static void
 652 feat_copy(feat_t * fcb, mfcc_t ** mfc, mfcc_t ** feat)
 653 {
 654     int32 win, i, j;
 655
 656     win = feat_window_size(fcb);
 657
 658     /* Concatenate input features */
 659     for (i = -win; i <= win; ++i) {
 660         uint32 spos = 0;
 661
 662         for (j = 0; j < feat_n_stream(fcb); ++j) {
 663             uint32 stream_len;
 664
 665             /* Unscale the stream length by the window. */
 666             stream_len = feat_stream_len(fcb, j) / (2 * win + 1);
 667             memcpy(feat[j] + ((i + win) * stream_len),
 668                    mfc[i] + spos,
 669                    stream_len * sizeof(mfcc_t));
 670             spos += stream_len;
 671         }
 672     }
 673 }
 674
 675 feat_t *
 676 feat_init(char const *type, cmn_type_t cmn, int32 varnorm,
 677           agc_type_t agc, int32 breport, int32 cepsize)
 678 {
 679     feat_t *fcb;
 680
 681     if (cepsize == 0)
 682         cepsize = 13;
 683     if (breport)
 684         E_INFO
 685             ("Initializing feature stream to type: '%s', ceplen=%d, CMN='%s', VARNORM='%s', AGC='%s'\n",
 686              type, cepsize, cmn_type_str[cmn], varnorm ? "yes" : "no", agc_type_str[agc]);
 687
 688     fcb = (feat_t *) ckd_calloc(1, sizeof(feat_t));
 689     fcb->refcount = 1;
 690     fcb->name = (char *) ckd_salloc(type);
 691     if (strcmp(type, "s2_4x") == 0) {
 692         /* Sphinx-II format 4-stream feature (Hack!! hardwired constants below) */
 693         if (cepsize != 13) {
 694             E_ERROR("s2_4x features require cepsize == 13\n");
 695             ckd_free(fcb);
 696             return NULL;
 697         }
 698         fcb->cepsize = 13;
 699         fcb->n_stream = 4;
 700         fcb->stream_len = (int32 *) ckd_calloc(4, sizeof(int32));
 701         fcb->stream_len[0] = 12;
 702         fcb->stream_len[1] = 24;
 703         fcb->stream_len[2] = 3;
 704         fcb->stream_len[3] = 12;
 705         fcb->out_dim = 51;
 706         fcb->window_size = 4;
 707         fcb->compute_feat = feat_s2_4x_cep2feat;
 708     }
 709     else if ((strcmp(type, "s3_1x39") == 0) || (strcmp(type, "1s_12c_12d_3p_12dd") == 0)) {
 710         /* 1-stream cep/dcep/pow/ddcep (Hack!! hardwired constants below) */
 711         if (cepsize != 13) {
 712             E_ERROR("s2_4x features require cepsize == 13\n");
 713             ckd_free(fcb);
 714             return NULL;
 715         }
 716         fcb->cepsize = 13;
 717         fcb->n_stream = 1;
 718         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 719         fcb->stream_len[0] = 39;
 720         fcb->out_dim = 39;
 721         fcb->window_size = 3;
 722         fcb->compute_feat = feat_s3_1x39_cep2feat;
 723     }
 724     else if (strncmp(type, "1s_c_d_dd", 9) == 0) {
 725         fcb->cepsize = cepsize;
 726         fcb->n_stream = 1;
 727         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 728         fcb->stream_len[0] = cepsize * 3;
 729         fcb->out_dim = cepsize * 3;
 730         fcb->window_size = FEAT_DCEP_WIN + 1; /* ddcep needs the extra 1 */
 731         fcb->compute_feat = feat_1s_c_d_dd_cep2feat;
 732     }
 733     else if (strncmp(type, "1s_c_d_ld_dd", 12) == 0) {
 734         fcb->cepsize = cepsize;
 735         fcb->n_stream = 1;
 736         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 737         fcb->stream_len[0] = cepsize * 4;
 738         fcb->out_dim = cepsize * 4;
 739         fcb->window_size = FEAT_DCEP_WIN * 2;
 740         fcb->compute_feat = feat_1s_c_d_ld_dd_cep2feat;
 741     }
 742     else if (strncmp(type, "cep_dcep", 8) == 0 || strncmp(type, "1s_c_d", 6) == 0) {
 743         /* 1-stream cep/dcep */
 744         fcb->cepsize = cepsize;
 745         fcb->n_stream = 1;
 746         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 747         fcb->stream_len[0] = feat_cepsize(fcb) * 2;
 748         fcb->out_dim = fcb->stream_len[0];
 749         fcb->window_size = 2;
 750         fcb->compute_feat = feat_s3_cep_dcep;
 751     }
 752     else if (strncmp(type, "cep", 3) == 0 || strncmp(type, "1s_c", 4) == 0) {
 753         /* 1-stream cep */
 754         fcb->cepsize = cepsize;
 755         fcb->n_stream = 1;
 756         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 757         fcb->stream_len[0] = feat_cepsize(fcb);
 758         fcb->out_dim = fcb->stream_len[0];
 759         fcb->window_size = 0;
 760         fcb->compute_feat = feat_s3_cep;
 761     }
 762     else if (strncmp(type, "1s_3c", 5) == 0 || strncmp(type, "1s_4c", 5) == 0) {
 763         /* 1-stream cep with frames concatenated, so called cepwin features */
 764         if (strncmp(type, "1s_3c", 5) == 0)
 765             fcb->window_size = 3;
 766         else
 767             fcb->window_size = 4;
 768
 769         fcb->cepsize = cepsize;
 770         fcb->n_stream = 1;
 771         fcb->stream_len = (int32 *) ckd_calloc(1, sizeof(int32));
 772         fcb->stream_len[0] = feat_cepsize(fcb) * (2 * fcb->window_size + 1);
 773         fcb->out_dim = fcb->stream_len[0];
 774         fcb->compute_feat = feat_copy;
 775     }
 776     else {
 777         int32 i, l, k;
 778         char *strp;
 779         char *mtype = ckd_salloc(type);
 780         char *wd = ckd_salloc(type);
 781         /*
 782          * Generic definition: Format should be %d,%d,%d,...,%d (i.e.,
 783          * comma separated list of feature stream widths; #items =
 784          * #streams).  An optional window size (frames will be
 785          * concatenated) is also allowed, which can be specified with
 786          * a colon after the list of feature streams.
 787          */
 788         l = strlen(mtype);
 789         k = 0;
 790         for (i = 1; i < l - 1; i++) {
 791             if (mtype[i] == ',') {
 792                 mtype[i] = ' ';
 793                 k++;
 794             }
 795             else if (mtype[i] == ':') {
 796                 mtype[i] = '\0';
 797                 fcb->window_size = atoi(mtype + i + 1);
 798                 break;
 799             }
 800         }
 801         k++;                    /* Presumably there are (#commas+1) streams */
 802         fcb->n_stream = k;
 803         fcb->stream_len = (int32 *) ckd_calloc(k, sizeof(int32));
 804
 805         /* Scan individual feature stream lengths */
 806         strp = mtype;
 807         i = 0;
 808         fcb->out_dim = 0;
 809         fcb->cepsize = 0;
 810         while (sscanf(strp, "%s%n", wd, &l) == 1) {
 811             strp += l;
 812             if ((i >= fcb->n_stream)
 813                 || (sscanf(wd, "%d", &(fcb->stream_len[i])) != 1)
 814                 || (fcb->stream_len[i] <= 0))
 815                 E_FATAL("Bad feature type argument\n");
 816             /* Input size before windowing */
 817             fcb->cepsize += fcb->stream_len[i];
 818             if (fcb->window_size > 0)
 819                 fcb->stream_len[i] *= (fcb->window_size * 2 + 1);
 820             /* Output size after windowing */
 821             fcb->out_dim += fcb->stream_len[i];
 822             i++;
 823         }
 824         if (i != fcb->n_stream)
 825             E_FATAL("Bad feature type argument\n");
 826         if (fcb->cepsize != cepsize)
 827             E_FATAL("Bad feature type argument\n");
 828
 829         /* Input is already the feature stream */
 830         fcb->compute_feat = feat_copy;
 831         ckd_free(mtype);
 832         ckd_free(wd);
 833     }
 834
 835     if (cmn != CMN_NONE)
 836         fcb->cmn_struct = cmn_init(feat_cepsize(fcb));
 837     fcb->cmn = cmn;
 838     fcb->varnorm = varnorm;
 839     if (agc != AGC_NONE) {
 840         fcb->agc_struct = agc_init();
 841         /*
 842          * No need to check if agc is set to EMAX; agc_emax_set() changes only emax related things
 843          * Moreover, if agc is not NONE and block mode is used, feat_agc() SILENTLY
 844          * switches to EMAX
 845          */
 846         /* HACK: hardwired initial estimates based on use of CMN (from Sphinx2) */
 847         agc_emax_set(fcb->agc_struct, (cmn != CMN_NONE) ? 5.0 : 10.0);
 848     }
 849     fcb->agc = agc;
 850     /*
 851      * Make sure this buffer is large enough to be used in feat_s2mfc2feat_block_utt()
 852      */
 853     fcb->cepbuf = (mfcc_t **) ckd_calloc_2d((LIVEBUFBLOCKSIZE < feat_window_size(fcb) * 2) ? feat_window_size(fcb) * 2 : LIVEBUFBLOCKSIZE,
 854                                             feat_cepsize(fcb),
 855                                             sizeof(mfcc_t));
 856     /* This one is actually just an array of pointers to "flatten out"
 857      * wraparounds. */
 858     fcb->tmpcepbuf = ckd_calloc(2 * feat_window_size(fcb) + 1,
 859                                 sizeof(*fcb->tmpcepbuf));
 860
 861     return fcb;
 862 }
 863
 864
 865 void
 866 feat_print(feat_t * fcb, mfcc_t *** feat, int32 nfr, FILE * fp)
 867 {
 868     int32 i, j, k;
 869
 870     for (i = 0; i < nfr; i++) {
 871         fprintf(fp, "%8d:\n", i);
 872
 873         for (j = 0; j < feat_dimension1(fcb); j++) {
 874             fprintf(fp, "\t%2d:", j);
 875
 876             for (k = 0; k < feat_dimension2(fcb, j); k++)
 877                 fprintf(fp, " %8.4f", MFCC2FLOAT(feat[i][j][k]));
 878             fprintf(fp, "\n");
 879         }
 880     }
 881
 882     fflush(fp);
 883 }
 884
 885 static void
 886 feat_cmn(feat_t *fcb, mfcc_t **mfc, int32 nfr, int32 beginutt, int32 endutt)
 887 {
 888     cmn_type_t cmn_type = fcb->cmn;
 889
 890     if (!(beginutt && endutt)
 891         && cmn_type != CMN_NONE) /* Only cmn_prior in block computation mode. */
 892         cmn_type = CMN_PRIOR;
 893
 894     switch (cmn_type) {
 895     case CMN_CURRENT:
 896         cmn(fcb->cmn_struct, mfc, fcb->varnorm, nfr);
 897         break;
 898     case CMN_PRIOR:
 899         cmn_prior(fcb->cmn_struct, mfc, fcb->varnorm, nfr);
 900         if (endutt)
 901             cmn_prior_update(fcb->cmn_struct);
 902         break;
 903     default:
 904         ;
 905     }
 906     cep_dump_dbg(fcb, mfc, nfr, "After CMN");
 907 }
 908
 909 static void
 910 feat_agc(feat_t *fcb, mfcc_t **mfc, int32 nfr, int32 beginutt, int32 endutt)
 911 {
 912     agc_type_t agc_type = fcb->agc;
 913
 914     if (!(beginutt && endutt)
 915         && agc_type != AGC_NONE) /* Only agc_emax in block computation mode. */
 916         agc_type = AGC_EMAX;
 917
 918     switch (agc_type) {
 919     case AGC_MAX:
 920         agc_max(fcb->agc_struct, mfc, nfr);
 921         break;
 922     case AGC_EMAX:
 923         agc_emax(fcb->agc_struct, mfc, nfr);
 924         if (endutt)
 925             agc_emax_update(fcb->agc_struct);
 926         break;
 927     case AGC_NOISE:
 928         agc_noise(fcb->agc_struct, mfc, nfr);
 929         break;
 930     default:
 931         ;
 932     }
 933     cep_dump_dbg(fcb, mfc, nfr, "After AGC");
 934 }
 935
 936 static void
 937 feat_compute_utt(feat_t *fcb, mfcc_t **mfc, int32 nfr, int32 win, mfcc_t ***feat)
 938 {
 939     int32 i;
 940
 941     cep_dump_dbg(fcb, mfc, nfr, "Incoming features (after padding)");
 942
 943     /* Create feature vectors */
 944     for (i = win; i < nfr - win; i++) {
 945         fcb->compute_feat(fcb, mfc + i, feat[i - win]);
 946     }
 947
 948     feat_print_dbg(fcb, feat, nfr - win * 2, "After dynamic feature computation");
 949
 950     if (fcb->lda) {
 951         feat_lda_transform(fcb, feat, nfr - win * 2);
 952         feat_print_dbg(fcb, feat, nfr - win * 2, "After LDA");
 953     }
 954
 955     if (fcb->subvecs) {
 956         feat_subvec_project(fcb, feat, nfr - win * 2);
 957         feat_print_dbg(fcb, feat, nfr - win * 2, "After subvector projection");
 958     }
 959 }
 960
 961
 962 /**
 963  * Read Sphinx-II format mfc file (s2mfc = Sphinx-II format MFC data).
 964  * If out_mfc is NULL, no actual reading will be done, and the number of
 965  * frames (plus padding) that would be read is returned.
 966  *
 967  * It's important that normalization is done before padding because
 968  * frames outside the data we are interested in shouldn't be taken
 969  * into normalization stats.
 970  *
 971  * @return # frames read (plus padding) if successful, -1 if
 972  * error (e.g., mfc array too small).
 973  */
 974 static int32
 975 feat_s2mfc_read_norm_pad(feat_t *fcb, char *file, int32 win,
 976                          int32 sf, int32 ef,
 977                          mfcc_t ***out_mfc,
 978                          int32 maxfr,
 979                          int32 cepsize)
 980 {
 981     FILE *fp;
 982     int32 n_float32;
 983     float32 *float_feat;
 984     struct stat statbuf;
 985     int32 i, n, byterev;
 986     int32 start_pad, end_pad;
 987     mfcc_t **mfc;
 988
 989     /* Initialize the output pointer to NULL, so that any attempts to
 990        free() it if we fail before allocating it will not segfault! */
 991     if (out_mfc)
 992         *out_mfc = NULL;
 993     E_INFO("Reading mfc file: '%s'[%d..%d]\n", file, sf, ef);
 994     if (ef >= 0 && ef <= sf) {
 995         E_ERROR("%s: End frame (%d) <= Start frame (%d)\n", file, ef, sf);
 996         return -1;
 997     }
 998
 999     /* Find filesize; HACK!! To get around intermittent NFS failures, use stat_retry */
1000     if ((stat_retry(file, &statbuf) < 0)
1001         || ((fp = fopen(file, "rb")) == NULL)) {
1002         E_ERROR("Failed to open file '%s' for reading: %s\n", file, strerror(errno));
1003         return -1;
1004     }
1005
1006     /* Read #floats in header */
1007     if (fread_retry(&n_float32, sizeof(int32), 1, fp) != 1) {
1008         E_ERROR("%s: fread(#floats) failed\n", file);
1009         fclose(fp);
1010         return -1;
1011     }
1012
1013     /* Check if n_float32 matches file size */
1014     byterev = 0;
1015     if ((int32) (n_float32 * sizeof(float32) + 4) != (int32) statbuf.st_size) { /* RAH, typecast both sides to remove compile warning */
1016         n = n_float32;
1017         SWAP_INT32(&n);
1018
1019         if ((int32) (n * sizeof(float32) + 4) != (int32) (statbuf.st_size)) {   /* RAH, typecast both sides to remove compile warning */
1020             E_ERROR
1021                 ("%s: Header size field: %d(%08x); filesize: %d(%08x)\n",
1022                  file, n_float32, n_float32, statbuf.st_size,
1023                  statbuf.st_size);
1024             fclose(fp);
1025             return -1;
1026         }
1027
1028         n_float32 = n;
1029         byterev = 1;
1030     }
1031     if (n_float32 <= 0) {
1032         E_ERROR("%s: Header size field (#floats) = %d\n", file, n_float32);
1033         fclose(fp);
1034         return -1;
1035     }
1036
1037     /* Convert n to #frames of input */
1038     n = n_float32 / cepsize;
1039     if (n * cepsize != n_float32) {
1040         E_ERROR("Header size field: %d; not multiple of %d\n", n_float32,
1041                 cepsize);
1042         fclose(fp);
1043         return -1;
1044     }
1045
1046     /* Check start and end frames */
1047     if (sf > 0) {
1048         if (sf >= n) {
1049             E_ERROR("%s: Start frame (%d) beyond file size (%d)\n", file,
1050                     sf, n);
1051             fclose(fp);
1052             return -1;
1053         }
1054     }
1055     if (ef < 0)
1056         ef = n-1;
1057     else if (ef >= n) {
1058         E_WARN("%s: End frame (%d) beyond file size (%d), will truncate\n",
1059                file, ef, n);
1060         ef = n-1;
1061     }
1062
1063     /* Add window to start and end frames */
1064     sf -= win;
1065     ef += win;
1066     if (sf < 0) {
1067         start_pad = -sf;
1068         sf = 0;
1069     }
1070     else
1071         start_pad = 0;
1072     if (ef >= n) {
1073         end_pad = ef - n + 1;
1074         ef = n - 1;
1075     }
1076     else
1077         end_pad = 0;
1078
1079     /* Limit n if indicated by [sf..ef] */
1080     if ((ef - sf + 1) < n)
1081         n = (ef - sf + 1);
1082     if (maxfr > 0 && n + start_pad + end_pad > maxfr) {
1083         E_ERROR("%s: Maximum output size(%d frames) < actual #frames(%d)\n",
1084                 file, maxfr, n + start_pad + end_pad);
1085         fclose(fp);
1086         return -1;
1087     }
1088
1089     /* If no output buffer was supplied, then skip the actual data reading. */
1090     if (out_mfc != NULL) {
1091         /* Position at desired start frame and read actual MFC data */
1092         mfc = (mfcc_t **)ckd_calloc_2d(n + start_pad + end_pad, cepsize, sizeof(mfcc_t));
1093         if (sf > 0)
1094             fseek(fp, sf * cepsize * sizeof(float32), SEEK_CUR);
1095         n_float32 = n * cepsize;
1096 #ifdef FIXED_POINT
1097         float_feat = ckd_calloc(n_float32, sizeof(float32));
1098 #else
1099         float_feat = mfc[start_pad];
1100 #endif
1101         if (fread_retry(float_feat, sizeof(float32), n_float32, fp) != n_float32) {
1102             E_ERROR("%s: fread(%dx%d) (MFC data) failed\n", file, n, cepsize);
1103             ckd_free_2d(mfc);
1104             fclose(fp);
1105             return -1;
1106         }
1107         if (byterev) {
1108             for (i = 0; i < n_float32; i++) {
1109                 SWAP_FLOAT32(&float_feat[i]);
1110             }
1111         }
1112 #ifdef FIXED_POINT
1113         for (i = 0; i < n_float32; ++i) {
1114             mfc[start_pad][i] = FLOAT2MFCC(float_feat[i]);
1115         }
1116         ckd_free(float_feat);
1117 #endif
1118
1119         /* Normalize */
1120         feat_cmn(fcb, mfc + start_pad, n, 1, 1);
1121         feat_agc(fcb, mfc + start_pad, n, 1, 1);
1122
1123         /* Replicate start and end frames if necessary. */
1124         for (i = 0; i < start_pad; ++i)
1125             memcpy(mfc[i], mfc[start_pad], cepsize * sizeof(mfcc_t));
1126         for (i = 0; i < end_pad; ++i)
1127             memcpy(mfc[start_pad + n + i], mfc[start_pad + n - 1],
1128                    cepsize * sizeof(mfcc_t));
1129
1130         *out_mfc = mfc;
1131     }
1132
1133     fclose(fp);
1134     return n + start_pad + end_pad;
1135 }
1136
1137
1138
1139 int32
1140 feat_s2mfc2feat(feat_t * fcb, const char *file, const char *dir, const char *cepext,
1141                 int32 sf, int32 ef, mfcc_t *** feat, int32 maxfr)
1142 {
1143     char *path;
1144     char *ps = "/";
1145     int32 win, nfr;
1146     int32 file_length, cepext_length, path_length = 0;
1147     mfcc_t **mfc;
1148
1149     if (fcb->cepsize <= 0) {
1150         E_ERROR("Bad cepsize: %d\n", fcb->cepsize);
1151         return -1;
1152     }
1153
1154     if (cepext == NULL)
1155         cepext = "";
1156
1157     /*
1158      * Create mfc filename, combining file, dir and extension if
1159      * necessary
1160      */
1161
1162     /*
1163      * First we decide about the path. If dir is defined, then use
1164      * it. Otherwise assume the filename already contains the path.
1165      */
1166     if (dir == NULL) {
1167         dir = "";
1168         ps = "";
1169         /*
1170          * This is not true but some 3rd party apps
1171          * may parse the output explicitly checking for this line
1172          */
1173         E_INFO("At directory . (current directory)\n");
1174     }
1175     else {
1176         E_INFO("At directory %s\n", dir);
1177         /*
1178          * Do not forget the path separator!
1179          */
1180         path_length += strlen(dir) + 1;
1181     }
1182
1183     /*
1184      * Include cepext, if it's not already part of the filename.
1185      */
1186     file_length = strlen(file);
1187     cepext_length = strlen(cepext);
1188     if ((file_length > cepext_length)
1189         && (strcmp(file + file_length - cepext_length, cepext) == 0)) {
1190         cepext = "";
1191         cepext_length = 0;
1192     }
1193
1194     /*
1195      * Do not forget the '\0'
1196      */
1197     path_length += file_length + cepext_length + 1;
1198     path = (char*) ckd_calloc(path_length, sizeof(char));
1199
1200 #ifdef HAVE_SNPRINTF
1201     /*
1202      * Paranoia is our best friend...
1203      */
1204     while ((file_length = snprintf(path, path_length, "%s%s%s%s", dir, ps, file, cepext)) > path_length) {
1205         path_length = file_length;
1206         path = (char*) ckd_realloc(path, path_length * sizeof(char));
1207     }
1208 #else
1209     sprintf(path, "%s%s%s%s", dir, ps, file, cepext);
1210 #endif
1211
1212     win = feat_window_size(fcb);
1213     /* Pad maxfr with win, so we read enough raw feature data to
1214      * calculate the requisite number of dynamic features. */
1215     if (maxfr >= 0)
1216         maxfr += win * 2;
1217
1218     if (feat != NULL) {
1219         /* Read mfc file including window or padding if necessary. */
1220         nfr = feat_s2mfc_read_norm_pad(fcb, path, win, sf, ef, &mfc, maxfr, fcb->cepsize);
1221         ckd_free(path);
1222         if (nfr < 0) {
1223             ckd_free_2d((void **) mfc);
1224             return -1;
1225         }
1226
1227         /* Actually compute the features */
1228         feat_compute_utt(fcb, mfc, nfr, win, feat);
1229
1230         ckd_free_2d((void **) mfc);
1231     }
1232     else {
1233         /* Just calculate the number of frames we would need. */
1234         nfr = feat_s2mfc_read_norm_pad(fcb, path, win, sf, ef, NULL, maxfr, fcb->cepsize);
1235         ckd_free(path);
1236         if (nfr < 0)
1237             return nfr;
1238     }
1239
1240
1241     return (nfr - win * 2);
1242 }
1243
1244 static int32
1245 feat_s2mfc2feat_block_utt(feat_t * fcb, mfcc_t ** uttcep,
1246                           int32 nfr, mfcc_t *** ofeat)
1247 {
1248     mfcc_t **cepbuf;
1249     int32 i, win, cepsize;
1250
1251     win = feat_window_size(fcb);
1252     cepsize = feat_cepsize(fcb);
1253
1254     /* Copy and pad out the utterance (this requires that the
1255      * feature computation functions always access the buffer via
1256      * the frame pointers, which they do)  */
1257     cepbuf = ckd_calloc(nfr + win * 2, sizeof(mfcc_t *));
1258     memcpy(cepbuf + win, uttcep, nfr * sizeof(mfcc_t *));
1259
1260     /* Do normalization before we interpolate on the boundary */
1261     feat_cmn(fcb, cepbuf + win, nfr, 1, 1);
1262     feat_agc(fcb, cepbuf + win, nfr, 1, 1);
1263
1264     /* Now interpolate */
1265     for (i = 0; i < win; ++i) {
1266         cepbuf[i] = fcb->cepbuf[i];
1267         memcpy(cepbuf[i], uttcep[0], cepsize * sizeof(mfcc_t));
1268         cepbuf[nfr + win + i] = fcb->cepbuf[win + i];
1269         memcpy(cepbuf[nfr + win + i], uttcep[nfr - 1], cepsize * sizeof(mfcc_t));
1270     }
1271     /* Compute as usual. */
1272     feat_compute_utt(fcb, cepbuf, nfr + win * 2, win, ofeat);
1273     ckd_free(cepbuf);
1274     return nfr;
1275 }
1276
1277 int32
1278 feat_s2mfc2feat_live(feat_t * fcb, mfcc_t ** uttcep, int32 *inout_ncep,
1279                      int32 beginutt, int32 endutt, mfcc_t *** ofeat)
1280 {
1281     int32 win, cepsize, nbufcep;
1282     int32 i, j, nfeatvec;
1283     int32 zero = 0;
1284
1285     /* Avoid having to check this everywhere. */
1286     if (inout_ncep == NULL) inout_ncep = &zero;
1287
1288     /* Special case for entire utterances. */
1289     if (beginutt && endutt && *inout_ncep > 0)
1290         return feat_s2mfc2feat_block_utt(fcb, uttcep, *inout_ncep, ofeat);
1291
1292     win = feat_window_size(fcb);
1293     cepsize = feat_cepsize(fcb);
1294
1295     /* Empty the input buffer on start of utterance. */
1296     if (beginutt)
1297         fcb->bufpos = fcb->curpos;
1298
1299     /* Calculate how much data is in the buffer already. */
1300     nbufcep = fcb->bufpos - fcb->curpos;
1301     if (nbufcep < 0)
1302         nbufcep = fcb->bufpos + LIVEBUFBLOCKSIZE - fcb->curpos;
1303     /* Add any data that we have to replicate. */
1304     if (beginutt && *inout_ncep > 0)
1305         nbufcep += win;
1306     if (endutt)
1307         nbufcep += win;
1308
1309     /* Only consume as much input as will fit in the buffer. */
1310     if (nbufcep + *inout_ncep > LIVEBUFBLOCKSIZE) {
1311         /* We also can't overwrite the trailing window, hence the
1312          * reason why win is subtracted here. */
1313         *inout_ncep = LIVEBUFBLOCKSIZE - nbufcep - win;
1314         /* Cancel end of utterance processing. */
1315         endutt = FALSE;
1316     }
1317
1318     /* FIXME: Don't modify the input! */
1319     feat_cmn(fcb, uttcep, *inout_ncep, beginutt, endutt);
1320     feat_agc(fcb, uttcep, *inout_ncep, beginutt, endutt);
1321
1322     /* Replicate first frame into the first win frames if we're at the
1323      * beginning of the utterance and there was some actual input to
1324      * deal with.  (FIXME: Not entirely sure why that condition) */
1325     if (beginutt && *inout_ncep > 0) {
1326         for (i = 0; i < win; i++) {
1327             memcpy(fcb->cepbuf[fcb->bufpos++], uttcep[0],
1328                    cepsize * sizeof(mfcc_t));
1329             fcb->bufpos %= LIVEBUFBLOCKSIZE;
1330         }
1331         /* Move the current pointer past this data. */
1332         fcb->curpos = fcb->bufpos;
1333         nbufcep -= win;
1334     }
1335
1336     /* Copy in frame data to the circular buffer. */
1337     for (i = 0; i < *inout_ncep; ++i) {
1338         memcpy(fcb->cepbuf[fcb->bufpos++], uttcep[i],
1339                cepsize * sizeof(mfcc_t));
1340         fcb->bufpos %= LIVEBUFBLOCKSIZE;
1341         ++nbufcep;
1342     }
1343
1344     /* Replicate last frame into the last win frames if we're at the
1345      * end of the utterance (even if there was no input, so we can
1346      * flush the output). */
1347     if (endutt) {
1348         int32 tpos; /* Index of last input frame. */
1349         if (fcb->bufpos == 0)
1350             tpos = LIVEBUFBLOCKSIZE - 1;
1351         else
1352             tpos = fcb->bufpos - 1;
1353         for (i = 0; i < win; ++i) {
1354             memcpy(fcb->cepbuf[fcb->bufpos++], fcb->cepbuf[tpos],
1355                    cepsize * sizeof(mfcc_t));
1356             fcb->bufpos %= LIVEBUFBLOCKSIZE;
1357         }
1358     }
1359
1360     /* We have to leave the trailing window of frames. */
1361     nfeatvec = nbufcep - win;
1362     if (nfeatvec <= 0)
1363         return 0; /* Do nothing. */
1364
1365     for (i = 0; i < nfeatvec; ++i) {
1366         /* Handle wraparound cases. */
1367         if (fcb->curpos - win < 0 || fcb->curpos + win >= LIVEBUFBLOCKSIZE) {
1368             /* Use tmpcepbuf for this case.  Actually, we just need the pointers. */
1369             for (j = -win; j <= win; ++j) {
1370                 int32 tmppos =
1371                     (fcb->curpos + j + LIVEBUFBLOCKSIZE) % LIVEBUFBLOCKSIZE;
1372                 fcb->tmpcepbuf[win + j] = fcb->cepbuf[tmppos];
1373             }
1374             fcb->compute_feat(fcb, fcb->tmpcepbuf + win, ofeat[i]);
1375         }
1376         else {
1377             fcb->compute_feat(fcb, fcb->cepbuf + fcb->curpos, ofeat[i]);
1378         }
1379         /* Move the read pointer forward. */
1380         ++fcb->curpos;
1381         fcb->curpos %= LIVEBUFBLOCKSIZE;
1382     }
1383
1384     if (fcb->lda)
1385         feat_lda_transform(fcb, ofeat, nfeatvec);
1386
1387     if (fcb->subvecs)
1388         feat_subvec_project(fcb, ofeat, nfeatvec);
1389
1390     return nfeatvec;
1391 }
1392
1393 feat_t *
1394 feat_retain(feat_t *f)
1395 {
1396     ++f->refcount;
1397     return f;
1398 }
1399
1400 int
1401 feat_free(feat_t * f)
1402 {
1403     if (f == NULL)
1404         return 0;
1405     if (--f->refcount > 0)
1406         return f->refcount;
1407
1408     if (f->cepbuf)
1409         ckd_free_2d((void **) f->cepbuf);
1410     ckd_free(f->tmpcepbuf);
1411
1412     if (f->name) {
1413         ckd_free((void *) f->name);
1414     }
1415     if (f->lda)
1416         ckd_free_3d((void ***) f->lda);
1417
1418     ckd_free(f->stream_len);
1419     ckd_free(f->sv_len);
1420     ckd_free(f->sv_buf);
1421     subvecs_free(f->subvecs);
1422
1423     cmn_free(f->cmn_struct);
1424     agc_free(f->agc_struct);
1425
1426     ckd_free(f);
1427     return 0;
1428 }
1429
1430
1431 void
1432 feat_report(feat_t * f)
1433 {
1434     int i;
1435     E_INFO_NOFN("Initialization of feat_t, report:\n");
1436     E_INFO_NOFN("Feature type         = %s\n", f->name);
1437     E_INFO_NOFN("Cepstral size        = %d\n", f->cepsize);
1438     E_INFO_NOFN("Number of streams    = %d\n", f->n_stream);
1439     for (i = 0; i < f->n_stream; i++) {
1440         E_INFO_NOFN("Vector size of stream[%d]: %d\n", i,
1441                     f->stream_len[i]);
1442     }
1443     E_INFO_NOFN("Number of subvectors = %d\n", f->n_sv);
1444     for (i = 0; i < f->n_sv; i++) {
1445         int32 *sv;
1446
1447         E_INFO_NOFN("Components of subvector[%d]:", i);
1448         for (sv = f->subvecs[i]; sv && *sv != -1; ++sv)
1449             E_INFOCONT(" %d", *sv);
1450         E_INFOCONT("\n");
1451     }
1452     E_INFO_NOFN("Whether CMN is used  = %d\n", f->cmn);
1453     E_INFO_NOFN("Whether AGC is used  = %d\n", f->agc);
1454     E_INFO_NOFN("Whether variance is normalized = %d\n", f->varnorm);
1455     E_INFO_NOFN("\n");
1456 }