Chris@366: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@366: /*
Chris@366:     Constant-Q library
Chris@366:     Copyright (c) 2013-2014 Queen Mary, University of London
Chris@366: 
Chris@366:     Permission is hereby granted, free of charge, to any person
Chris@366:     obtaining a copy of this software and associated documentation
Chris@366:     files (the "Software"), to deal in the Software without
Chris@366:     restriction, including without limitation the rights to use, copy,
Chris@366:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@366:     of the Software, and to permit persons to whom the Software is
Chris@366:     furnished to do so, subject to the following conditions:
Chris@366: 
Chris@366:     The above copyright notice and this permission notice shall be
Chris@366:     included in all copies or substantial portions of the Software.
Chris@366: 
Chris@366:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@366:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@366:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@366:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
Chris@366:     CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@366:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@366:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@366: 
Chris@366:     Except as contained in this notice, the names of the Centre for
Chris@366:     Digital Music; Queen Mary, University of London; and Chris Cannam
Chris@366:     shall not be used in advertising or otherwise to promote the sale,
Chris@366:     use or other dealings in this Software without prior written
Chris@366:     authorization.
Chris@366: */
Chris@366: 
Chris@366: #include "CQSpectrogram.h"
Chris@366: 
Chris@366: #include <iostream>
Chris@366: #include <stdexcept>
Chris@366: 
Chris@366: using std::cerr;
Chris@366: using std::endl;
Chris@366: 
Chris@366: //#define DEBUG_CQSPECTROGRAM 1
Chris@366: 
Chris@366: CQSpectrogram::CQSpectrogram(CQParameters params,
Chris@366:                              Interpolation interpolation) :
Chris@366:     m_cq(params),
Chris@366:     m_interpolation(interpolation)
Chris@366: {
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::~CQSpectrogram()
Chris@366: {
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::process(const RealSequence &td)
Chris@366: {
Chris@366:     return postProcess(m_cq.process(td), false);
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::getRemainingOutput()
Chris@366: {
Chris@366:     return postProcess(m_cq.getRemainingOutput(), true);
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::postProcess(const ComplexBlock &cq, bool insist)
Chris@366: {
Chris@366:     int width = cq.size();
Chris@366: 
Chris@366:     // convert to magnitudes
Chris@366:     RealBlock spec;
Chris@366:     for (int i = 0; i < width; ++i) {
Chris@366:         int height = cq[i].size();
Chris@366:         RealColumn col(height, 0);
Chris@366:         for (int j = 0; j < height; ++j) {
Chris@366: #ifdef DEBUG_CQSPECTROGRAM
Chris@366:             if (isnan(cq[i][j].real())) {
Chris@366:                 cerr << "WARNING: NaN in real at (" << i << "," << j << ")" << endl;
Chris@366:             }
Chris@366:             if (isnan(cq[i][j].imag())) {
Chris@366:                 cerr << "WARNING: NaN in imag at (" << i << "," << j << ")" << endl;
Chris@366:             }
Chris@366: #endif
Chris@366:             col[j] = abs(cq[i][j]);
Chris@366:         }
Chris@366:         spec.push_back(col);
Chris@366:     }
Chris@366: 
Chris@366:     if (m_interpolation == InterpolateZeros) {
Chris@366:         for (int i = 0; i < width; ++i) {
Chris@366:             int sh = spec[i].size();
Chris@366:             int fh = getTotalBins();
Chris@366:             for (int j = sh; j < fh; ++j) {
Chris@366:                 spec[i].push_back(0);
Chris@366:             }
Chris@366:         }
Chris@366: 	return spec;
Chris@366:     }
Chris@366: 
Chris@366:     for (int i = 0; i < width; ++i) {
Chris@366: 	m_buffer.push_back(spec[i]);
Chris@366:     }
Chris@366:     
Chris@366:     if (m_interpolation == InterpolateHold) {
Chris@366: 	return fetchHold(insist);
Chris@366:     } else {
Chris@366: 	return fetchLinear(insist);
Chris@366:     }
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::fetchHold(bool)
Chris@366: {
Chris@366:     RealBlock out;
Chris@366:     
Chris@366:     int width = m_buffer.size();
Chris@366:     int height = getTotalBins();
Chris@366: 
Chris@366:     for (int i = 0; i < width; ++i) {
Chris@366: 	
Chris@366: 	RealColumn col = m_buffer[i];
Chris@366: 
Chris@366: 	int thisHeight = col.size();
Chris@366: 	int prevHeight = m_prevColumn.size();
Chris@366: 
Chris@366: 	for (int j = thisHeight; j < height; ++j) {
Chris@366: 	    if (j < prevHeight) {
Chris@366: 		col.push_back(m_prevColumn[j]);
Chris@366: 	    } else {
Chris@366: 		col.push_back(0.0);
Chris@366: 	    }
Chris@366: 	}
Chris@366: 
Chris@366: 	m_prevColumn = col;
Chris@366: 	out.push_back(col);
Chris@366:     }
Chris@366: 
Chris@366:     m_buffer.clear();
Chris@366: 
Chris@366:     return out;
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::fetchLinear(bool insist)
Chris@366: {
Chris@366:     RealBlock out;
Chris@366: 
Chris@366:     //!!! This is surprisingly messy. I must be missing something.
Chris@366: 
Chris@366:     // We can only return any data when we have at least one column
Chris@366:     // that has the full height in the buffer, that is not the first
Chris@366:     // column.
Chris@366:     //
Chris@366:     // If the first col has full height, and there is another one
Chris@366:     // later that also does, then we can interpolate between those, up
Chris@366:     // to but not including the second full height column.  Then we
Chris@366:     // drop and return the columns we interpolated, leaving the second
Chris@366:     // full-height col as the first col in the buffer. And repeat as
Chris@366:     // long as enough columns are available.
Chris@366:     //
Chris@366:     // If the first col does not have full height, then (so long as
Chris@366:     // we're following the logic above) we must simply have not yet
Chris@366:     // reached the first full-height column in the CQ output, and we
Chris@366:     // can interpolate nothing.
Chris@366:     
Chris@366:     int width = m_buffer.size();
Chris@366:     int height = getTotalBins();
Chris@366: 
Chris@366:     if (width == 0) return out;
Chris@366:     
Chris@366:     int firstFullHeight = -1;
Chris@366:     int secondFullHeight = -1;
Chris@366: 
Chris@366:     for (int i = 0; i < width; ++i) {
Chris@366: 	if ((int)m_buffer[i].size() == height) {
Chris@366: 	    if (firstFullHeight == -1) {
Chris@366: 		firstFullHeight = i;
Chris@366: 	    } else if (secondFullHeight == -1) {
Chris@366: 		secondFullHeight = i;
Chris@366: 		break;
Chris@366: 	    }
Chris@366: 	}
Chris@366:     }
Chris@366: 
Chris@366: //    cerr << "fetchLinear: firstFullHeight = " << firstFullHeight << ", secondFullHeight = " << secondFullHeight << endl;
Chris@366: 
Chris@366:     if (firstFullHeight < 0) {
Chris@366: 	if (insist) {
Chris@366:             return fetchHold(true);
Chris@366: 	} else {
Chris@366: 	    return out;
Chris@366: 	}
Chris@366:     } else if (firstFullHeight > 0) {
Chris@366: 	// can interpolate nothing, stash up to first full height & recurse
Chris@366: 	out = RealBlock(m_buffer.begin(), m_buffer.begin() + firstFullHeight);
Chris@366: 	m_buffer = RealBlock(m_buffer.begin() + firstFullHeight, m_buffer.end());
Chris@366: 	RealBlock more = fetchLinear(insist);
Chris@366: 	out.insert(out.end(), more.begin(), more.end());
Chris@366: 	return out;
Chris@366:     } else if (secondFullHeight < 0) {
Chris@366: 	// firstFullHeight == 0, but there is no second full height --
Chris@366: 	// wait for it unless insist flag is set
Chris@366: 	if (insist) {
Chris@366:             return fetchHold(true);
Chris@366: 	} else {
Chris@366: 	    return out;
Chris@366: 	}
Chris@366:     } else {
Chris@366: 	// firstFullHeight == 0 and secondFullHeight also valid. Can interpolate
Chris@366: 	out = linearInterpolated(m_buffer, 0, secondFullHeight);
Chris@366: 	m_buffer = RealBlock(m_buffer.begin() + secondFullHeight, m_buffer.end());
Chris@366: 	RealBlock more = fetchLinear(insist);
Chris@366: 	out.insert(out.end(), more.begin(), more.end());
Chris@366: 	return out;
Chris@366:     }
Chris@366: }
Chris@366: 
Chris@366: CQSpectrogram::RealBlock
Chris@366: CQSpectrogram::linearInterpolated(const RealBlock &g, int x0, int x1)
Chris@366: {
Chris@366:     // g must be a grid with full-height columns at x0 and x1
Chris@366: 
Chris@366:     if (x0 >= x1) {
Chris@366: 	throw std::logic_error("x0 >= x1");
Chris@366:     }
Chris@366:     if (x1 >= (int)g.size()) {
Chris@366: 	throw std::logic_error("x1 >= g.size()");
Chris@366:     }
Chris@366:     if (g[x0].size() != g[x1].size()) {
Chris@366: 	throw std::logic_error("x0 and x1 are not the same height");
Chris@366:     }
Chris@366: 
Chris@366:     int height = g[x0].size();
Chris@366:     int width = x1 - x0;
Chris@366: 
Chris@366:     RealBlock out(g.begin() + x0, g.begin() + x1);
Chris@366: 
Chris@366:     for (int y = 0; y < height; ++y) {
Chris@366: 
Chris@366: 	int spacing = width;
Chris@366: 	for (int i = 1; i < width; ++i) {
Chris@366: 	    int thisHeight = g[x0 + i].size();
Chris@366: 	    if (thisHeight > height) {
Chris@366: 		throw std::logic_error("First column not full-height");
Chris@366: 	    }
Chris@366: 	    if (thisHeight > y) {
Chris@366: 		spacing = i;
Chris@366: 		break;
Chris@366: 	    }
Chris@366: 	}
Chris@366: 
Chris@366: 	if (spacing < 2) continue;
Chris@366: 
Chris@366: 	for (int i = 0; i + spacing <= width; i += spacing) {
Chris@366: 	    for (int j = 1; j < spacing; ++j) {
Chris@366: 		double proportion = double(j)/double(spacing);
Chris@366: 		double interpolated = 
Chris@366: 		    g[x0 + i][y] * (1.0 - proportion) +
Chris@366: 		    g[x0 + i + spacing][y] * proportion;
Chris@366: 		out[i + j].push_back(interpolated);
Chris@366: 	    }
Chris@366: 	}
Chris@366:     }
Chris@366: 
Chris@366:     return out;
Chris@366: }
Chris@366: 
Chris@366: 
Chris@366: