view LowFreq.cpp @ 5:bf742ae09443

Frequency range, and resample to higher frequency to avoid transition band -- and some comments on what is still wrong with it

#include "LowFreq.h"

#include "dsp/rateconversion/Resampler.h"
#include "dsp/transforms/FFT.h"

#include <cmath>
#include <cstdio>

using std::cerr;
using std::endl;
using std::vector;

static float minFmin = 0;
static float maxFmin = 500;
static float defaultFmin = 0.1;

static float minFmax = 0.1;
static float maxFmax = 500;
static float defaultFmax = 10;

static int minN = 1;
static int maxN = 16384;
static int defaultN = 64;

static float minOverlap = 0.0;
static float maxOverlap = 87.5;
static float overlapStep = 12.5;
static float defaultOverlap = 50.0;

LowFreq::LowFreq(float inputSampleRate) :
    Plugin(inputSampleRate),
    m_fmin(defaultFmin),
    m_fmax(defaultFmax),
    m_n(defaultN),
    m_overlap(defaultOverlap),
    m_blockSize(0),
    m_resampler(0),
    m_fft(0),
    m_window(0)
{
}

LowFreq::~LowFreq()
{
    delete m_resampler;
    delete m_fft;
    delete m_window;
}

string
LowFreq::getIdentifier() const
{
    return "lowfreq";
}

string
LowFreq::getName() const
{
    return "Low-frequency Spectrogram";
}

string
LowFreq::getDescription() const
{
    //!!! Return something helpful here!
    return "";
}

string
LowFreq::getMaker() const
{
    return "Queen Mary, University of London";
}

int
LowFreq::getPluginVersion() const
{
    return 1;
}

string
LowFreq::getCopyright() const
{
    return "GPL";
}

LowFreq::InputDomain
LowFreq::getInputDomain() const
{
    return TimeDomain;
}

size_t
LowFreq::getPreferredBlockSize() const
{
    return 1024;
}

size_t 
LowFreq::getPreferredStepSize() const
{
    return 1024;
}

size_t
LowFreq::getMinChannelCount() const
{
    return 1;
}

size_t
LowFreq::getMaxChannelCount() const
{
    return 1;
}

LowFreq::ParameterList
LowFreq::getParameterDescriptors() const
{
    ParameterList list;

    ParameterDescriptor d;
    d.identifier = "fmin";
    d.name = "Minimum frequency";
    d.description = "Low bound for frequencies to include in the spectrogram. This will be rounded up to the next highest bin frequency based on the highest frequency and number of bins.";
    d.unit = "Hz";
    d.minValue = minFmin;
    d.maxValue = maxFmin;
    d.defaultValue = defaultFmin;
    d.isQuantized = false;
    list.push_back(d);

    d.identifier = "fmax";
    d.name = "Maximum frequency";
    d.description = "Highest frequency component to include in the spectrogram.";
    d.unit = "Hz";
    d.minValue = minFmax;
    d.maxValue = maxFmax;
    d.defaultValue = defaultFmax;
    d.isQuantized = false;
    list.push_back(d);

    d.identifier = "n";
    d.name = "Number of bins";
    d.description = "Number of spectrogram bins to return.";
    d.unit = "";
    d.minValue = minN;
    d.maxValue = maxN;
    d.defaultValue = defaultN;
    d.isQuantized = true;
    d.quantizeStep = 1;
    list.push_back(d);

    d.identifier = "overlap";
    d.name = "Frame overlap";
    d.description = "Amount of overlap between one Hann-windowed frame and the next. Should be at least 50% to avoid gaps due to windowing.";
    d.unit = "%";
    d.minValue = minOverlap;
    d.maxValue = maxOverlap;
    d.defaultValue = defaultOverlap;
    d.isQuantized = true;
    d.quantizeStep = overlapStep;
    list.push_back(d);

    return list;
}

float
LowFreq::getParameter(string identifier) const
{
    if (identifier == "fmin") {
	return m_fmin;
    } else if (identifier == "fmax") {
	return m_fmax;
    } else if (identifier == "n") {
	return m_n;
    } else if (identifier == "overlap") {
	return m_overlap;
    }
    return 0;
}

