1/*---------------------------------------------------------------------------*
2 *  chelmel4.c  *
3 *                                                                           *
4 *  Copyright 2007, 2008 Nuance Communciations, Inc.                               *
5 *                                                                           *
6 *  Licensed under the Apache License, Version 2.0 (the 'License');          *
7 *  you may not use this file except in compliance with the License.         *
8 *                                                                           *
9 *  You may obtain a copy of the License at                                  *
10 *      http://www.apache.org/licenses/LICENSE-2.0                           *
11 *                                                                           *
12 *  Unless required by applicable law or agreed to in writing, software      *
13 *  distributed under the License is distributed on an 'AS IS' BASIS,        *
14 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
15 *  See the License for the specific language governing permissions and      *
16 *  limitations under the License.                                           *
17 *                                                                           *
18 *---------------------------------------------------------------------------*/
19
20
21
22#include <stdlib.h>
23#ifndef _RTT
24#include <stdio.h>
25#endif
26#include <string.h>
27#include <math.h>
28#include <limits.h>
29#include <assert.h>
30
31#include "hmm_desc.h"
32#include "front.h"
33#include "portable.h"
34
35#define DEBUG           0
36#define LOG_AS_PRINT 0
37
38#define LPCMAX 100
39
40
41#if LOG_AS_PRINT /* BP temp hack */
42#define log_report printf
43#endif
44
45#include "sh_down.h"
46
47
48/* cepstrum_params has been broken into three functions:
49 filterbank_emulation    - does preemp, window, fft and filtbank
50 gain_adjustment                 - estimates gain
51 cepstrum_params                 - does gain adj.(if on), spec corr and cos transform
52   This enables us to bypass gain adjustment.
53*/
54
55//static void mel_cuberoot_offset(cepdata *fbo, cepdata *ch_off, int nf);
56#if SPEC_CORRECT
57static void mel_spectrum_correction(cepdata *fbo, cepdata *ch_gain, int nf);
58#endif
59static void mel_loglookup_with_offset(front_cep *cepobj,
60                                      front_channel *channel);
61//static void mel_exp(cepdata *fbo, int nf);
62//static void durbin(cepdata *a, cepdata *r, int n);
63//static void lpc_to_cepstral_recursion(cepdata *c, cepdata *a, int nc, int n);
64static void icostrans(cepdata *cs, cepdata fb[], cepdata cep[], int nf, int nc);
65
66
67void cepstrum_params(front_channel *channel, front_wave *waveobj,
68                     front_freq *freqobj, front_cep *cepobj)
69{
70#if SPEC_CORRECT
71  /*  2.30 Apply a spectrum correction */
72  if (cepobj->mel_loop)
73    mel_spectrum_correction(freqobj->filterbank, cepobj->mel_loop, channel->num_freq);
74#endif
75  /*  2.33 Calculate log dB energy values */
76    mel_loglookup_with_offset(cepobj, channel);
77#if DEBUG
78  log_report("Filterbank output: ");
79  write_scaled_frames(freqobj->nf, 1, channel->filterbank, D_FIXED, 1 / (float)LOG_SCALE);
80#endif
81
82  /*  2.34 Cosine transformation */
83  icostrans(cepobj->cs, channel->filterbank, channel->cep,
84            channel->num_freq, cepobj->mel_dim);
85
86#if DEBUG
87  log_report("Cepstrum coefficients: ");
88  write_scaled_frames((cepobj->mel_dim + 1), 1, channel->cep, D_FIXED, (float)1 / (0x01 << (LOG_SCALE_SHIFT + COSINE_TABLE_SHIFT)));
89#endif
90  return;
91}
92
93static void icostrans(cepdata *cs, cepdata fb[], cepdata cep[], int nf, int nc)
94/*
95**  inv rotated-cosine transform
96**  ref Davis and Mermelstein, ASSP 1980 */
97{
98  int   i, j;
99  cepdata *cp;
100  cepdata temp;
101
102  for (i = 0; i <= nc; i++)
103  {
104    cp = &cs[i*nf];
105    for (j = 0, temp = 0; j < nf; j++)
106      temp += fb[j] * cp[j];
107    cep[i] = temp;
108  }
109  return;
110}
111
112#if SPEC_CORRECT
113static void mel_spectrum_correction(cepdata *fbo, cepdata *ch_gain, int nf)
114/*
115**  pwr spect -> filter bank output */
116{
117  int i;
118
119  for (i = 0;i < nf; i++)
120    fbo[i] = fbo[i] * ch_gain[i]; /* TODO: Fixedpt scale up and down */
121  return;
122}
123#endif
124
125static void mel_loglookup_with_offset(front_cep *cepobj,
126                                      front_channel *channel)
127/*
128**  pwr spect -> filter bank output */
129{
130  int ii;
131
132  if (channel->shift > 0)
133    for (ii = 0; ii < channel->num_freq; ii++)
134    {
135      channel->filterbank[ii] = (cepdata) log_lookup(&cepobj->logtab,
136                                (int)(channel->filterbank[ii] +
137                                      SHIFT_DOWN(cepobj->mel_offset[ii], channel->shift)),
138                                channel->shift);
139    }
140  else
141    for (ii = 0; ii < channel->num_freq; ii++)
142    {
143      channel->filterbank[ii] = (cepdata) log_lookup(&cepobj->logtab,
144                                (int)(channel->filterbank[ii] +
145                                      SHIFT_UP(cepobj->mel_offset[ii], -channel->shift)),
146                                channel->shift);
147    }
148
149  return;
150}
151
152//static void mel_exp(cepdata *fbo, int nf)
153///*
154//**  pwr spect -> filter bank output */
155//{
156//  int i;
157//  for (i = 0; i < nf; i++)
158//  {
159//    fbo[i] = (cepdata) exp((double) fbo[i]);
160//  }
161//  return;
162//}
163//
164//static void durbin(
165//  cepdata *a,   /* lpc coefficients           */
166//  cepdata *r,   /* autocorrelation coefficients       */
167//  int n)     /* order of lpc analysis        */
168//{
169//  int i, j;
170//  cepdata A[LPCMAX+1][LPCMAX+1], sum;
171//  cepdata k[LPCMAX+1];
172//  cepdata e[LPCMAX+1];
173//
174//  e[0] = r[0];
175//  for (i = 1; i <= n; i++)
176//  {
177//    sum = 0;
178//    for (j = 1; j <= (i - 1); j++)
179//    {
180//      sum += A[j][i-1] * r[i-j];
181//    }
182//    k[i] = -(r[i] + sum) / e[i-1];
183//    A[i][i] = k[i] ;
184//    for (j = 1; j <= (i - 1); j++)
185//    {
186//      A[j][i] = A[j][i-1] + k[i] * A[i-j][i-1];
187//    }
188//    e[i] = (1 - k[i] * k[i]) * e[i-1];
189//  }
190//  for (j = 1 ; j <= n; j++)
191//  {
192//    a[j] = A[j][n];
193//  }
194//
195//  a[0] = (cepdata) 1.0 ;
196//  return;
197//}
198//
199//static void lpc_to_cepstral_recursion(
200//  cepdata *c,     /* cepstral coefficients        */
201//  cepdata *a,     /* lpc coeffiecients            */
202//  int nc,         /* order of cepstra             */
203//  int n)          /* order of lpc                 */
204//{
205//  int k, i;
206//  cepdata sum;
207//
208//  ASSERT(nc < LPCMAX);
209//
210//  for (i = n + 1; i <= nc; i++)
211//  {
212//    a[i] = (cepdata) 0.0;
213//  }
214//  /* if lpc order less    */
215//  /* than cepstral order  */
216//  /* define higher lpccos */
217//
218//  for (i = 1; i <= nc; i++)
219//  {
220//    sum = (cepdata) 0.0;
221//    for (k = 1; k <= (i - 1); k++)
222//    {
223//      sum = sum + k * c[k] * a[i-k]; /* TODO: fixedpt range for mult */
224//    }
225//    c[i] = -a[i] - (sum / i);               /* cepstral calculated  */
226//    /* to <=nc in icostrans */
227//    /* so I shall do the    */                                                      /* same here            */
228//  }
229//}
230
231