tabledesign.c 6.73 KB
#include "tabledesign.h"

char usage[] = "[-o order -s bits -t thresh -i refine_iter -f frame_size] aifcfile";

main(int argc, char *argv[])

{
    extern char	*optarg;
    extern int	optind, opterr;	
    int 	c;

    char *progname = argv[0];
    
    double thresh = DEFAULT_THRESH;
    int order = DEFAULT_ORDER;
    int table_size = DEFAULT_SIZE;
    int refine_iter = DEFAULT_REFINE;
    int frame_size = DEFAULT_FSIZE;

    short *in_buffer;
    double *ac, *ref, e2, **a, dist;
    double **codebook, **training;
    double *dir;
    
    int nframes = 0;
    int i, j, d;
    int actual_size, n_entries;
    int *indx;
    int overflow = 0;
    
    double srate;
    long ntracks, nsamples, chans, afsamform, afsamwidth;
    
    AFfilehandle afinfile;
    AFfilesetup afsetup;

    if (argc<2){
        fprintf(stderr,"%s %s\n",progname,usage);
        exit(1);
    }

    while ((c=getopt(argc,argv,"o:s:t:i:f:"))!=-1){
        switch (c) {
            case 'o':
                if (sscanf(optarg,"%d",&order)!=1)
                    order = DEFAULT_ORDER;
                break;
            case 's':
                if (sscanf(optarg,"%d",&table_size)!=1)
                    table_size = DEFAULT_SIZE;
                break;
            case 'f':
                if (sscanf(optarg,"%d",&frame_size)!=1)
                    frame_size = DEFAULT_FSIZE;
                break;
            case 'i':
                if (sscanf(optarg,"%d",&refine_iter)!=1)
                    refine_iter = DEFAULT_REFINE;
                break;
            case 't':
                if (sscanf(optarg,"%lf",&thresh)!=1)
                    thresh = DEFAULT_THRESH;
                break;
        }
    }
    argv += optind-1;
    
    if ((afinfile = AFopenfile(argv[1], "r", NULL)) == AF_NULL_FILEHANDLE){
        fprintf(stderr,"%s: input AIFC file [%s] could not be opened.\n",progname, argv[1]);
        exit(1);
    }

    /*
     * Check for a valid AIFC file
     */
    if ((chans = AFgetchannels(afinfile, AF_DEFAULT_TRACK))!=1){
        fprintf(stderr,"%s: file [%s] contains %d channels, only 1 channel supported.\n",
                progname, argv[1], chans);
        exit(1);
    }

    if ((ntracks = AFgettrackids(afinfile, NULL)) != 1){
        fprintf(stderr,"%s: file [%s] contains %d tracks, only 1 track supported.\n",
                progname, argv[1], ntracks);
        exit(1);
    }

    AFgetsampfmt(afinfile, AF_DEFAULT_TRACK, &afsamform, &afsamwidth);
    if (afsamwidth != 16){
        fprintf(stderr,"%s: file [%s] contains %d bit samples, only 16 bit samples supported.\n",
                progname, argv[1], afsamwidth);
        exit(1);
    }

    /*
     * Initialize codebook storage
     */
    codebook = (double **) malloc((1<<table_size)*sizeof(double *));
    for (i=0; i<(1<<table_size); i++)
        codebook[i] = (double *) malloc((order+1)*sizeof(double));
    
    /*
     * Splitting direction
     */
    dir = (double *) malloc((order+1)*sizeof(double));
    
    in_buffer = (short *) malloc(2*frame_size*sizeof(short));
    for (i=0; i<2*frame_size; i++)
        in_buffer[i] = 0;
    ac = (double *) malloc((order+1)*sizeof(double));
    ref = (double *) malloc((order+1)*sizeof(double));
    
    /* For matrix method */
    a = (double **) malloc((order+1)*sizeof(double *));
    for (i=0; i<=order; i++)
        *(a+i) = (double *) malloc((order+1)*sizeof(double));
    indx = (int *) malloc((order+1)*sizeof(int));
    
    nsamples = AFgetframecnt(afinfile, AF_DEFAULT_TRACK);
    srate = AFgetrate(afinfile, AF_DEFAULT_TRACK);

    /*
     * Reserve storage for the training data
     */
    training = (double **) malloc(nsamples*sizeof(double *));
    nframes = 0;
    while(AFreadframes(afinfile, AF_DEFAULT_TRACK, in_buffer+frame_size, frame_size) == frame_size){
        acvect(in_buffer+frame_size, order, frame_size, ac);
        if (fabs(ac[0]) > thresh){

            acmat(in_buffer+frame_size, order, frame_size, a);

            /*
             * Lower-upper decomposition
             */
            if (lud(a, order, indx, &d)==0){
                
                /*
                 * Put solution in ac
                 */
                lubksb(a, order, indx, ac);
                ac[0] = 1.0;

                /*
                 * Convert to reflection coefficients - reject unstable vectors
                 */
                if (!kfroma(ac, ref, order)){

                    /*
                     * The training data is stored as tap values
                     */
                    training[nframes] = (double *) malloc((order+1)*sizeof(double));

                    training[nframes][0] = 1.0;
                    for (i=1; i<=order; i++){
                        /*
                         * Stabilize the filter
                         */
                        if (ref[i]>=1.0)
                            ref[i] = 1.0 - TINY;
                        if (ref[i]<=-1.0)
                            ref[i] = -1.0 + TINY;
                    }
                    afromk(ref, training[nframes], order);
                    nframes++;
                }
            }
        }
        for (i=0; i<frame_size; i++)
            in_buffer[i] = in_buffer[i+frame_size];
    }

    /*
     * To start things off find the average auto-correlation over
     * the complete data set.
     */
    ac[0] = 1.0;
    for (j=1; j<=order; j++)
        ac[j] = 0;
    for (i=0; i<nframes; i++){
        rfroma(training[i], order, codebook[0]);
        for (j=1; j<=order; j++)
            ac[j] += codebook[0][j];
    }
    for (j=1; j<=order; j++)
        ac[j] = ac[j]/nframes;

    /*
     * The average model
     */
    durbin(ac, order, ref, codebook[0], &e2);
    /*
     * Stabilize - could put this in durbin
     */
    for (j=1; j<=order; j++){
        if (ref[j]>=1.0)
            ref[j] = 1.0 - TINY;
        if (ref[j]<=-1.0)
            ref[j] = -1.0 + TINY;
    }
    afromk(ref, codebook[0], order);
        
    actual_size = 0;
    while (actual_size<table_size){

        n_entries = 1<<actual_size;
        
        /*
         * Split each codebook template into
         * two - the original and a shifted version
         */
        for (i=0; i<=order; i++)
            dir[i] = 0;
        dir[order-1] = -1.0;
        
        split(codebook,dir,order,n_entries,0.01);

        /*
         * Iterative refinement of templates
         */
        actual_size++;
        refine(codebook, order, 1<<actual_size, training, nframes, refine_iter, 0);
    }
    n_entries = 1<<actual_size;
    
    /*
     * Let's see what it looks like
     */
    fprintf(stdout,"%d\n%d\n",order,n_entries);
    for (i=0; i<n_entries; i++)
        overflow += print_entry(stdout, codebook[i], order);
    
    /*
     *
     */
    if (overflow>0)
        fprintf(stderr,"There was overflow - check the table\n");
}