Chris@27: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@27: 
Chris@27: /*
Chris@27:     Sonic Visualiser
Chris@27:     An audio file viewer and annotation editor.
Chris@27:     Centre for Digital Music, Queen Mary, University of London.
Chris@27:     This file copyright 2006 Chris Cannam.
Chris@27:     
Chris@27:     This program is free software; you can redistribute it and/or
Chris@27:     modify it under the terms of the GNU General Public License as
Chris@27:     published by the Free Software Foundation; either version 2 of the
Chris@27:     License, or (at your option) any later version.  See the file
Chris@27:     COPYING included with this distribution for more information.
Chris@27: */
Chris@27: 
Chris@27: #include "IntegerTimeStretcher.h"
Chris@27: 
Chris@27: #include <iostream>
Chris@27: #include <cassert>
Chris@27: 
Chris@27: //#define DEBUG_INTEGER_TIME_STRETCHER 1
Chris@27: 
Chris@27: IntegerTimeStretcher::IntegerTimeStretcher(size_t ratio,
Chris@27: 					   size_t maxProcessInputBlockSize,
Chris@27: 					   size_t inputIncrement,
Chris@27: 					   size_t windowSize,
Chris@27: 					   WindowType windowType) :
Chris@27:     m_ratio(ratio),
Chris@27:     m_n1(inputIncrement),
Chris@27:     m_n2(m_n1 * ratio),
Chris@27:     m_wlen(std::max(windowSize, m_n2 * 2)),
Chris@27:     m_inbuf(m_wlen),
Chris@27:     m_outbuf(maxProcessInputBlockSize * ratio)
Chris@27: {
Chris@27:     m_window = new Window<double>(windowType, m_wlen),
Chris@27: 
Chris@27:     m_time = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_wlen);
Chris@27:     m_freq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_wlen);
Chris@27:     m_dbuf = (double *)fftw_malloc(sizeof(double) * m_wlen);
Chris@27: 
Chris@27:     m_plan = fftw_plan_dft_1d(m_wlen, m_time, m_freq, FFTW_FORWARD, FFTW_ESTIMATE);
Chris@27:     m_iplan = fftw_plan_dft_c2r_1d(m_wlen, m_freq, m_dbuf, FFTW_ESTIMATE);
Chris@27: 
Chris@27:     m_mashbuf = new double[m_wlen];
Chris@27:     for (int i = 0; i < m_wlen; ++i) {
Chris@27: 	m_mashbuf[i] = 0.0;
Chris@27:     }
Chris@27: }
Chris@27: 
Chris@27: IntegerTimeStretcher::~IntegerTimeStretcher()
Chris@27: {
Chris@27:     std::cerr << "IntegerTimeStretcher::~IntegerTimeStretcher" << std::endl;
Chris@27: 
Chris@27:     fftw_destroy_plan(m_plan);
Chris@27:     fftw_destroy_plan(m_iplan);
Chris@27: 
Chris@27:     fftw_free(m_time);
Chris@27:     fftw_free(m_freq);
Chris@27:     fftw_free(m_dbuf);
Chris@27: 
Chris@27:     delete m_window;
Chris@27:     delete m_mashbuf;
Chris@27: }	
Chris@27: 
Chris@27: size_t
Chris@27: IntegerTimeStretcher::getProcessingLatency() const
Chris@27: {
Chris@27:     return getWindowSize() - getInputIncrement();
Chris@27: }
Chris@27: 
Chris@27: void
Chris@27: IntegerTimeStretcher::process(double *input, double *output, size_t samples)
Chris@27: {
Chris@27:     // We need to add samples from input to our internal buffer.  When
Chris@27:     // we have m_windowSize samples in the buffer, we can process it,
Chris@27:     // move the samples back by m_n1 and write the output onto our
Chris@27:     // internal output buffer.  If we have (samples * ratio) samples
Chris@27:     // in that, we can write m_n2 of them back to output and return
Chris@27:     // (otherwise we have to write zeroes).
Chris@27: 
Chris@27:     // When we process, we write m_wlen to our fixed output buffer
Chris@27:     // (m_mashbuf).  We then pull out the first m_n2 samples from that
Chris@27:     // buffer, push them into the output ring buffer, and shift
Chris@27:     // m_mashbuf left by that amount.
Chris@27: 
Chris@27:     // The processing latency is then m_wlen - m_n2.
Chris@27: 
Chris@27:     size_t consumed = 0;
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27:     std::cerr << "IntegerTimeStretcher::process(" << samples << ", consumed = " << consumed << "), writable " << m_inbuf.getWriteSpace() <<", readable "<< m_outbuf.getReadSpace() << std::endl;
Chris@27: #endif
Chris@27: 
Chris@27:     while (consumed < samples) {
Chris@27: 
Chris@27: 	size_t writable = m_inbuf.getWriteSpace();
Chris@27: 	writable = std::min(writable, samples - consumed);
Chris@27: 
Chris@27: 	if (writable == 0) {
Chris@27: 	    //!!! then what? I don't think this should happen, but
Chris@27: 	    std::cerr << "WARNING: IntegerTimeStretcher::process: writable == 0" << std::endl;
Chris@27: 	    break;
Chris@27: 	}
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27: 	std::cerr << "writing " << writable << " from index " << consumed << " to inbuf, consumed will be " << consumed + writable << std::endl;
Chris@27: #endif
Chris@27: 	m_inbuf.write(input + consumed, writable);
Chris@27: 	consumed += writable;
Chris@27: 
Chris@27: 	while (m_inbuf.getReadSpace() >= m_wlen &&
Chris@27: 	       m_outbuf.getWriteSpace() >= m_n2) {
Chris@27: 
Chris@27: 	    // We know we have at least m_wlen samples available
Chris@27: 	    // in m_inbuf.  We need to peek m_wlen of them for
Chris@27: 	    // processing, and then read m_n1 to advance the read
Chris@27: 	    // pointer.
Chris@27: 
Chris@27: 	    size_t got = m_inbuf.peek(m_dbuf, m_wlen);
Chris@27: 	    assert(got == m_wlen);
Chris@27: 		
Chris@27: 	    processBlock(m_dbuf, m_mashbuf);
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27: 	    std::cerr << "writing first " << m_n2 << " from mashbuf, skipping " << m_n1 << " on inbuf " << std::endl;
Chris@27: #endif
Chris@27: 	    m_inbuf.skip(m_n1);
Chris@27: 	    m_outbuf.write(m_mashbuf, m_n2);
Chris@27: 
Chris@27: 	    for (size_t i = 0; i < m_wlen - m_n2; ++i) {
Chris@27: 		m_mashbuf[i] = m_mashbuf[i + m_n2];
Chris@27: 	    }
Chris@27: 	    for (size_t i = m_wlen - m_n2; i < m_wlen; ++i) {
Chris@27: 		m_mashbuf[i] = 0.0f;
Chris@27: 	    }
Chris@27: 	}
Chris@27: 
Chris@27: //	std::cerr << "WARNING: IntegerTimeStretcher::process: writespace not enough for output increment (" << m_outbuf.getWriteSpace() << " < " << m_n2 << ")" << std::endl;
Chris@27: //	}
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27: 	std::cerr << "loop ended: inbuf read space " << m_inbuf.getReadSpace() << ", outbuf write space " << m_outbuf.getWriteSpace() << std::endl;
Chris@27: #endif
Chris@27:     }
Chris@27: 
Chris@27:     if (m_outbuf.getReadSpace() < samples * m_ratio) {
Chris@27: 	std::cerr << "WARNING: IntegerTimeStretcher::process: not enough data (yet?) (" << m_outbuf.getReadSpace() << " < " << (samples * m_ratio) << ")" << std::endl;
Chris@27: 	size_t fill = samples * m_ratio - m_outbuf.getReadSpace();
Chris@27: 	for (size_t i = 0; i < fill; ++i) {
Chris@27: 	    output[i] = 0.0;
Chris@27: 	}
Chris@27: 	m_outbuf.read(output + fill, m_outbuf.getReadSpace());
Chris@27:     } else {
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27: 	std::cerr << "enough data - writing " << samples * m_ratio << " from outbuf" << std::endl;
Chris@27: #endif
Chris@27: 	m_outbuf.read(output, samples * m_ratio);
Chris@27:     }
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27:     std::cerr << "IntegerTimeStretcher::process returning" << std::endl;
Chris@27: #endif
Chris@27: }
Chris@27: 
Chris@27: void
Chris@27: IntegerTimeStretcher::processBlock(double *buf, double *out)
Chris@27: {
Chris@27:     size_t i;
Chris@27: 
Chris@27:     // buf contains m_wlen samples; out contains enough space for
Chris@27:     // m_wlen * ratio samples (we mix into out, rather than replacing)
Chris@27: 
Chris@27: #ifdef DEBUG_INTEGER_TIME_STRETCHER
Chris@27:     std::cerr << "IntegerTimeStretcher::processBlock" << std::endl;
Chris@27: #endif
Chris@27: 
Chris@27:     m_window->cut(buf);
Chris@27: 
Chris@27:     for (i = 0; i < m_wlen/2; ++i) {
Chris@27: 	double temp = buf[i];
Chris@27: 	buf[i] = buf[i + m_wlen/2];
Chris@27: 	buf[i + m_wlen/2] = temp;
Chris@27:     }
Chris@27:     
Chris@27:     for (i = 0; i < m_wlen; ++i) {
Chris@27: 	m_time[i][0] = buf[i];
Chris@27: 	m_time[i][1] = 0.0;
Chris@27:     }
Chris@27: 
Chris@27:     fftw_execute(m_plan); // m_time -> m_freq
Chris@27: 
Chris@27:     for (i = 0; i < m_wlen; ++i) {
Chris@27: 	
Chris@27: 	double mag = sqrt(m_freq[i][0] * m_freq[i][0] +
Chris@27: 			  m_freq[i][1] * m_freq[i][1]);
Chris@27: 		
Chris@27: 	double phase = atan2(m_freq[i][1], m_freq[i][0]);
Chris@27: 	
Chris@27: 	phase = phase * m_ratio;
Chris@27: 	
Chris@27: 	double real = mag * cos(phase);
Chris@27: 	double imag = mag * sin(phase);
Chris@27: 	m_freq[i][0] = real;
Chris@27: 	m_freq[i][1] = imag;
Chris@27:     }
Chris@27:     
Chris@27:     fftw_execute(m_iplan); // m_freq -> in, inverse fft
Chris@27:     
Chris@27:     for (i = 0; i < m_wlen/2; ++i) {
Chris@27: 	double temp = buf[i] / m_wlen;
Chris@27: 	buf[i] = buf[i + m_wlen/2] / m_wlen;
Chris@27: 	buf[i + m_wlen/2] = temp;
Chris@27:     }
Chris@27:     
Chris@27:     m_window->cut(buf);
Chris@27:     
Chris@27:     int div = m_wlen / m_n2;
Chris@27:     if (div > 1) div /= 2;
Chris@27:     for (i = 0; i < m_wlen; ++i) {
Chris@27: 	buf[i] /= div;
Chris@27:     }
Chris@27: 
Chris@27:     for (i = 0; i < m_wlen; ++i) {
Chris@27: 	out[i] += buf[i];
Chris@27:     }
Chris@27: }