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