This is fitlpc.c in view mode; [Download] [Up]
/* fitlpc routine by Perry R. Cook
Stanford CCRMA 1991.
This performs linear predictive coding
to a file. Read lpcInstructions for
details.
*/
#import <stdio.h>
#import <fcntl.h>
#import <string.h>
#import <math.h>
#define MAX_BLOCK 4096
#define MAX_ORDER 100
float autocorr(int size,float *data,float *result)
{
int i,j,k;
float temp,norm;
for (i=0;i<size/2;i++) {
result[i] = 0.0;
for (j=0;j<size-i-1;j++) {
result[i] += data[i+j] * data[j];
}
}
temp = result[0];
j = size*0.02;
while (result[j]<temp && j < size/2) {
temp = result[j];
j += 1;
}
temp = 0.0;
for (i=j;i<size*0.5;i++) {
if (result[i]>temp) {
j = i;
temp = result[i];
}
}
norm = 1.0 / size;
k = size/2;
for (i=0;i<size/2;i++) result[i] *= (k - i) * norm;
if ((result[j] / result[0]) < 0.4) j = 0;
if (j > size/4) j = 0;
return (float) j;
}
int minvert(int size,float mat[][MAX_ORDER])
{
int item,row,col,rank,t2;
float temp,res[MAX_ORDER][MAX_ORDER];
int ok,zerorow;
for (row=1;row<=size;row++) {
for (col=1;col<=size;col++) {
// fprintf(stdout," %f ",mat[row][col]);
if (row==col)
res[row][col] = 1.0;
else
res[row][col] = 0.0;
}
// fprintf(stdout,"\n");
}
for (item=1;item<=size;item++) {
if (mat[item][item]==0) {
for (row=item;row<=size;row++) {
for (col=1;col<=size;col++) {
mat[item][col] = mat[item][col] + mat[row][col];
res[item][col] = res[item][col] + res[row][col];
}
}
}
for (row=item;row<=size;row++) {
temp=mat[row][item];
if (temp!=0) {
for (col=1;col<=size;col++) {
mat[row][col] = mat[row][col] / temp;
res[row][col] = res[row][col] / temp;
}
}
}
if (item!=size) {
for (row=item+1;row<=size;row++) {
temp=mat[row][item];
if (temp!=0) {
for (col=1;col<=size;col++) {
mat[row][col] = mat[row][col] - mat[item][col];
res[row][col] = res[row][col] - res[item][col];
}
}
}
}
}
for (item=2;item<=size;item++) {
for (row=1;row<item;row++) {
temp = mat[row][item];
for (col=1;col<=size;col++) {
mat[row][col] = mat[row][col] - temp * mat[item][col];
res[row][col] = res[row][col] - temp * res[item][col];
}
}
}
/* ok = TRUE;
rank = 0;
for (row=1;row<=size;row++) {
zerorow = TRUE;
for (col=1;col<=size;col++) {
if (mat[row][col]!=0) zerorow = FALSE;
t2 = (mat[row][col] + 0.5);
if (row==col&&t2!=1) ok = FALSE;
t2 = fabs(mat[row][col]*100.0);
if (row!=col&&t2!=0) ok = FALSE;
}
if (!zerorow) rank += 1;
}
if (!ok) {
fprintf(stdout,"Matrix Not Invertible\n");
fprintf(stdout,"Rank is Only %i of %i\n",rank,size);
} */
for (row=1;row<=size;row++) {
for (col=1;col<=size;col++) {
mat[row][col] = res[row][col];
}
}
return rank;
}
float lpc_from_data(int order, int size, float *data, float *coeffs)
{
float r_mat[MAX_ORDER][MAX_ORDER];
int i,j;
float pitch;
float corr[MAX_BLOCK];
pitch = autocorr(size, data,corr);
for (i=1;i<order;i++) {
for (j=1;j<order;j++) r_mat[i][j] = corr[abs(i-j)];
}
minvert(order-1,r_mat);
for (i=0;i<order-1;i++) {
coeffs[i] = 0.0;
for (j=0;j<order-1;j++) {
coeffs[i] += r_mat[i+1][j+1] * corr[1+j];
}
}
return pitch;
}
float predict(int order,int hop_size,float *data,float *coeffs, int resFile)
{
int i,j,n_write;
float power=0.0,error,tmp;
float Zs[MAX_ORDER];
for (i=0;i<order;i++) Zs[i] = data[order-i-1];
for (i=order;i<=hop_size+order;i++) {
tmp = 0.0;
for (j=0;j<order;j++) tmp += Zs[j]*coeffs[j];
for (j=order-1;j>0;j--) Zs[j] = Zs[j-1];
Zs[0] = data[i];
// }
// for (i=0;i<hop_size;i++) {
// tmp = coeffs[0] + data[i+order-1];
// for (j=1;j<order;j++) tmp += coeffs[j] * data[i+order-1-j];
error = data[i] - tmp;
n_write = write(resFile,&error,4);
power += error * error;
}
return sqrt(power) / hop_size;
}
int check_stable(int order,float *coeffs)
{
int i,j;
float output,Zs[MAX_ORDER];
for (i=0;i<order;i++) Zs[i] = 0.0;
Zs[0] = 1.0;
for (i=0;i<2*order;i++) {
output = 0.0;
for (j=0;j<order;j++) output += Zs[j]*coeffs[j];
for (j=order-1;j>0;j--) Zs[j] = Zs[j-1];
Zs[0] = output;
}
j = 1;
if (output>1.0) j = 0;
return j;
}
main(ac,av)
int ac;
char *av[];
{
char *file_name_in,*file_name_out,residual_file[200];
float data[MAX_BLOCK],out_data[MAX_ORDER];
int fd,fd2,fdr,n_read,n_write;
int i,time=0,block_size,hop_size,order;
float pitch,power;
if (ac==6) {
order = atoi(av[1]);
block_size = atoi(av[2]);
hop_size = atoi(av[3]);
if (hop_size>block_size) {
hop_size = block_size - order;
fprintf(stderr,"Setting hop size to block size - order.\n");
}
file_name_in = av[4];
file_name_out = av[5];
fd = open(file_name_in,0);
if (fd!=-1) {
if (!strcmp(file_name_in,file_name_out)) {
fprintf(stderr,"Input and output files can't have same name!!\n");
exit(0);
}
fd2 = creat(file_name_out, 0666);
n_write = write(fd2,&order,4);
n_write = write(fd2,&hop_size,4);
n_read = read(fd,data,block_size*4);
n_read = 10000;
strcpy(residual_file,file_name_in);
strcat(residual_file,".residual");
fdr = creat(residual_file, 0666);
while (n_read>0) {
pitch = lpc_from_data(order+1,block_size,data,out_data);
if (!check_stable(order,out_data)) {
fprintf(stderr,"This set is not stable!!\n");
// stabilize(order,out_data);
}
power = predict(order,hop_size,data,out_data,fdr);
time+= hop_size;
printf("%i %f %f\n",time,pitch,power);
for (i=0;i<block_size-hop_size;i++) data[i] = data[hop_size+i];
n_read = read(fd,&data[block_size-hop_size],hop_size*4);
n_write = write(fd2,&pitch,4);
n_write = write(fd2,&power,4);
n_write = write(fd2,&out_data,order*4);
}
close(fd);
close(fd2);
}
else
printf("I couldn't find your input file!!!\n");
}
else
printf("Format is fitlpc order blocksize hopsize infile outfile.lpc\n");
return;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.