view qm-dsp-decimate/decimate.cpp @ 17:f8efad075df0

Better output, to match resample

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

#include "qm-dsp/dsp/rateconversion/Decimator.h"
#include "qm-dsp/maths/MathUtilities.h"

#include <iostream>
#include <sndfile.h>
#include <time.h>
#include <string>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <getopt.h>
#include <unistd.h>
#include <sys/time.h>

using namespace std;

int main(int argc, char **argv)
{
    int factor = 0;
    bool version = false;
    bool help = false;
    int c = 0;

    while (1) {
        int optionIndex = 0;

        static struct option longOpts[] = {
            { "help",          0, 0, 'h' },
            { "version",       0, 0, 'V' },
            { "by",            1, 0, 'b' },
            { 0, 0, 0, 0 }
        };

        c = getopt_long(argc, argv,
                        "b:hV",
                        longOpts, &optionIndex);
        if (c == -1) break;

        switch (c) {
        case 'h': help = true; break;
        case 'V': version = true; break;
        case 'b': factor = atoi(optarg); break;
        default:  help = true; break;
        }
    }

    if (version) {
	cerr << "decimate v0.1" << endl; //!!!
	return 0;
    }

    if (help || factor == 0 || optind + 2 != argc) {
	cerr << endl;
	cerr << "usage: " << argv[0] << " --by <factor> <infile.wav> <outfile.wav>" << endl;
	cerr << endl;
	return 2;
    }

    if (!MathUtilities::isPowerOfTwo(factor)) {
        cerr << "ERROR: Factor must be a power of two" << endl;
        return 1;
    }
    if (factor < 2 || factor > Decimator::getHighestSupportedFactor()) {
        cerr << "ERROR: Only factors between 2 and "
             << Decimator::getHighestSupportedFactor() 
             << " inclusive are supported" <<endl;
        return 1;
    }
    
    timeval tv;
    (void)gettimeofday(&tv, 0);

    char *fileName = strdup(argv[optind++]);
    char *fileNameOut = strdup(argv[optind++]);

    SNDFILE *sndfile;
    SNDFILE *sndfileOut;
    SF_INFO sfinfo;
    SF_INFO sfinfoOut;
    memset(&sfinfo, 0, sizeof(SF_INFO));

    sndfile = sf_open(fileName, SFM_READ, &sfinfo);
    if (!sndfile) {
	cerr << "ERROR: Failed to open input file \"" << fileName << "\": "
	     << sf_strerror(sndfile) << endl;
	return 1;
    }

    int sourceRate = sfinfo.samplerate;
    int targetRate = sourceRate / factor;

    cerr << "Decimating by factor " << factor << ", from " << sourceRate
         << " to " << targetRate << " Hz" << endl;

    if (sourceRate % factor != 0) {
        cerr << "WARNING: Decimation factor " << factor << " from source rate "
             << sourceRate << " gives non-integral target rate (not supported "
             << "by output file format), recording as " << targetRate << endl;
    }

    sfinfoOut.channels = sfinfo.channels;
    sfinfoOut.format = sfinfo.format;
    sfinfoOut.frames = sfinfo.frames / factor;
    sfinfoOut.samplerate = targetRate;
    sfinfoOut.sections = sfinfo.sections;
    sfinfoOut.seekable = sfinfo.seekable;

    sndfileOut = sf_open(fileNameOut, SFM_WRITE, &sfinfoOut);
    if (!sndfileOut) {
	cerr << "ERROR: Failed to open output file \"" << fileNameOut << "\" for writing: "
	     << sf_strerror(sndfileOut) << endl;
	return 1;
    }

    int channels = sfinfo.channels;
    vector<Decimator *> decimators; // one per channel

    int ibs = 1024;
    int obs = ibs / factor;

    for (int c = 0; c < channels; ++c) {
	decimators.push_back(new Decimator(ibs, factor));
    }

    float *ibuf = new float[channels * ibs];
    float *obuf = new float[channels * obs];
    
    double **prebuf = new double*[channels];
    double **postbuf = new double*[channels];

    for (int c = 0; c < channels; ++c) {
	prebuf[c] = new double[ibs];
	postbuf[c] = new double[obs];
    }
    
    int n = 0;
    
    while (true) {

	int count = sf_readf_float(sndfile, ibuf, ibs);
	if (count <= 0) break;
	
	for (int c = 0; c < channels; ++c) {

	    for (int i = 0; i < count; ++i) {
		prebuf[c][i] = ibuf[i * channels + c];
	    }
            for (int i = count; i < ibs; ++i) {
                prebuf[c][i] = 0.0;
            }

	    decimators[c]->process(prebuf[c], postbuf[c]);
	    
            for (int i = 0; i < obs; ++i) {
                obuf[i * channels + c] = postbuf[c][i];
            }
	}

        int ocount = obs;
        if (count < ibs) {
            ocount = count / factor;
        }

        sf_writef_float(sndfileOut, obuf, ocount);

	n += ocount;
    }

    sf_close(sndfile);
    sf_close(sndfileOut);

    for (int c = 0; c < channels; ++c) {
	delete[] prebuf[c];
	delete[] postbuf[c];
	delete decimators[c];
    }

    delete[] prebuf;
    delete[] postbuf;
    delete[] ibuf;
    delete[] obuf;

    timeval etv;
    (void)gettimeofday(&etv, 0);
        
    etv.tv_sec -= tv.tv_sec;
    if (etv.tv_usec < tv.tv_usec) {
	etv.tv_usec += 1000000;
            etv.tv_sec -= 1;
    }
    etv.tv_usec -= tv.tv_usec;
        
    double sec = double(etv.tv_sec) + (double(etv.tv_usec) / 1000000.0);
    cerr << sfinfo.frames << " frames in, " << n << " frames out"
	 << ", nominal factor " << (1.0/factor)
	 << ", actual " << double(n)/double(sfinfo.frames)
	 << endl << "Elapsed time: " << sec << " sec, in frames/sec: "
         << int(sfinfo.frames/sec) << ", out frames/sec: " << int(n/sec) << endl;

    return 0;
}