I have encountered a problem with overflow. It can be seen on the picture with the flat tops.

It should have looked like the picture with the sharp tops.

The contexts is this: I am doing speechrecognition and need to calculate cepstrum coefficients of input speech:

Matlab code: Cepstrum = real(ifft(log(abs(fft(x)))));

So far I have made...

1) fft with the rfft_fr16 function,

2) abs with the cabs_fr16 function

3) Then i needed to do the log. BUT i could only find a floating point version of this one? Does a fract16 version exist?

i could´nt find one - therefore i used the fr16_to_float function and used the logf function.

4) Everything untill this point has given quite accurate results, but then I had to convert from float_to_fr16 to be able to use the

ifft_fr16 function. Here is my problem.

5) Because the floating point values are higher (on the tops and lower in the valleys) than +1 and -1 it overflows! I know that I probably need to divide the log(fft) with some factor and then convert, but I also want to have a nice precision when doing the ifft afterwards.

Should I loop through the log(fft) output and find the highest value and divide everything with this value, so that it is normalized to [-1 1]?

Then I could use static scaling on the ifft afterwards! Is there a more efficient way to do the normalization?

By the way: This discussion is only posted here.

Thank you

Benjamin

Here is a snippet of my code:

#define NUMPOINTS 256

complex_fract16 w[NUMPOINTS]; //twiddle sequence

complex_fract16 fftout[NUMPOINTS];

int block_exponent;

int n = NUMPOINTS;

float fft_scaled[NUMPOINTS];

float log_of_fft[NUMPOINTS];

//Make w twiddle table for FFT once.

twidfftrad2_fr16(w, n);

void cepstrum_test(const fract16 samples256[]) {

//cepstrum = real(ifft(log(abs(fft(x)))));

int i = 0;

//FFT

rfft_fr16(samples256, fftout, w, 1, NUMPOINTS, &block_exponent, 2);

float scaling_fft = ldexpf(1.0F,block_exponent);

for(i=0;i<NUMPOINTS;i++) //Tak abs(fft) and scale fft and take logarithm and convert back to fr16:

{

fft_scaled[i] = fr16_to_float(cabs_fr16(fftout[i]))*scaling_fft; // WOULD IT BE MORE EFFICIENT TO SHIFT THE FR16 INSTEAD OF MULTIPLYING // WITH THE SCALING FACTOR IN THE FLOATING POINT DOMAIN?

log_of_fft[i] = logf(fft_scaled[i]);

fftout[i].re = float_to_fr16(log_of_fft[i]); //HERE IS THE PROBLEM.

fftout[i].im = 0;

}

}

I found a possible soliution to my problem!!!

I think it is a fairly good implementation of the matlab cepstrum calculation: cepstrum = real(ifft(log(abs(fft(x)))));

"CORRECTION: I had wrote that it gave an error less than 1 %:

This was only true for one special welldefined signal consisting of two sinusoids. It gave a much greater error when trying a real speech signal and therefore I decided to move to floating point or fract32 instead of fract16. Floating point of course gave precise results. I have not tried the fract32 yet. If one is looking for the best FFT calculation with either fract16, fract32 or floating point - It depends on the degree of accuracy required, as I see it."

BUT I would still like to get feedback on it, because I´m not sure this

is the most efficient way to handle my problem. Thank you for your help.

Best regards

Benjamin

Here´s the solution I made:

void cepstrum_test(const fract16 samples256[]) {

//cepstrum = real(ifft(log(abs(fft(x)))));

int i = 0;

//FFT

rfft_fr16(samples256, fftout, w, 1, NUMPOINTS, &block_exponent, 1);

float scaling_fft = ldexpf(1.0F,block_exponent);

for(i=0;i<NUMPOINTS;i++) //Take absolute value, convert to float, scale it and take logarithm.

{

fft_scaled[i] = fr16_to_float(cabs_fr16(fftout[i]))*scaling_fft;

log_of_fft[i] = logf(fft_scaled[i]);

}

// There is a need to normalize the floating point values before converting to

// fract16 or else there will be overflow. The fract16 has the range [-1; 0.99997]. Thats almost from -1 to 1.

// This includes some steps:

// (1) Finding the maximum or minimum value of the vector.

// (2) This value will be converted to the nearest power of two (To avoid dividing by a skew number like 3.971)

// (3) Normalize by dividing with the value from (2), convert to fract16 and continue to ifft.

max_value = vecmaxf(log_of_fft,NUMPOINTS);

min_value = vecminf(log_of_fft,NUMPOINTS);

unsigned int vv;

unsigned int next_pow2;

if (max_value > fabsf(min_value)) //Check wheither it is min or max which is greatest.

{

vv = ceilf(max_value); //cast from floating point to unsigned int.

next_pow2 = power_of_2(vv); //Find next power of two.

}

else

{

vv = ceilf(fabsf(min_value)); //Check if cast is ok

next_pow2 = power_of_2(vv);

}

for(i=0;i<NUMPOINTS;i++) //Normalize and convert to fract16

{

//log_of_fft[i]/max_value;

fftout[i].re = float_to_fr16(log_of_fft[i]/next_pow2);

fftout[i].im = 0;

//For plotting purposes: new_log_of_fft[i] = fr16_to_float(fftout[i].re);

}

ifft_fr16(fftout, fftout, w, 1, NUMPOINTS, &block_exponent2, 1);

for(i=0;i<NUMPOINTS;i++)

{

inversefftout[i] = fr16_to_float(fftout[i].re)*next_pow2;

}

for(i=0;i<14;i++) //Take real of the first 14 cepstrum coeff without c[0].

{

featureVec[i] = inversefftout[i+1];

}

}

// Thanks to: http://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2:

unsigned int power_of_2(unsigned int v)

{

//unsigned int v; // compute the next highest power of 2 of 32-bit v

v--;

v |= v >> 1;

v |= v >> 2;

v |= v >> 4;

v |= v >> 8;

v |= v >> 16;

v++;

return v;

}