Loading...
Searching...
No Matches
sigpr_frame.cc
/*************************************************************************/
/* */
/* Centre for Speech Technology Research */
/* University of Edinburgh, UK */
/* Copyright (c) 1996 */
/* All Rights Reserved. */
/* */
/* Permission is hereby granted, free of charge, to use and distribute */
/* this software and its documentation without restriction, including */
/* without limitation the rights to use, copy, modify, merge, publish, */
/* distribute, sublicense, and/or sell copies of this work, and to */
/* permit persons to whom this work is furnished to do so, subject to */
/* the following conditions: */
/* 1. The code must retain the above copyright notice, this list of */
/* conditions and the following disclaimer. */
/* 2. Any modifications must be clearly marked as such. */
/* 3. Original authors' names are not deleted. */
/* 4. The authors' names are not used to endorse or promote products */
/* derived from this software without specific prior written */
/* permission. */
/* */
/* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
/* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
/* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
/* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
/* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
/* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
/* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
/* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
/* THIS SOFTWARE. */
/* */
/*************************************************************************/
/* Authors: Paul Taylor and Simon King */
/* Date : March 1998 */
/*-----------------------------------------------------------------------*/
/* Functions operating on a single frame */
/* */
/*=======================================================================*/
#include "sigpr/EST_sigpr_frame.h"
#include "sigpr/EST_fft.h"
#include "EST_inline_utils.h"
#include "EST_math.h"
#include "EST_error.h"
#include "EST_TBuffer.h"
#define ALMOST1 0.99999
#define MAX_ABS_CEPS 4.0
float Hz2Mel(float frequency_in_Hertz);
float Mel2Hz(float frequency_in_Mel);
float lpredict(float *adc, int wsize,
float *acf, float *ref, float *lpc,
int order);
float ogi_lpc(float *adc, int wsize, int order,
float *acf, float *ref, float *lpc);
/*
void lpc2ref(float const *a, float *b, int c)
{
(void) a;
(void) b;
(void) c;
}
*/
void convert2lpc(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame)
{
if (in_type == "sig")
sig2lpc(in_frame, out_frame);
else if (in_type == "lsf")
lsf2lpc(in_frame, out_frame);
else if (in_type == "ref")
ref2lpc(in_frame, out_frame);
else
EST_error("Cannot convert coefficient type %s to lpc\n",
(const char *) in_type);
}
void convert2ref(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame)
{
EST_FVector tmp;
if (in_type == "lpc")
lpc2ref(in_frame, out_frame);
else if (in_type == "sig")
{
tmp.resize(out_frame.length());
sig2lpc(in_frame, tmp);
lpc2ref(tmp, out_frame);
}
else if (in_type == "lsf")
{
tmp.resize(out_frame.length());
lsf2lpc(in_frame, tmp);
lpc2ref(tmp, out_frame);
}
else
EST_error("Cannot convert coefficient type %s to reflection coefs\n",
(const char *) in_type);
}
void convert2area(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame)
{
EST_FVector tmp;
if (in_type == "lpc")
lpc2ref(in_frame, out_frame);
else if (in_type == "sig")
{
tmp.resize(out_frame.length());
sig2lpc(in_frame, tmp);
lpc2ref(tmp, out_frame);
}
else if (in_type == "lsf")
{
tmp.resize(out_frame.length());
lsf2lpc(in_frame, tmp);
lpc2ref(tmp, out_frame);
}
else
EST_error("Cannot convert coefficient type %s to reflection coefs\n",
(const char *) in_type);
}
void convert2cep(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame)
{
EST_FVector tmp;
if (in_type == "lpc")
lpc2cep(in_frame, out_frame);
else if (in_type == "sig")
{
tmp.resize(out_frame.length());
sig2lpc(in_frame, tmp);
lpc2cep(tmp, out_frame);
}
else if (in_type == "lsf")
{
tmp.resize(out_frame.length());
lsf2lpc(in_frame, tmp);
lpc2cep(tmp, out_frame);
}
else if (in_type == "ref")
{
tmp.resize(out_frame.length());
ref2lpc(in_frame, tmp);
lpc2cep(tmp, out_frame);
}
else
EST_error("Cannot convert coefficient type %s to cepstrum coefs\n",
(const char *) in_type);
}
// void convert2melcep(const EST_FVector &in_frame, const EST_String &in_type,
// EST_FVector &out_frame, int num_mfccs, int fbank_order,
// float liftering_parameter)
// {
// EST_FVector tmp;
// if (in_type == "fbank")
// // fbank2melcep(in_frame, out_frame);
// return;
// else if (in_type == "sig")
// {
// tmp.resize(out_frame.length());
// // incomplete !
// // sig2melcep(in_frame, outframe,num_mfccs,fbank_order,liftering_parameter);
// }
// else
// EST_error("Cannot convert coefficient type %s to Mel cepstrum coefs\n",
// (const char *) in_type);
// }
void convert2lsf(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame)
{
EST_FVector tmp;
if (in_type == "lpc")
lpc2lsf(in_frame, out_frame);
else if (in_type == "sig")
{
tmp.resize(out_frame.length());
sig2lpc(in_frame, tmp);
lpc2lsf(tmp, out_frame);
}
else if (in_type == "ref")
{
tmp.resize(out_frame.length());
ref2lpc(in_frame, tmp);
lpc2lsf(tmp, out_frame);
}
else
EST_error("Cannot convert coefficient type %s to reflection coefs\n",
(const char *) in_type);
}
void frame_convert(const EST_FVector &in_frame, const EST_String &in_type,
EST_FVector &out_frame, const EST_String &out_type)
{
if (out_type == "lpc")
convert2lpc(in_frame, in_type, out_frame);
else if (out_type == "lsf")
convert2lsf(in_frame, in_type, out_frame);
else if (out_type == "ref")
convert2ref(in_frame, in_type, out_frame);
else if (out_type == "cep")
convert2cep(in_frame, in_type, out_frame);
else if (out_type == "area")
convert2area(in_frame, in_type, out_frame);
else
EST_error("Cannot convert coefficients to type %s\n",
(const char *) out_type);
}
void sig2lpc(const EST_FVector &sig, EST_FVector &acf,
EST_FVector &ref, EST_FVector &xlpc);
float lpredict2(EST_FVector &adc, int wsize,
EST_FVector &acf, float *ref, float *lpc,
int order);
void sig2lpc(const EST_FVector &sig, EST_FVector &lpc)
{
EST_FVector acf(lpc.length()), ref(lpc.length());
/* float *fadc, *facf, *fref, *flpc;
fadc = new float[sig.length()];
facf = new float[lpc.length()];
fref = new float[lpc.length()];
flpc = new float[lpc.length()];
// for (int i = 0; i < sig.length(); ++i)
// adc[i] = sig(i);
lpredict2(sig, sig.length(), acf, fref, flpc, lpc.length()-1);
for (int i = 0; i < lpc.length(); ++i)
lpc.a(i) = flpc[i];
*/
sig2lpc(sig, acf, ref, lpc);
}
void sig2ref(const EST_FVector &sig, EST_FVector &ref)
{
(void)sig;
EST_FVector acf(ref.length()), lpc(ref.length());
// sig2lpc(sig, acf, ref, lpc);
}
void ref2area(const EST_FVector &ref, EST_FVector &area)
{
for (int i = 1; i < ref.length(); i++)
area[i] = (1.0 - ref(i)) / (1.0 + ref(i));
}
void ref2logarea(const EST_FVector &ref, EST_FVector &logarea)
{
int order = ref.length() -1;
for(int i = 1; i <= order; i++)
{
if (ref(i) > ALMOST1)
logarea[i] = log((1.0 - ALMOST1) / (1.0 + ALMOST1));
else if (ref(i) < -ALMOST1)
logarea[i] = log((1.0 + ALMOST1)/(1.0 - ALMOST1));
else
logarea[i] = log((1.0 - ref(i)) / (1.0 + ref(i)));
}
}
void ref2truearea(const EST_FVector &ref, EST_FVector &area)
{
int order = ref.length() -1;
area[1] = (1.0 - ref(1)) / (1.0 + ref(1));
for(int i = 2; i <= order; i++)
area[i] = area[i - 1] * (1.0 - ref(i)) / (1.0 + ref(i));
}
void lpc2cep(const EST_FVector &lpc, EST_FVector &cep)
{
int n, k;
float sum;
int order = lpc.length() - 1;
for (n = 1; n <= order && n <= cep.length(); n++)
{
sum = 0.0;
for (k = 1; k < n; k++)
sum += k * cep(k-1) * lpc(n - k);
cep[n-1] = lpc(n) + sum / n;
}
/* be wary of these interpolated values */
for(n = order + 1; n <= cep.length(); n++)
{
sum = 0.0;
for (k = n - (order - 1); k < n; k++)
sum += k * cep(k-1) * lpc(n - k);
cep[n-1] = sum / n;
}
/* very occasionally the above can go unstable, fudge if this happens */
for (n = 0; n < cep.length(); n++)
{
// check if NaN -- happens on some frames of silence
if (isnanf(cep[n]) ) cep[n] = 0.0;
if (cep[n] > MAX_ABS_CEPS){
cerr << "WARNING : cepstral coeff " << n << " was " <<
cep[n] << endl;
cerr << "lpc coeff " << n << " = " << lpc(n + 1) << endl;
cep[n] = MAX_ABS_CEPS;
}
if (cep[n] < -MAX_ABS_CEPS) {
cerr << "WARNING : cepstral coeff " << n << " was " <<
cep[n] << endl;
cep[n] = -MAX_ABS_CEPS;
}
}
}
// REORG - test this!!
void lpc2ref(const EST_FVector &lpc, EST_FVector &ref)
{
// seem to get weird output from this - best not to use it !
EST_error("lpc2ref Code unfinished\n");
// LPC to reflection coefficients
// from code from Borja Etxebarria
// This code does clever things with pointer and so has been
// left using float * arrays.
// simonk (May 99) : fixed because lpc coeffs always have energy at
// coeff 0 - the code here would need changing is lpc coeff 0 was
// ever made optional.
int lpc_offset=1;
int order = lpc.length() - 1;
int i,j;
float f,ai;
float *vo,*vx;
float *vn = new float[order];
i = order - 1;
ref[i] = lpc(i+lpc_offset);
ai = lpc(i+lpc_offset);
f = 1-ai*ai;
i--;
for (j=0; j<=i; j++)
ref[j] = (lpc(j+lpc_offset)+((ai*lpc(i-j+lpc_offset))))/f;
/* vn=vtmp in previous #define */
// Check whether this should really be a pointer
vo = new float[order];
for (i = 0; i < order; ++i)
vo[i] = ref(i);
for ( ;i>0; )
{
ai=vo[i];
f = 1-ai*ai;
i--;
for (j=0; j<=i; j++)
vn[j] = (vo[j]+((ai*vo[i-j])))/f;
ref[i]=vn[i];
vx = vn;
vn = vo;
vo = vx;
}
delete [] vn;
}
void ref2lpc(const EST_FVector &ref, EST_FVector &lpc)
{
// Here we use Christopher Longet Higgin's algorithm converted to
// an equivalent by awb. It doesn't have the reverse order or
// negation requirement.
int order = ref.length() - 1;
float a, b;
int n, k;
for (n=0; n < order; n++)
{
lpc[n] = ref(n);
for (k=0; 2 * (k+1) <= n + 1; k++)
{
a = lpc[k];
b = lpc[n-(k + 1)];
lpc[k] = a-b * lpc[n];
lpc[n-(k+1)] = b-a * lpc[n];
}
}
}
/************************************************************
** LPC_TO_LSF -
** pass the LENGTH of the LPC vector - this is the LPC
** order plus 1. Must pre-allocate lsfs to length+1
************************************************************/
void lpc2lsf(const EST_FVector &lpc, EST_FVector &lsf)
{
(void) lpc;
(void) lsf;
EST_error("LSF Code unfinished\n");
}
void lsf2lpc(const EST_FVector &lpc, EST_FVector &lsf)
{
(void) lpc;
(void) lsf;
EST_error("LSF Code unfinished\n");
}
void sig2lpc(const EST_FVector &sig, EST_FVector &acf,
EST_FVector &ref, EST_FVector &lpc)
{
int i, j;
float e, ci, sum;
int order = lpc.length() -1;
EST_FVector tmp(order);
int stableorder=-1;
if ((acf.length() != ref.length()) ||
(acf.length() != lpc.length()))
EST_error("sig2lpc: acf, ref are not of lpc's order");
//cerr << "sig2lpc order " << order << endl;
for (i = 0; i <= order; i++)
{
sum = 0.0;
for(j = 0; j < sig.length() - i; j++)
sum += sig.a_no_check(j) * sig.a_no_check(j + i);
acf.a_no_check(i) = sum;
}
// find lpc coefficients
e = acf.a_no_check(0);
lpc.a_no_check(0) = 1.0;
for (i = 1; i <= order; i++)
{
ci = 0.0;
for(j = 1; j < i; j++)
ci += lpc.a_no_check(j) * acf.a_no_check(i-j);
if (e == 0)
ref.a_no_check(i) = ci = 0.0;
else
ref.a_no_check(i) = ci = (acf.a_no_check(i) - ci) / e;
//Check stability of the recursion
if (absval(ci) < 1.000000)
{
lpc.a_no_check(i) = ci;
for (j = 1; j < i; j++)
tmp.a_no_check(j) = lpc.a_no_check(j) -
(ci * lpc.a_no_check(i-j));
for( j = 1; j < i; j++)
lpc.a_no_check(j) = tmp.a_no_check(j);
e = (1 - ci * ci) * e;
stableorder = i;
}
else break;
}
if (stableorder != order)
{
fprintf(stderr,
"warning:levinson instability, order restricted to %d\n",
stableorder);
for (; i <= order; i++)
lpc.a_no_check(i) = 0.0;
}
// normalisation for frame length
lpc.a_no_check(0) = e / sig.length();
}
void sig2pow(EST_FVector &frame, float &power)
{
power = 0.0;
for (int i = 0; i < frame.length(); i++)
power += pow(frame(i), float(2.0));
power /= frame.length();
}
void sig2rms(EST_FVector &frame, float &rms_energy)
{
sig2pow(frame, rms_energy);
rms_energy = sqrt(rms_energy);
}
float lpredict2(EST_FVector &adc, int wsize,
EST_FVector &acf, float *ref, float *lpc,
int order)
{
int i, j;
float e, ci, sum;
EST_TBuffer<float> tmp(order);
int stableorder=-1;
EST_FVector vref(order + 1), vlpc(order + 1);
for (i = 0; i <= order; i++)
{
sum = 0.0;
for (j = 0; j < wsize - i; j++)
sum += adc[j] * adc[j + i];
acf[i] = sum;
}
/* find lpc coefficients */
e = acf[0];
lpc[0] = 1.0;
for(i = 1; i <= order; i++) {
ci = 0.0;
for(j = 1; j < i; j++)
ci += lpc[j] * acf[i-j];
ref[i] = ci = (acf[i] - ci) / e;
//Check stability of the recursion
if (absval(ci) < 1.000000) {
lpc[i] = ci;
for(j = 1; j < i; j++)
tmp[j] = lpc[j] - ci * lpc[i-j];
for(j = 1; j < i; j++)
lpc[j] = tmp[j];
e = (1 - ci * ci) * e;
stableorder = i;
}
else break;
}
if (stableorder != order) {
fprintf(stderr,
"warning:levinson instability, order restricted to %d\n",
stableorder);
for (;i<=order;i++)
lpc[i]=0.0;
}
return(e);
}
void sig2fbank(const EST_FVector &sig,
EST_FVector &fbank_frame,
const float sample_rate,
const bool use_power_rather_than_energy,
const bool take_log)
{
EST_FVector fft_frame;
int i,fbank_order;
float Hz_per_fft_coeff;
// upper and lower limits of filter bank
// where the upper limit depends on the sampling frequency
// TO DO : add low/high pass filtering HERE
float mel_low=0;
float mel_high=Hz2Mel(sample_rate / 2);
// FFT this frame. FFT order will be computed by sig2fft
// FFT frame returned will be half length of actual FFT performed
sig2fft(sig,fft_frame,use_power_rather_than_energy);
// this is more easily understood as half the sampling
// frequency over half the fft order, but fft_frame_length()
// is already halved
Hz_per_fft_coeff = 0.5 * sample_rate / fft_frame.length();
fbank_order = fbank_frame.length();
// store the list of centre frequencies and lower and upper bounds of
// the triangular filters
EST_FVector mel_fbank_centre_frequencies(fbank_order+2);
mel_fbank_centre_frequencies[0]=mel_low;
for(i=1;i<=fbank_order;i++)
mel_fbank_centre_frequencies[i] = mel_low +
(float)(i) * (mel_high - mel_low) / (fbank_order+1);
mel_fbank_centre_frequencies[fbank_order+1]=mel_high;
// bin FFT in Mel filters
fft2fbank(fft_frame,
fbank_frame,
Hz_per_fft_coeff,
mel_fbank_centre_frequencies);
if(take_log)
for(i=0;i<fbank_frame.length();i++)
fbank_frame[i] = safe_log(fbank_frame[i]);
}
void sig2fft(const EST_FVector &sig,
EST_FVector &fft_vec,
const bool use_power_rather_than_energy)
{
int i,half_fft_order;
float real,imag;
float window_size = sig.length();
int fft_order = fft_vec.length();
// work out FFT order required
fft_order = 2;
while (window_size > fft_order)
fft_order *= 2;
fft_vec = sig;
// pad with zeros
fft_vec.resize(fft_order);
// in place FFT
(void)fastFFT(fft_vec);
// of course, we only need the lower half of the fft
half_fft_order = fft_order/2;
for(i=0;i<half_fft_order;i++)
{
real = fft_vec(i*2);
imag = fft_vec(i*2 + 1);
fft_vec[i] = real * real + imag * imag;
if(!use_power_rather_than_energy)
fft_vec[i] = sqrt(fft_vec(i));
}
// discard mirror image, retaining energy/power spectrum
fft_vec.resize(half_fft_order);
}
void fft2fbank(const EST_FVector &fft_frame,
EST_FVector &fbank_vec,
const float Hz_per_fft_coeff,
const EST_FVector &mel_fbank_frequencies)
{
// expects "half length" FFT - i.e. energy or power spectrum
// energy is magnitude; power is squared magnitude
// mel_fbank_frequencies is a vector of centre frequencies
// BUT : first element is lower bound of first filter
// last element is upper bound of final filter
// i.e. length = num filters + 2
int i,k;
float this_mel_centre,this_mel_low,this_mel_high;
EST_FVector filter;
int fft_index_start;
// check that fbank_vec and mel_fbank_frequencies lengths match
if(mel_fbank_frequencies.length() != fbank_vec.length() + 2)
{
EST_error("Filter centre frequencies length (%i) is not equal to fbank order (%i) plus 2\n",mel_fbank_frequencies.length(),
fbank_vec.length());
return;
}
// filters are computed on the fly
for(i=0;i<fbank_vec.length();i++)
{
// work out shape of the i'th filter
this_mel_low=mel_fbank_frequencies(i);
this_mel_centre=mel_fbank_frequencies(i+1);
this_mel_high=mel_fbank_frequencies(i+2);
make_mel_triangular_filter(this_mel_centre,this_mel_low,this_mel_high,
Hz_per_fft_coeff,fft_frame.length(),
fft_index_start,filter);
// do filtering
fbank_vec[i]=0.0;
for(k=0;k<filter.length();k++)
fbank_vec[i] += fft_frame(fft_index_start + k) * filter(k);
}
}
void fbank2melcep(const EST_FVector &fbank_vec,
EST_FVector &mfcc_vec,
const float liftering_parameter,
const bool include_c0)
{
// a cosine transform of the fbank output
// remember to pass LOG fbank params (energy or power)
int i,j,actual_mfcc_index;
float pi_i_over_N,cos_xform_order,const_factor;
float PI_over_liftering_parameter;
if(liftering_parameter != 0.0)
PI_over_liftering_parameter = PI / liftering_parameter;
else
PI_over_liftering_parameter = PI; // since sin(n.PI) == 0
// if we are not including cepstral coeff 0 (c0) then we need
// to do a cosine transform 1 longer than otherwise
cos_xform_order = include_c0 ? mfcc_vec.length() : mfcc_vec.length() + 1;
const_factor = sqrt(2 / (float)(fbank_vec.length()));
for(i=0;i<mfcc_vec.length();i++)
{
actual_mfcc_index = include_c0 ? i : i+1;
pi_i_over_N =
PI * (float)(actual_mfcc_index) / (float)(fbank_vec.length());
for(j=0;j<fbank_vec.length();j++)
// j + 0.5 is because we want (j+1) - 0.5
mfcc_vec[i] += fbank_vec(j) * cos(pi_i_over_N * ((float)j + 0.5));
mfcc_vec[i] *= const_factor;
// liftering
mfcc_vec[i] *= 1 + (0.5 * liftering_parameter
* sin(PI_over_liftering_parameter * (float)(actual_mfcc_index)));
}
}
void make_mel_triangular_filter(const float this_mel_centre,
const float this_mel_low,
const float this_mel_high,
const float Hz_per_fft_coeff,
const int half_fft_order,
int &fft_index_start,
EST_FVector &filter)
{
// makes a triangular (on a Mel scale) filter and creates
// a weight vector to apply to FFT coefficients
int i,filter_vector_length,fft_index_stop;
float rise_slope,fall_slope,this_mel;
// slopes are in units per Mel
// this is important - slope is linear in MEl domain, not Hz
rise_slope = 1/(this_mel_centre - this_mel_low);
fall_slope = 1/(this_mel_centre - this_mel_high);
// care with rounding - we want FFT indices **guaranteed**
// to be within filter so we get no negative filter gains
// (irint gives the _nearest_ integer)
// round up
if(this_mel_low == 0)
fft_index_start=0;
else
fft_index_start = irint(0.5 + (Mel2Hz(this_mel_low) / Hz_per_fft_coeff));
// round down
fft_index_stop = irint((Mel2Hz(this_mel_high) / Hz_per_fft_coeff) - 0.5);
if(fft_index_stop > half_fft_order-1)
fft_index_stop = half_fft_order-1;
filter_vector_length = fft_index_stop - fft_index_start + 1;
filter.resize(filter_vector_length);
for(i=0;i<filter_vector_length;i++)
{
this_mel = Hz2Mel( (i + fft_index_start) * Hz_per_fft_coeff );
if(this_mel <= this_mel_centre)
{
filter[i] = rise_slope * (this_mel - this_mel_low);
}
else
{
filter[i] = 1 + fall_slope * (this_mel - this_mel_centre);
}
}
}
float Hz2Mel(float frequency_in_Hertz)
{
return 1127 * log(1 + frequency_in_Hertz/700.0);
}
float Mel2Hz(float frequency_in_Mel)
{
return (exp(frequency_in_Mel / 1127) - 1) * 700;
}
EST_FVector
Definition EST_FMatrix.h:114
EST_String
Definition EST_String.h:70
EST_TBuffer
Definition EST_TBuffer.h:86
EST_TSimpleVector::resize
void resize(int n, int set=1)
resize vector
Definition EST_TSimpleVector.cc:70
EST_TVector::length
INLINE int length() const
number of items in vector.
Definition EST_TVector.h:252
EST_TVector::a_no_check
INLINE const T & a_no_check(int n) const
read-only const access operator: without bounds checking
Definition EST_TVector.h:257