//Calculate LPC coefficients coefficients based on autocorrelation
void corref(int ip,float cor[],float alf[])

{
int i,j;
float refi,alfsave,resid;
alf[1]=(-cor[1]);
alf[0]=1;
resid=1-cor[1]*cor[1];
for(i=2;i<=ip;i++)
{
refi=cor[i];
for(j=1;j<i;j++)
refi+=alf[j]*cor[i-j];
refi/=resid;
alf[i]=(-refi);
for(j=1;2*j<=i;j++)
{
alfsave=alf[j];
alf[j]=alfsave-refi*alf[i-j];
if(2*j!=i)
alf[i-j]-=refi*alfsave;
}
}
}

//Calculate LPC cepstrum based on alf
void alfcep(int ip,float alf[],float cep[],int n)

{
int m,i;
float ss;
m=0;
goto start;
next:
m+=1;
ss=0;
i=0;
goto step3;
step2:
ss-=alf[i]*cep[m-i];
step3:
i+=1;
if(i>ip) goto step4;
if(i<m) goto step2;
ss-=alf[m]*m;
step4:
cep[m]=ss;
start:
if(m<n) goto next;
for(m=2;m<=n;m++)
cep[m]/=m;
}

//Calculate mel-warped LPC cepstrum
void cepmel(float cep[],float mel[])
{
int k,n;
float u0[ORDER+2],u1[ORDER+2];
float alpha=ALPHA;
//Calculation mel[0]*
u0[ORDER+1]=0;
for(n=ORDER;n>=0;--n)
u0[n]=cep[n]+alpha*u0[n+1];
//mel[0]=u0[0];
//Calculation mel[1]*
u1[ORDER+1]=0;
for(n=ORDER;n>=0;--n)
u1[n]=((1-alpha*alpha)*u0[n+1]+alpha*u1[n+1]);
mel[1]=u1[0];
//Calculation mel[2],mel[3],......,mel[ORDER] *
for(k=2;k<=ORDER;++k)
{
for(n=0;n<=ORDER+1;++n)
u0[n]=u1[n];
for(n=ORDER;n>=0;--n)
u1[n]=u0[n+1]+alpha*(u1[n+1]-u0[n]);
mel[k]=u1[0];
}
}

main()
{
int rad[][]; //端点检测后的语音数据
float cor[ORDER+1]; //autocorrelation coefficients ORDER＝14
float alf[ORDER+1]; //LPC
float cep[ORDER+1]; //LPC cepstrum
float mel[ORDER+1];
float mel_data[length][ORDER];
//mel-warped LPC cepstrum，length为端点检测后的每个语音数据的帧长
int mel_v=0;

//Auto-Correlation Calculation
correl(rad[count],cor,W_SIZE,ORDER);
//LPC Analysis Program
corref(ORDER,cor,alf);
//LPC Cepstrum Coefficient Calculation
alfcep(ORDER,alf,cep,ORDER);
alf[0]=cep[0]=cor[0]; //cor[0] is power
//MEL Cepstrum Coefficient Calculation
cepmel(cep,mel);
//printf("\n Frame %d :",count);

for(i=0;i<ORDER;i++)
{
mel_data[mel_v][i]=mel[i+1];
//printf("\nFrame %d: Mel[%d][%d]= %g",count,mel_v,i,mel_data[mel_v][i]);
}
}