a tip about Nan problem

[YangangCao]

Hi All

When we calculate lpc, auto correlation series is needed. However, if the element in auto correlation series all are zero, nan will appears. (Test code is in offered in the end)
PercepNet calculates pitch using auto correlation. In function "_celt_autocorr" in celp_lpc.cpp, if the absolute value of ac[0]、ac[1]、ac[2] and ac[3] all are less than 1e-10, Nan maybe appear. (I don't know why Nan isn't must appear, but I think it's must be wrong feature in this frame)
I try to add a decision about it, if less than 1e-10, just recursive return -1 and add continue to calculate the next frame. By do this, I avoid Nan.

// This is my test code, just copy into *.c and gcc *.c to compile
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


static int autoCorrelation(const float *in, float *out, size_t sampleCount, size_t p)
{
    if ((in == NULL) || (out == NULL))
        return -1;

    size_t times = sampleCount;
    size_t k = 0;
    for (size_t i = 0; i < p + 1; ++i)
    {
        float result = .0f;
        for (size_t j = 0; j < times; ++j)
        {
            result += in[j + k] * in[j];
        }
        k++;
        times--;
        out[i] = result;
    }
    return 0;
}


// r: length is p+1, r(0),...r(p)
// a: length is p+1, a(0),...a(p)
static float levinson (const float *r, size_t p, float *a)
{
    if ((r == NULL) || (a == NULL))
        return -1;

    size_t i, k;
    float q,v=0;
    float *b = (float *)malloc(sizeof(float)*(p+1));

    a[0] = 1.f;
    a[1] = -r[1] / r[0];
    v = (1 - a[1] * a[1]) * r[0];

    for (k = 2; k <= p; ++k)
    {
        q = 0.f;
        for (i = 0; i < k; ++i)
            q += a[i] * r[k - i];
        a[k] = -q / v;
        for (i = 1; i < k; ++i)
            b[i] = a[i] + a[k] * a[k - i];
        for (i = 1; i < k; ++i)
            a[i] = b[i];
        v = (1 - a[k] * a[k]) * v;
    }
    float gain = 0.f;
    for(size_t i = 0; i < p+1; i++)
        gain += r[i]*a[i];
    free(b);
    return gain;

}

// out: the length is p+1
float lpc(const float * in, float * out, size_t sampleCount, size_t p)
{
    if ((in == NULL) || (out == NULL))
        return -1;
    float *autoCorr = (float *)malloc(sizeof(float)*(p+1));
    autoCorrelation(in, autoCorr, sampleCount, p);
    float v = levinson(autoCorr, p, out)/ sampleCount;
    free(autoCorr);
    return v; // v is error power
}

int main(int argc, char *argv[] )
{
    const int len = 100;
    const int p = 10;
    float in[len] = {0};
    float out[p+1]={0};
    float c = lpc(&in[0], &out[0],len,p);
    for (int i = 0; i < p+1; ++i)
    {
        if isnan(out[i])
            printf("out[%d] is nan\n",i);
        else
            printf("out[%d]:%f\n",i,out[i]);
    }

    float a = 1.f/0.f;
    if (isinf(a))
    {
        printf("1.f/0.f is inf\n");
    }

    if (isnan(a/a))
    {
        printf("inf/inf is nan\n");
    }

    if (isnan(0.f/0.f))
    {
        printf("0.f/0.f is nan\n");
    }
    return 0;
}

[jzi040941]

Hi, @YangangCao

Thank you for let us know this solution.
I've try to understand what you are saying but not fully understand how to add exceptional decision ( return -1, & continue to next frame)
I want to ask about this

I was getting in to deeply how pitch estimation is working in PercepNet
pitch.cpp : pitch_downsample function calls _celt_autocorr
according to your tip when ac[0]~ac[4] is very small like 1e-10 almost zero
So do I have to add decision inside pitch_downsample function after calling _celt_autocorr like below?

//pitch.cpp : L183
   _celt_autocorr(x_lp, ac, NULL, 0,
                  4, len>>1);

  if(ac[0]<1e-10 && ac[1]<1e-10 && ac[2]< 1e-10 && ac[3] < 1e-10)
     contiune; //how can I do continue ..?

I don't think my understanding is correct. could you tell me which part you added for exceptional decision on code?
Thanks!

[YangangCao]

Hi @jzi040941
Maybe my words made you confused, I will explain it with code. recursive return -1 in above mean that all function in call stack return -1 until to the while (1) in main function in denoise.cpp, my code like this:

 if (fabsf(ac[0]) < 1e-10 && fabsf(ac[1]) < 1e-10 && fabsf(ac[2]) < 1e-10 && fabsf(ac[3]) < 1e-10)
 {
    return -1;
 } // _celt_autocorr

if(_celt_autocorr(x_lp, ac, NULL, 0,  4, len>>1)==-1)
{
   return -1;
} // pitch_downsample

if (pitch_downsample(pre, pitch_buf, PITCH_BUF_SIZE, 1) == -1)  
{
    return -1;
} // compute_frame_features

int silence1 = compute_frame_features(noisy, Y, Phat/*not use*/, Ey, Ephat/*not use*/, Ephaty, features, xn);
int silence2 = compute_frame_features(st, X, P, Ex, Ep, Exp, features, x);
if (silence==-1 || silence2==-1)
{
  continue;
} // while(1) in main

As to the pitch calculate, I recommend you to read this doc https://jmvalin.ca/papers/aes135_opus_silk.pdf#:~:text=Vos%20et%20al.%20Voice%20Coding%20with%20Opus%20carded,residual%20taking%20into%20account%20the%20LSF%20quantization%20e%0Bects.

[jzi040941]

Thanks for your kind explanation and reference doc
I will try it soon