lpctest.c

#include <sp/spBaseLib.h>
#include <sp/spLib.h>
#include <sp/spWave.h>
#include <sp/spInputPlugin.h>
 
#include <sp/fft.h>
#include <sp/lpc.h>
 
#include <sp/spMain.h>
 
#define LPCTEST_BUF_LENGTH 2048
 
long buf_length = LPCTEST_BUF_LENGTH;
 
long order;
long frame_start;
long frame_length;
 
int debug_level;
spBool help_flag;
 
static spOptions options;
static spOption option[] = {
    {"-order", NULL, "analysis order", NULL,
         SP_TYPE_LONG, &order, "10"},
    {"-st", "-start", "analysis start point", NULL,
         SP_TYPE_LONG, &frame_start, "0"},
    {"-frame", NULL, "frame length", NULL,
         SP_TYPE_LONG, &frame_length, "256"},
    {"-debug", NULL, "debug level", NULL,
         SP_TYPE_INT, &debug_level, NULL},
    {"-h", "-help", "display this message", NULL,
         SP_TYPE_BOOLEAN, &help_flag, SP_FALSE_STRING},
};
 
static char *filelabel[] = {
    "<input file>",
};
 
int spMain(int argc, char *argv[])
{
    long k;
    const char *input_file;
    spPlugin *i_plugin;
    spWaveInfo wave_info;
    long read_length;
    long fftl;
    double resid;
    spDVector ana_buf, window;
    spDVector cor, cor2, alpha, parcor, psd, psdfft, vtar, logvtar, vta;
    spDVector ovtar, ovta, oparcor, oalpha, ocor;
    spFFTRec fftrec;
            
    /* get command line options */
    spSetHelpMessage(&help_flag, "LPC test");
    options = spGetOptions(argc, argv, option, filelabel);
    spGetOptionsValue(argc, argv, options);
    spSetDebugLevel(debug_level);
    spCheckNumFile(options);
 
    if (frame_length <= 0) {
        spError(1, "Frame length must be larger than 0.");
    }
    if (order <= 0) {
        spError(1, "Analysis order must be larger than 0.");
    }
    frame_start = MAX(frame_start, 0);
 
    input_file = spGetFile(options);
    spDebug(-1, "main", "input_file = %s\n", input_file);
 
    spInitWaveInfo(&wave_info);
 
    if ((i_plugin = spOpenFilePluginArg(NULL, input_file, "r",
                                      SP_PLUGIN_DEVICE_FILE,
                                      &wave_info, NULL, 0, NULL, NULL)) != NULL) {
        ana_buf = xdvalloc(frame_length);
 
        spSeekPlugin(i_plugin, frame_start);
        
        if ((read_length = spReadPluginDouble(i_plugin, ana_buf->data, ana_buf->length)) > 0) {
            spDebug(-1, "main", "read_length = %ld, frame_length = %ld\n", read_length, frame_length);
            for (k = read_length; k < ana_buf->length; k++) {
                ana_buf->data[k] = 0.0;
            }
 
#if 0
            for (k = 0; k < read_length; k++) {
                printf("%f\n", ana_buf->data[k]);
            }
#endif
            fftl = POW2(spNextPow2(frame_length));
            
            cor = xdvalloc(order + 1);
            alpha = xdvalloc(order);
            psd = xdvalloc(fftl / 2 + 1);
            psdfft = xdvalloc(fftl);
            parcor = xdvalloc(order);
            vtar = xdvalloc(order);
            logvtar = xdvalloc(order);
            vta = xdvalloc(order);
            ovtar = xdvalloc(order);
            ovta = xdvalloc(order);
            oparcor = xdvalloc(order);
            oalpha = xdvalloc(order);
            ocor = xdvalloc(order + 1);
            
            window = xdvhamming(frame_length);
            dvoper(ana_buf, "*", window);
 
            spSigToCor(ana_buf->data, frame_length, cor->data, order + 1);
            
#if 1
            if ((fftrec = spInitFFT(spNextPow2(frame_length), SP_FFT_DOUBLE_PRECISION)) != NULL) {
                cor2 = xdvalloc(order + 1);
                
#if 0
                spSigToCorFFT(fftrec, ana_buf->data, frame_length, cor2->data, order + 1);
#else
                {
                    spDVector spec;
                    
                    spec = xdvalloc(spGetFFTLength(fftrec));
                    spCorToSpecFFT(fftrec, cor->data, cor->length, spec->data, spec->length, 0.0);
#if 0
                    for (k = 0; k < spec->length; k++) {
                        printf("%f\n", spec->data[k]);
                    }
#endif
                    spSpecToCorFFT(fftrec, spec->data, spec->length, cor2->data, cor2->length, 0.0);
 
                    xdvfree(spec);
                }
#endif
 
#if 0
                for (k = 0; k < cor->length; k++) {
                    printf("%f %f\n", cor->data[k], cor2->data[k]);
                }
#endif
                
                xdvfree(cor2);
                spFreeFFT(fftrec);
            }
#endif
 
            /* auto-correlation --> PARCOR coefficients */
            spCorToParcor(cor->data, alpha->data, parcor->data, order, &resid);
            spDebug(-1, "main", "resid = %f\n", resid);
 
            /* alpha --> power spectrum density */
            spAlphaToPSD(alpha->data, order, /*resid*/1.0, psd->data, psd->length);
            spAlphaToPSDFFT(alpha->data, order, /*resid*/1.0, psdfft->data, psdfft->length);
#if 0
            for (k = 0; k < psd->length; k++) {
                printf("%f %f\n", psd->data[k], psdfft->data[k]);
            }
#endif
 
#if 0
            /* PARCOR coefficients --> vocal tract area ratio function */
            spParcorToVTAreaRatio(parcor->data, vtar->data, order, SP_TRUE);
            spParcorToVTAreaRatio(parcor->data, logvtar->data, order, SP_FALSE);
 
            /* PARCOR coefficients --> vocal tract area function */
            spParcorToVTArea(parcor->data, ovta->data, order, 1.0);
            
            /* vocal tract area ratio function --> vocal tract area function */
            spVTRatioToArea(vtar->data, vta->data, order, 1.0, SP_TRUE);
 
            /* vocal tract area function --> vocal tract area ratio function */
            spVTAreaToRatio(vta->data, ovtar->data, order, 1.0, SP_TRUE);
            
            /* vocal tract area ratio function --> PARCOR coefficients */
            spVTAreaRatioToParcor(ovtar->data, oparcor->data, order, SP_TRUE);
#else
            /* PARCOR coefficients --> vocal tract area ratio function */
            spParcorToVTAreaRatio(parcor->data, vtar->data, order, SP_TRUE);
            
            /* PARCOR coefficients --> vocal tract area function */
            spParcorToVTArea(parcor->data, vta->data, order, 1.0);
            
            /* vocal tract area function --> vocal tract area ratio function */
            spVTAreaToRatio(vta->data, ovtar->data, order, 1.0, SP_TRUE);
            spVTAreaToRatio(vta->data, logvtar->data, order, 1.0, SP_FALSE);
            
            /* vocal tract area ratio function --> vocal tract area function */
            spVTRatioToArea(ovtar->data, ovta->data, order, 1.0, SP_TRUE);
            
            /* vocal tract area function --> PARCOR coefficients */
            spVTAreaToParcor(ovta->data, oparcor->data, order, 1.0);
 
#endif
            
            /* PARCOR coefficients --> alpha parameter & auto-correlation */
            spParcorToCor(oparcor->data, oalpha->data, order, cor->data[0], ocor->data, order + 1, NULL);
#if 1
            for (k = 0; k < order; k++) {
                printf("%f %f %f %f %f %f %f %f %f %f %f\n",
                       cor->data[k], ocor->data[k],
                       alpha->data[k], oalpha->data[k],
                       parcor->data[k], oparcor->data[k],
                       vtar->data[k], ovtar->data[k],
                       vta->data[k], ovta->data[k], 
                       logvtar->data[k]);
            }
            printf("%f %f\n", cor->data[k], ocor->data[k]);
#endif
            
 
            xdvfree(window);
            
            xdvfree(cor);
            xdvfree(alpha);
            xdvfree(psd);
            xdvfree(psdfft);
            xdvfree(parcor);
            xdvfree(vtar);
            xdvfree(logvtar);
            xdvfree(vta);
            xdvfree(ovtar);
            xdvfree(ovta);
            xdvfree(oparcor);
            xdvfree(ocor);
            xdvfree(oalpha);
        }
        
        xdvfree(ana_buf);
    }
    
    return 0;
}