/*Introducing Speech and Language Processing
http://www.islp.org.uk*/


//lpc_spectrum.c
===============================================
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "slputils.c"
#include "four1.c"

/* LPC_SPECTRUM.C                                                                */
/* Reads a signal from a disk file into a variable, x_in.                        */
/* Computes the LPC coefficients of sample_number, pads it with zeroes to 512        */
/* and writes the resulting signal to outfile.                                */

int main(int argc, char *argv[])
{
        short int *x_in, *signal_in();
        char *infile, *sampleno, *endcp;
        int sample, i, *length=0;
        float frame[81], *xms, d[17], *x_in_f, logpsd[512], x_frame[512];
        void memcof();
        float fftdata[1024];                        /* 512 complex data points        */
        void four1();

        if (argc != 3) {
                printf("usage: lpc_spectrum input_file sample_number [> output_file]\n");
                exit(1);
        }

        infile = argv[1];
        sampleno = argv[2];

        x_in = signal_in(infile,length);
        sample = (int) strtoul(sampleno,&endcp,10);

/* Make a floating-point version of x_in, called x_in_f                        */
        x_in_f = (float *) calloc(*length,sizeof(float));
        if (!x_in_f) {
                printf("Can't allocate space for x_in_f.\n");
                exit(1);
        }

        for (i=0;i<*length;i++) {
                x_in_f[i] = ((float) x_in[i]);
        }

        for (i=0;i<=79;i++) {
                frame[i+1] = x_in_f[sample+i];
        }

        memcof(frame,80,10,&xms,d);

        for(i=0;i<251;i++) x_frame[i] = 0; 
        for(i=0;i<10;i++) x_frame[i+248] = 15000.0*d[i+1]; 
        for(i=261;i<512;i++) x_frame[i] = 0;

        for (i=0;i<=511;i++) {
                fftdata[2*i] = (float) x_frame[i];
                fftdata[2*i+1] = 0.0;
        }

        four1(fftdata-1,512,1);

/* In the log power spectral density, magnitudes are in dB                */
/* in steps of SampleRate/512 Hz (15.625Hz at 8000 samples/s).        */

   printf("f (Hz)\tAmplitude (dB)\n");

   for (i=0;i<=256;i++) {
      logpsd[i] = fftdata[2*i];
      logpsd[i] *= logpsd[i];
      logpsd[i] += SQR(fftdata[2*i+1]);
      logpsd[i] = -10*log10(logpsd[i]);
      printf("%d\t%.2f\n",(int) (i*15.625),logpsd[i]);
   }
   return 0;
}


four1.c
===============================================
#include <math.h>
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

/* four1.c                                                                */
/* Discrete Fourier Transform                                                */
/* This routine is from the book Numerical Recipes in C.                */
/* (Cambridge University Press), Copyright (C)  1987-1992 by            */
/* Numerical Recipes Software.  Used by permission.  Use of this        */
/* routine other than as an integral part of the software collection    */
/* provided in accompaniment with the book Introducing Speech and       */
/* Language Processing (Cambridge University Press) requires an         */
/* additional license from Numerical Recipes Software ([url]www.nr.com).[/url]     */
/* Further distribution in any form is prohibited.                      */

void four1(float data[], unsigned long nn, int isign)
/* Replaces data[1..2*nn] by its discrete Fourier transform, if isign        */
/* is input as 1; or replaces data[1..2*nn] by nn times its inverse        */
/* discrete Fourier transform, if isign is input as -1. data is a         */
/* complex array of length nn or, equivalently, a real array of length        */
/* 2**nn. nn MUST be an integer power of 2 (this is not checked for!).        */
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        float tempr,tempi;

        n=nn << 1;
        j=1;
        for(i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        
        mmax=2;
        while(n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318503717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for(m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1]+= tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}



slputils.c
===============================================
#include <stdio.h>
#include <stdlib.h>
#include "nrutil.h"

/* SLPUTILS.C                                                                        */

short int *signal_in(char *infile, int *length)

/* Reads a signal from a disk file into a variable, x_in.                        */
/* The length of the signal (in shorts) is returned in the variable *length.        */

{
        FILE *fid;
        int j, status;
        short int input, *x_in;

/* Open the file */

        fid = fopen(infile,"rb");
        if(fid == NULL) {
                printf("Error opening %s\n",infile);
                exit(1);
        }

/* First read in shorts, until the end of the file is reached, just to measure the file */
        j = 0;
        while(!feof(fid)) {
                status = fread(&input,sizeof(short int),1,fid);
                if(status != 1) {
                        break;
                }
                j++;
        }
        (*length) = j;
        fclose(fid);

/* Allocate memory, and read in the file again        */

        x_in = (short int *) calloc(*length,sizeof(short int));

/* Open the file again */

        fid = fopen(infile,"rb");
        if(fid == NULL) {
                printf("Error re-opening %s\n",infile);
                exit(1);
        }

        for(j=0;j<*length;j++) {
                status = fread(&input,sizeof(short int),1,fid);
                x_in[j] = input;
                if(status != 1) {
                        printf("Read error at sample %d\n",j);
                        exit(1);
                }
        }
        *length = j;
        fclose(fid);
        return(x_in);
}

void signal_out(int *length, short int *x_out, char *outfile)
{
        FILE *fid;
        int status;

        fid = fopen(outfile,"wb");
        if(fid == NULL) {
                printf("Can't open %s\n",outfile);
                printf("It may be in use by another application.\n");
                exit(1);
        }
        status = fwrite(x_out,sizeof(short int),*length,fid);
        if(status < *length) {
                printf("Error writing %s\n",outfile);        
                exit(1);
        }
        fclose(fid);
}



nrutil.h
===============================================
/* nrutil.h                                                                */

/* From W. H.Press, S. A. Teukolsky, W. T. Vettering and B. P. Flannery */
/* (1992) Numerical Recipes in C: The Art of Scientific Computing        */
/* (2nd edition). Cambridge University Press.                                */

/* Non-copyright                                                         */

/* pp. 941-2: This file only contains selected portions, however.        */

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))

void nrerror(char error_text[]);
float *vector(long nl, long nh);
void free_vector(float *v, long nl, long nh);