void
LowFreq::setParameter(string identifier, float value) 
{
    if (identifier == "fmin") {
	m_fmin = value;
    } else if (identifier == "fmax") {
	m_fmax = value;
    } else if (identifier == "n") {
	m_n = int(value + 0.01);
    } else if (identifier == "overlap") {
	m_overlap = value;
    }
}

LowFreq::ProgramList
LowFreq::getPrograms() const
{
    ProgramList list;
    return list;
}

string
LowFreq::getCurrentProgram() const
{
    return ""; // no programs
}

void
LowFreq::selectProgram(string name)
{
}

LowFreq::OutputList
LowFreq::getOutputDescriptors() const
{
    OutputList list;

    OutputDescriptor d;
    d.identifier = "spectrogram";
    d.name = "Spectrogram";
    d.description = "";
    d.unit = "";
    d.hasFixedBinCount = true;
    d.binCount = m_n/2 + 1;

    d.hasKnownExtents = false;
    d.isQuantized = false;
    d.sampleType = OutputDescriptor::FixedSampleRate;

    d.sampleRate = float(getTargetSampleRate()) / float(getTargetStepSize());
    
    char namebuf[50];
    for (int i = 0; i <= m_n/2; ++i) {
	sprintf(namebuf, "%.3gHz", double(getOutputBinFrequency(i)));
	d.binNames.push_back(namebuf);
    }

    cerr << "output descriptor effective sample rate = " << d.sampleRate << endl;

    d.hasDuration = false;
    list.push_back(d);

    return list;
}

bool
LowFreq::initialise(size_t channels, size_t stepSize, size_t blockSize)
{
    if (channels < getMinChannelCount() ||
	channels > getMaxChannelCount()) {
	cerr << "LowFreq::initialise: ERROR: channels " << channels
	     << " out of acceptable range " << getMinChannelCount()
	     << " -> " << getMaxChannelCount() << endl;
	return false;
    }

    if (stepSize != blockSize) {
	// We don't actually care what the block size is, but there
	// must be no overlap
	cerr << "LowFreq::initialise: ERROR: step size and block size must be "
	     << "equal (" << stepSize << " != " << blockSize << ")" << endl;
	return false;
    }

    if (m_n < minN || m_n > maxN) {
	cerr << "LowFreq::initialise: ERROR: bin count " << m_n
	     << " out of acceptable range " << minN
	     << " -> " << maxN << endl;
	return false;
    }

    if (m_fmin < minFmin || m_fmin > maxFmin) {
	cerr << "LowFreq::initialise: ERROR: lowest frequency " << m_fmin
	     << " out of acceptable range " << minFmin
	     << " -> " << maxFmin << endl;
	return false;
    }

    if (m_fmax < minFmax || m_fmax > maxFmax) {
	cerr << "LowFreq::initialise: ERROR: highest frequency " << m_fmax
	     << " out of acceptable range " << minFmax
	     << " -> " << maxFmax << endl;
	return false;
    }

    if (m_overlap < minOverlap || m_overlap > maxOverlap) {
	cerr << "LowFreq::initialise: ERROR: overlap " << m_overlap
	     << " out of acceptable range " << minOverlap
	     << " -> " << maxOverlap << endl;
	return false;
    }

    if (fabs(m_inputSampleRate - round(m_inputSampleRate)) > 1e-6) {
	cerr << "LowFreq::initialise: WARNING: input sample rate " 
	     << m_inputSampleRate << " is non-integral, output frequencies "
	     << "will be skewed by rounding it to nearest integer" << endl;
    }

    m_blockSize = blockSize;

    reset();

    return true;
}

