src/libsphinxbase/fe/fe_warp_piecewise_linear.c

   1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
   2 /* ====================================================================
   3  * Copyright (c) 2006 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  *
  39  * File: fe_warp_piecewise_linear.c
  40  *
  41  * Description:
  42  *
  43  *      Warp the frequency axis according to an piecewise linear
  44  *      function. The function is linear up to a frequency F, where
  45  *      the slope changes so that the Nyquist frequency in the warped
  46  *      axis maps to the Nyquist frequency in the unwarped.
  47  *
  48  *              w' = a * w, w < F
  49  *              w' = a' * w + b, W > F
  50  *              w'(0) = 0
  51  *              w'(F) = F
  52  *              w'(Nyq) = Nyq
  53  *
  54  *********************************************************************/
  55
  56 /* static char rcsid[] = "@(#)$Id: fe_warp_piecewise_linear.c,v 1.2 2006/02/17 00:31:34 egouvea Exp $"; */
  57
  58 #include <stdio.h>
  59 #include <stdlib.h>
  60 #include <math.h>
  61 #include <string.h>
  62
  63 #ifdef _MSC_VER
  64 #pragma warning (disable: 4996)
  65 #endif
  66
  67 #include "sphinxbase/strfuncs.h"
  68 #include "sphinxbase/err.h"
  69
  70 #include "fe_warp.h"
  71 #include "fe_warp_piecewise_linear.h"
  72
  73 #define N_PARAM         2
  74 #define YES             1
  75 #define NO              0
  76
  77 /*
  78  * params[0] : a
  79  * params[1] : F (the non-differentiable point)
  80  */
  81 static float params[N_PARAM] = { 1.0f, 6800.0f };
  82 static float final_piece[2];
  83 static int32 is_neutral = YES;
  84 static char p_str[256] = "";
  85 static float nyquist_frequency = 0.0f;
  86 \f
  87
  88 const char *
  89 fe_warp_piecewise_linear_doc()
  90 {
  91     return "piecewise_linear :== < w' = a * w, w < F >";
  92 }
  93 \f
  94 uint32
  95 fe_warp_piecewise_linear_id()
  96 {
  97     return FE_WARP_ID_PIECEWISE_LINEAR;
  98 }
  99 \f
 100 uint32
 101 fe_warp_piecewise_linear_n_param()
 102 {
 103     return N_PARAM;
 104 }
 105 \f
 106 void
 107 fe_warp_piecewise_linear_set_parameters(char const *param_str,
 108                                         float sampling_rate)
 109 {
 110     char *tok;
 111     char *seps = " \t";
 112     char temp_param_str[256];
 113     int param_index = 0;
 114
 115     nyquist_frequency = sampling_rate / 2;
 116     if (param_str == NULL) {
 117         is_neutral = YES;
 118         return;
 119     }
 120     /* The new parameters are the same as the current ones, so do nothing. */
 121     if (strcmp(param_str, p_str) == 0) {
 122         return;
 123     }
 124     is_neutral = NO;
 125     strcpy(temp_param_str, param_str);
 126     memset(params, 0, N_PARAM * sizeof(float));
 127     memset(final_piece, 0, 2 * sizeof(float));
 128     strcpy(p_str, param_str);
 129     /* FIXME: strtok() is not re-entrant... */
 130     tok = strtok(temp_param_str, seps);
 131     while (tok != NULL) {
 132         params[param_index++] = (float) atof_c(tok);
 133         tok = strtok(NULL, seps);
 134         if (param_index >= N_PARAM) {
 135             break;
 136         }
 137     }
 138     if (tok != NULL) {
 139         E_INFO
 140             ("Piecewise linear warping takes up to two arguments, %s ignored.\n",
 141              tok);
 142     }
 143     if (params[1] < sampling_rate) {
 144         /* Precompute these. These are the coefficients of a
 145          * straight line that contains the points (F, aF) and (N,
 146          * N), where a = params[0], F = params[1], N = Nyquist
 147          * frequency.
 148          */
 149         if (params[1] == 0) {
 150             params[1] = sampling_rate * 0.85f;
 151         }
 152         final_piece[0] =
 153             (nyquist_frequency -
 154              params[0] * params[1]) / (nyquist_frequency - params[1]);
 155         final_piece[1] =
 156             nyquist_frequency * params[1] * (params[0] -
 157                                          1.0f) / (nyquist_frequency -
 158                                                   params[1]);
 159     }
 160     else {
 161         memset(final_piece, 0, 2 * sizeof(float));
 162     }
 163     if (params[0] == 0) {
 164         is_neutral = YES;
 165         E_INFO
 166             ("Piecewise linear warping cannot have slope zero, warping not applied.\n");
 167     }
 168 }
 169 \f
 170 float
 171 fe_warp_piecewise_linear_warped_to_unwarped(float nonlinear)
 172 {
 173     if (is_neutral) {
 174         return nonlinear;
 175     }
 176     else {
 177         /* linear = (nonlinear - b) / a */
 178         float temp;
 179         if (nonlinear < params[0] * params[1]) {
 180             temp = nonlinear / params[0];
 181         }
 182         else {
 183             temp = nonlinear - final_piece[1];
 184             temp /= final_piece[0];
 185         }
 186         if (temp > nyquist_frequency) {
 187             E_WARN
 188                 ("Warp factor %g results in frequency (%.1f) higher than Nyquist (%.1f)\n",
 189                  params[0], temp, nyquist_frequency);
 190         }
 191         return temp;
 192     }
 193 }
 194 \f
 195 float
 196 fe_warp_piecewise_linear_unwarped_to_warped(float linear)
 197 {
 198     if (is_neutral) {
 199         return linear;
 200     }
 201     else {
 202         float temp;
 203         /* nonlinear = a * linear - b */
 204         if (linear < params[1]) {
 205             temp = linear * params[0];
 206         }
 207         else {
 208             temp = final_piece[0] * linear + final_piece[1];
 209         }
 210         return temp;
 211     }
 212 }
 213 \f
 214 void
 215 fe_warp_piecewise_linear_print(const char *label)
 216 {
 217     uint32 i;
 218
 219     for (i = 0; i < N_PARAM; i++) {
 220         printf("%s[%04u]: %6.3f ", label, i, params[i]);
 221     }
 222     printf("\n");
 223 }
 224 \f
 225 /*
 226  * Log record.  Maintained by RCS.
 227  *
 228  * $Log: fe_warp_piecewise_linear.c,v $
 229  * Revision 1.2  2006/02/17 00:31:34  egouvea
 230  * Removed switch -melwarp. Changed the default for window length to
 231  * 0.025625 from 0.256 (so that a window at 16kHz sampling rate has
 232  * exactly 410 samples). Cleaned up include's. Replaced some E_FATAL()
 233  * with E_WARN() and return.
 234  *
 235  * Revision 1.1  2006/02/16 00:18:26  egouvea
 236  * Implemented flexible warping function. The user can specify at run
 237  * time which of several shapes they want to use. Currently implemented
 238  * are an affine function (y = ax + b), an inverse linear (y = a/x) and a
 239  * piecewise linear (y = ax, up to a frequency F, and then it "breaks" so
 240  * Nyquist frequency matches in both scales.
 241  *
 242  * Added two switches, -warp_type and -warp_params. The first specifies
 243  * the type, which valid values:
 244  *
 245  * -inverse or inverse_linear
 246  * -linear or affine
 247  * -piecewise or piecewise_linear
 248  *
 249  * The inverse_linear is the same as implemented by EHT. The -mel_warp
 250  * switch was kept for compatibility (maybe remove it in the
 251  * future?). The code is compatible with EHT's changes: cepstra created
 252  * from code after his changes should be the same as now. Scripts that
 253  * worked with his changes should work now without changes. Tested a few
 254  * cases, same results.
 255  *
 256  */