void
LowFreq::reset()
{
    delete m_resampler;
    delete m_fft;
    delete m_window;

    // We want to return a number of bins n between frequencies fmin
    // and fmax.

    // The bin frequency for bin n of an n-point fft is fs (the sample
    // rate), with the top half aliasing the bottom half. So in theory
    // we can downsample to a sample rate of fmax*2 and then run an
    // m-point FFT, where m is sufficient to give us n bins remaining
    // between fmin and fmax. For m then we calculate fmax /
    // ((fmax-fmin)/n).

    // However, our resampler introduces artifacts close to its cutoff
    // frequency within the filter transition band so we probably
    // shouldn't aim to resample to exactly fmax*2. Instead we
    // resample to fmax*4 and carry out an FFT of twice the expected
    // number of bins.

    m_resampler = new Resampler(int(round(m_inputSampleRate)),
				getTargetSampleRate());
    m_fft = new FFT(getFFTSize());
    m_window = new Window<double>(HanningWindow, getFFTSize());
    m_buffer = std::vector<double>();

    cerr << "LowFreq::reset: freq range " << m_fmin << " to " << m_fmax << endl;
    cerr << "LowFreq::reset: block size " << m_blockSize
	 << ", input sample rate " << m_inputSampleRate << endl;
    cerr << "LowFreq::reset: target rate " << getTargetSampleRate()
	 << ", target step " << getTargetStepSize()
	 << ", resampler latency " << m_resampler->getLatency()
	 << ", fft size " << getFFTSize() << endl;

    //!!! output has wrong number of bins

    //!!! not handling resampler latency

    //!!! how do we know the fft size is long enough to capture the low frequencies we want?
}

int
LowFreq::getTargetSampleRate() const
{
    // See notes in reset()
    return int(ceil(m_fmax)) * 4;
}

int
LowFreq::getTargetStepSize() const
{
    return int(round(m_n * (1.0 - (m_overlap / 100.0))));
}    

int
LowFreq::getFFTSize() const
{
    int fftSize = 2 * int(ceil(m_fmax / ((m_fmax - m_fmin) / m_n)));
    cerr << "LowFreq::getFFTSize: for range " << m_fmin << " -> " << m_fmax
	 << " and n = " << m_n << ", fft size is " << fftSize << endl;
    return fftSize;
}

int
LowFreq::getFirstOutputBin() const
{
    int first = int(floor(m_fmin / getFFTSize()));
    cerr << "LowFreq::getFirstOutputBin: for range " << m_fmin << " -> " << m_fmax
	 << " and n = " << m_n << ", bin is " << first << endl;
    return first;
}

float
LowFreq::getOutputBinFrequency(int i) const
{
    return (float(getTargetSampleRate()) / getFFTSize())
	* (i + getFirstOutputBin());
}

LowFreq::FeatureSet
LowFreq::process(const float *const *inputBuffers, Vamp::RealTime timestamp)
{
    double *data = new double[m_blockSize];
    for (int i = 0; i < m_blockSize; ++i) {
	data[i] = inputBuffers[0][i];
    }

    vector<double> resampled = m_resampler->process(data, m_blockSize);
    m_buffer.insert(m_buffer.end(), resampled.begin(), resampled.end());

    delete[] data;

    FeatureSet fs;

    while (int(m_buffer.size()) >= getFFTSize()) {
	Feature f = processColumn();
	fs[0].push_back(f);
	advance();
    }

    return fs;
}

LowFreq::FeatureSet
LowFreq::getRemainingFeatures()
{
    FeatureSet fs;

    while (!m_buffer.empty() && int(m_buffer.size()) < m_n) {
	m_buffer.push_back(0.0);
    }

    while (int(m_buffer.size()) >= getFFTSize()) {
	Feature f = processColumn();
	fs[0].push_back(f);
	advance();
    }

    return fs;
}

LowFreq::Feature
LowFreq::processColumn()
{
    Feature f;
    
    int sz = getFFTSize();

    double *realOut = new double[sz];
    double *imagOut = new double[sz];

    double *windowed = new double[sz];
    m_window->cut(m_buffer.data(), windowed);	

    m_fft->process(false, windowed, 0, realOut, imagOut);

    for (int i = 0; i <= m_n/2; ++i) {
	int ix = getFirstOutputBin() + i;
	f.values.push_back(realOut[ix] * realOut[ix] +
			   imagOut[ix] * imagOut[ix]);
    }

    return f;
}

void
LowFreq::advance()
{
    cerr << "step = " << getTargetStepSize() << endl;
    
    std::vector<double> advanced(m_buffer.data() + getTargetStepSize(),
				 m_buffer.data() + m_buffer.size());
    
    m_buffer = advanced;
}