Mercurial > hg > constant-q-cpp
changeset 123:e5c36c054c48
Merge from branch standalone-reorg
author | Chris Cannam <c.cannam@qmul.ac.uk> |
---|---|
date | Thu, 15 May 2014 14:42:11 +0100 |
parents | 93be4aa255e5 (current diff) edbec47f4a3d (diff) |
children | fcf7e33aa56b |
files | cpp-qm-dsp/.keep cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQInverse.cpp cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQKernel.cpp cpp-qm-dsp/CQKernel.h cpp-qm-dsp/CQSpectrogram.cpp cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.cpp cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/processfile.cpp cpp-qm-dsp/test.cpp cpp-standalone/.keep paper/.keep paper/smc2010.pdf yeti/.keep yeti/Makefile yeti/build.xml yeti/cqt.yeti yeti/cqtkernel.yeti yeti/experiment.yeti yeti/icqt.yeti yeti/nbproject/ide-file-targets.xml yeti/nbproject/project.xml yeti/plotfile.yeti yeti/test.yeti yeti/test_cqtkernel.yeti yeti/test_frequency.yeti |
diffstat | 67 files changed, 5383 insertions(+), 3334 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgsub Thu May 15 11:59:11 2014 +0100 +++ b/.hgsub Thu May 15 14:42:11 2014 +0100 @@ -1,1 +1,2 @@ -matlab = https://code.soundsoftware.ac.uk/hg/constant-q-toolbox +src/ext/kissfft = http://hg.code.sf.net/p/kissfft/code +misc/matlab = https://code.soundsoftware.ac.uk/hg/constant-q-toolbox
--- a/.hgsubstate Thu May 15 11:59:11 2014 +0100 +++ b/.hgsubstate Thu May 15 14:42:11 2014 +0100 @@ -1,1 +1,3 @@ 2b03ca77abcc91140dbb37827f06034769c53547 matlab +2b03ca77abcc91140dbb37827f06034769c53547 misc/matlab +fbe1bb0bc7b94ec252842b8b7e3f3347ec75d92f src/ext/kissfft
--- a/COPYING Thu May 15 11:59:11 2014 +0100 +++ b/COPYING Thu May 15 14:42:11 2014 +0100 @@ -1,12 +1,6 @@ -The code in this library is provided under an MIT-style licence as -follows. However, note that other code which it depends on may be -licensed differently (for example, qm-dsp is under the GPL). This -licence applies only to the source files actually contained in this -repository, and not to any of its subrepositories. See individual -files for more details. +This library is provided to you under the following licence: - Constant-Q library Copyright (c) 2013-2014 Queen Mary, University of London Permission is hereby granted, free of charge, to any person @@ -34,4 +28,38 @@ use or other dealings in this Software without prior written authorization. +The code in src/ext/kissfft is under the following licence: + Copyright (c) 2003-2010 Mark Borgerding + + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the author nor the names of any contributors may be used + to endorse or promote products derived from this software + without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND + CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, + INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS + BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED + TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR + TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF + THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + SUCH DAMAGE.
--- a/Makefile.inc Thu May 15 11:59:11 2014 +0100 +++ b/Makefile.inc Thu May 15 14:42:11 2014 +0100 @@ -1,8 +1,10 @@ -LIB_DIR := cpp-qm-dsp +LIB_DIR := src +INC_DIR := cq +TEST_DIR := test VAMP_DIR := vamp +KFFT_DIR := src/ext/kissfft -QMDSP_DIR ?= ../qm-dsp VAMPSDK_DIR ?= ../vamp-plugin-sdk PLUGIN_EXT ?= .so @@ -10,8 +12,9 @@ CXX ?= g++ CC ?= gcc -CFLAGS := $(CFLAGS) -CXXFLAGS := -I. -I$(VAMPSDK_DIR) -I$(QMDSP_DIR) $(CXXFLAGS) +GENERAL_FLAGS := -I. -I$(VAMPSDK_DIR) -I$(INC_DIR) -I$(LIB_DIR) -I$(KFFT_DIR) -I$(KFFT_DIR)/tools -Dkiss_fft_scalar=double +CFLAGS := $(GENERAL_FLAGS) $(CFLAGS) +CXXFLAGS := $(GENERAL_FLAGS) $(CXXFLAGS) LDFLAGS := $(LDFLAGS) PLUGIN_LDFLAGS := $(LDFLAGS) $(PLUGIN_LDFLAGS) @@ -19,26 +22,27 @@ PF_LDFLAGS := $(LDFLAGS) -lsndfile PLUGIN := cqvamp$(PLUGIN_EXT) -TEST := $(LIB_DIR)/test -PF := $(LIB_DIR)/processfile +TEST := $(TEST_DIR)/test +PF := $(TEST_DIR)/processfile -LIB_HEADERS := $(LIB_DIR)/CQBase.h $(LIB_DIR)/CQKernel.h $(LIB_DIR)/ConstantQ.h $(LIB_DIR)/CQSpectrogram.h $(LIB_DIR)/CQInverse.h -LIB_SOURCES := $(LIB_DIR)/CQKernel.cpp $(LIB_DIR)/ConstantQ.cpp $(LIB_DIR)/CQSpectrogram.cpp $(LIB_DIR)/CQInverse.cpp +LIB_HEADERS := $(INC_DIR)/CQBase.h $(LIB_DIR)/CQKernel.h $(INC_DIR)/ConstantQ.h $(INC_DIR)/CQSpectrogram.h $(INC_DIR)/CQInverse.h $(LIB_DIR)/dsp/FFT.h $(LIB_DIR)/dsp/KaiserWindow.h $(LIB_DIR)/dsp/MathUtilities.h $(LIB_DIR)/dsp/nan-inf.h $(LIB_DIR)/dsp/Resampler.h $(LIB_DIR)/dsp/SincWindow.h $(LIB_DIR)/dsp/Window.h $(KFFT_DIR)/kiss_fft.h $(KFFT_DIR)/tools/kiss_fftr.h +LIB_SOURCES := $(LIB_DIR)/CQKernel.cpp $(LIB_DIR)/ConstantQ.cpp $(LIB_DIR)/CQSpectrogram.cpp $(LIB_DIR)/CQInverse.cpp $(LIB_DIR)/dsp/FFT.cpp $(LIB_DIR)/dsp/KaiserWindow.cpp $(LIB_DIR)/dsp/MathUtilities.cpp $(LIB_DIR)/dsp/Resampler.cpp $(LIB_DIR)/dsp/SincWindow.cpp $(KFFT_DIR)/kiss_fft.c $(KFFT_DIR)/tools/kiss_fftr.c -VAMP_HEADERS := $(VAMP_DIR)/CQVamp.h $(VAMP_DIR)/CQChromaVamp.h -VAMP_SOURCES := $(VAMP_DIR)/CQVamp.cpp $(VAMP_DIR)/CQChromaVamp.cpp $(VAMP_DIR)/libmain.cpp +VAMP_HEADERS := $(VAMP_DIR)/CQVamp.h $(VAMP_DIR)/CQChromaVamp.h $(VAMP_DIR)/Pitch.h +VAMP_SOURCES := $(VAMP_DIR)/CQVamp.cpp $(VAMP_DIR)/CQChromaVamp.cpp $(VAMP_DIR)/libmain.cpp $(VAMP_DIR)/Pitch.cpp HEADERS := $(LIB_HEADERS) $(VAMP_HEADERS) SOURCES := $(LIB_SOURCES) $(VAMP_SOURCES) OBJECTS := $(SOURCES:.cpp=.o) +OBJECTS := $(OBJECTS:.c=.o) -TEST_SOURCES := $(LIB_DIR)/test.cpp +TEST_SOURCES := $(TEST_DIR)/test.cpp TEST_OBJECTS := $(TEST_SOURCES:.cpp=.o) $(OBJECTS) -PF_SOURCES := $(LIB_DIR)/processfile.cpp +PF_SOURCES := $(TEST_DIR)/processfile.cpp PF_OBJECTS := $(PF_SOURCES:.cpp=.o) $(OBJECTS) -LIBS := $(QMDSP_DIR)/libqm-dsp.a $(VAMPSDK_DIR)/libvamp-sdk.a -lpthread +LIBS := $(VAMPSDK_DIR)/libvamp-sdk.a -lpthread all: $(PLUGIN) $(TEST) $(PF) @@ -62,29 +66,15 @@ # DO NOT DELETE -cpp-qm-dsp/CQKernel.o: cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h -cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQInverse.h cpp-qm-dsp/CQBase.h -cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQKernel.h -vamp/CQVamp.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h -vamp/CQVamp.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -vamp/CQVamp.o: cpp-qm-dsp/CQKernel.h -vamp/CQChromaVamp.o: vamp/CQChromaVamp.h cpp-qm-dsp/CQSpectrogram.h -vamp/CQChromaVamp.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -vamp/CQChromaVamp.o: cpp-qm-dsp/CQKernel.h -vamp/libmain.o: vamp/CQVamp.h cpp-qm-dsp/CQSpectrogram.h -vamp/libmain.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -vamp/libmain.o: cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/test.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h -cpp-qm-dsp/test.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/processfile.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -cpp-qm-dsp/processfile.o: cpp-qm-dsp/CQKernel.h cpp-qm-dsp/CQInverse.h -cpp-qm-dsp/ConstantQ.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/ConstantQ.h cpp-qm-dsp/CQBase.h -cpp-qm-dsp/CQSpectrogram.o: cpp-qm-dsp/CQKernel.h -cpp-qm-dsp/CQInverse.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h -vamp/CQVamp.o: cpp-qm-dsp/CQSpectrogram.h cpp-qm-dsp/ConstantQ.h -vamp/CQVamp.o: cpp-qm-dsp/CQBase.h cpp-qm-dsp/CQKernel.h +src/CQKernel.o: src/CQKernel.h src/dsp/MathUtilities.h src/dsp/nan-inf.h +src/CQKernel.o: src/dsp/FFT.h src/dsp/Window.h +src/ConstantQ.o: src/CQKernel.h src/dsp/Resampler.h src/dsp/MathUtilities.h +src/ConstantQ.o: src/dsp/nan-inf.h src/dsp/FFT.h +src/CQInverse.o: src/dsp/Resampler.h src/dsp/MathUtilities.h +src/CQInverse.o: src/dsp/nan-inf.h src/dsp/FFT.h +vamp/CQVamp.o: vamp/CQVamp.h +vamp/CQChromaVamp.o: vamp/CQChromaVamp.h +vamp/libmain.o: vamp/CQVamp.h vamp/CQChromaVamp.h +cq/ConstantQ.o: cq/CQBase.h src/CQKernel.h +cq/CQSpectrogram.o: cq/ConstantQ.h cq/CQBase.h src/CQKernel.h +cq/CQInverse.o: cq/CQBase.h src/CQKernel.h
--- a/cpp-qm-dsp/CQBase.h Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#ifndef CQBASE_H -#define CQBASE_H - -#include <vector> -#include <complex> - -class CQBase // interface class -{ -public: - typedef std::complex<double> Complex; - typedef std::vector<double> RealSequence; - typedef std::vector<double> RealColumn; - typedef std::vector<Complex> ComplexSequence; - typedef std::vector<Complex> ComplexColumn; - typedef std::vector<RealColumn> RealBlock; - typedef std::vector<ComplexColumn> ComplexBlock; - - virtual double getSampleRate() const = 0; - virtual int getBinsPerOctave() const = 0; - virtual int getOctaves() const = 0; - virtual int getTotalBins() const = 0; - virtual int getColumnHop() const = 0; - virtual int getLatency() const = 0; - virtual double getMaxFrequency() const = 0; - virtual double getMinFrequency() const = 0; // actual min, not that passed to ctor - virtual double getBinFrequency(int bin) const = 0; -}; - -#endif
--- a/cpp-qm-dsp/CQInverse.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,411 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#include "CQInverse.h" - -#include "dsp/rateconversion/Resampler.h" -#include "maths/MathUtilities.h" -#include "dsp/transforms/FFT.h" - -#include <algorithm> -#include <iostream> -#include <stdexcept> - -using std::vector; -using std::cerr; -using std::endl; - -//#define DEBUG_CQ 1 - -CQInverse::CQInverse(double sampleRate, - double minFreq, - double maxFreq, - int binsPerOctave) : - m_sampleRate(sampleRate), - m_maxFrequency(maxFreq), - m_minFrequency(minFreq), - m_binsPerOctave(binsPerOctave), - m_fft(0) -{ - if (minFreq <= 0.0 || maxFreq <= 0.0) { - throw std::invalid_argument("Frequency extents must be positive"); - } - - initialise(); -} - -CQInverse::~CQInverse() -{ - delete m_fft; - for (int i = 0; i < (int)m_upsamplers.size(); ++i) { - delete m_upsamplers[i]; - } - delete m_kernel; -} - -double -CQInverse::getMinFrequency() const -{ - return m_p.minFrequency / pow(2.0, m_octaves - 1); -} - -double -CQInverse::getBinFrequency(int bin) const -{ - return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave())); -} - -void -CQInverse::initialise() -{ - m_octaves = int(ceil(log2(m_maxFrequency / m_minFrequency))); - m_kernel = new CQKernel(m_sampleRate, m_maxFrequency, m_binsPerOctave); - m_p = m_kernel->getProperties(); - - // Use exact powers of two for resampling rates. They don't have - // to be related to our actual samplerate: the resampler only - // cares about the ratio, but it only accepts integer source and - // target rates, and if we start from the actual samplerate we - // risk getting non-integer rates for lower octaves - - int sourceRate = pow(2, m_octaves); - vector<int> latencies; - - // top octave, no resampling - latencies.push_back(0); - m_upsamplers.push_back(0); - - for (int i = 1; i < m_octaves; ++i) { - - int factor = pow(2, i); - - Resampler *r = new Resampler - (sourceRate / factor, sourceRate, 50, 0.05); - -#ifdef DEBUG_CQ - cerr << "inverse: octave " << i << ": resample from " << sourceRate/factor << " to " << sourceRate << endl; -#endif - - // See ConstantQ.cpp for discussion on latency -- output - // latency here is at target rate which, this way around, is - // what we want - - latencies.push_back(r->getLatency()); - m_upsamplers.push_back(r); - } - - // additionally we will have fftHop latency at individual octave - // rate (before upsampling) for the overlap-add in each octave - for (int i = 0; i < m_octaves; ++i) { - latencies[i] += m_p.fftHop * pow(2, i); - } - - // Now reverse the drop adjustment made in ConstantQ to align the - // atom centres across different octaves (but this time at output - // sample rate) - - int emptyHops = m_p.firstCentre / m_p.atomSpacing; - - vector<int> pushes; - for (int i = 0; i < m_octaves; ++i) { - int factor = pow(2, i); - int pushHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops; - int push = ((pushHops * m_p.fftHop) * factor) / m_p.atomsPerFrame; - pushes.push_back(push); - } - - int maxLatLessPush = 0; - for (int i = 0; i < m_octaves; ++i) { - int latLessPush = latencies[i] - pushes[i]; - if (latLessPush > maxLatLessPush) maxLatLessPush = latLessPush; - } - - int totalLatency = maxLatLessPush + 10; - if (totalLatency < 0) totalLatency = 0; - - m_outputLatency = totalLatency + m_p.firstCentre * pow(2, m_octaves-1); - -#ifdef DEBUG_CQ - cerr << "totalLatency = " << totalLatency << ", m_outputLatency = " << m_outputLatency << endl; -#endif - - for (int i = 0; i < m_octaves; ++i) { - - // Calculate the difference between the total latency applied - // across all octaves, and the existing latency due to the - // upsampler for this octave. - - int latencyPadding = totalLatency - latencies[i] + pushes[i]; - -#ifdef DEBUG_CQ - cerr << "octave " << i << ": push " << pushes[i] << ", resampler latency inc overlap space " << latencies[i] << ", latencyPadding = " << latencyPadding << " (/factor = " << latencyPadding / pow(2, i) << ")" << endl; -#endif - - m_buffers.push_back(RealSequence(latencyPadding, 0.0)); - } - - for (int i = 0; i < m_octaves; ++i) { - // Fixed-size buffer for IFFT overlap-add - m_olaBufs.push_back(RealSequence(m_p.fftSize, 0.0)); - } - - m_fft = new FFTReal(m_p.fftSize); -} - -CQInverse::RealSequence -CQInverse::process(const ComplexBlock &block) -{ - // The input data is of the form produced by ConstantQ::process -- - // an unknown number N of columns of varying height. We assert - // that N is a multiple of atomsPerFrame * 2^(octaves-1), as must - // be the case for data that came directly from our ConstantQ - // implementation. - - int widthProvided = block.size(); - - if (widthProvided == 0) { - return drawFromBuffers(); - } - - int blockWidth = m_p.atomsPerFrame * int(pow(2, m_octaves - 1)); - - if (widthProvided % blockWidth != 0) { - cerr << "ERROR: CQInverse::process: Input block size (" - << widthProvided - << ") must be a multiple of processing block width " - << "(atoms-per-frame * 2^(octaves-1) = " - << m_p.atomsPerFrame << " * 2^(" << m_octaves << "-1) = " - << blockWidth << ")" << endl; - throw std::invalid_argument - ("Input block size must be a multiple of processing block width"); - } - - // Procedure: - // - // 1. Slice the list of columns into a set of lists of columns, - // one per octave, each of width N / (2^octave-1) and height - // binsPerOctave, containing the values present in that octave - // - // 2. Group each octave list by atomsPerFrame columns at a time, - // and stack these so as to achieve a list, for each octave, of - // taller columns of height binsPerOctave * atomsPerFrame - // - // 3. For each taller column, take the product with the inverse CQ - // kernel (which is the conjugate of the forward kernel) and - // perform an inverse FFT - // - // 4. Overlap-add each octave's resynthesised blocks (unwindowed) - // - // 5. Resample each octave's overlap-add stream to the original - // rate - // - // 6. Sum the resampled streams and return - - for (int i = 0; i < m_octaves; ++i) { - - // Step 1 - - ComplexBlock oct; - - for (int j = 0; j < widthProvided; ++j) { - int h = block[j].size(); - if (h < m_binsPerOctave * (i+1)) { - continue; - } - ComplexColumn col(block[j].begin() + m_binsPerOctave * i, - block[j].begin() + m_binsPerOctave * (i+1)); - oct.push_back(col); - } - - // Steps 2, 3, 4, 5 - processOctave(i, oct); - } - - // Step 6 - return drawFromBuffers(); -} - -CQInverse::RealSequence -CQInverse::drawFromBuffers() -{ - // 6. Sum the resampled streams and return - - int available = 0; - - for (int i = 0; i < m_octaves; ++i) { - if (i == 0 || int(m_buffers[i].size()) < available) { - available = m_buffers[i].size(); - } - } - - RealSequence result(available, 0); - - if (available == 0) { - return result; - } - - for (int i = 0; i < m_octaves; ++i) { - for (int j = 0; j < available; ++j) { - result[j] += m_buffers[i][j]; - } - m_buffers[i] = RealSequence(m_buffers[i].begin() + available, - m_buffers[i].end()); - } - - return result; -} - -CQInverse::RealSequence -CQInverse::getRemainingOutput() -{ - for (int j = 0; j < m_octaves; ++j) { - int factor = pow(2, j); - int latency = (j > 0 ? m_upsamplers[j]->getLatency() : 0) / factor; - for (int i = 0; i < (latency + m_p.fftSize) / m_p.fftHop; ++i) { - overlapAddAndResample(j, RealSequence(m_olaBufs[j].size(), 0)); - } - } - - return drawFromBuffers(); -} - -void -CQInverse::processOctave(int octave, const ComplexBlock &columns) -{ - // 2. Group each octave list by atomsPerFrame columns at a time, - // and stack these so as to achieve a list, for each octave, of - // taller columns of height binsPerOctave * atomsPerFrame - - int ncols = columns.size(); - - if (ncols % m_p.atomsPerFrame != 0) { - cerr << "ERROR: CQInverse::process: Number of columns (" - << ncols - << ") in octave " << octave - << " must be a multiple of atoms-per-frame (" - << m_p.atomsPerFrame << ")" << endl; - throw std::invalid_argument - ("Columns in octave must be a multiple of atoms per frame"); - } - - for (int i = 0; i < ncols; i += m_p.atomsPerFrame) { - - ComplexColumn tallcol; - for (int b = 0; b < m_binsPerOctave; ++b) { - for (int a = 0; a < m_p.atomsPerFrame; ++a) { - tallcol.push_back(columns[i + a][m_binsPerOctave - b - 1]); - } - } - - processOctaveColumn(octave, tallcol); - } -} - -void -CQInverse::processOctaveColumn(int octave, const ComplexColumn &column) -{ - // 3. For each taller column, take the product with the inverse CQ - // kernel (which is the conjugate of the forward kernel) and - // perform an inverse FFT - - if ((int)column.size() != m_p.atomsPerFrame * m_binsPerOctave) { - cerr << "ERROR: CQInverse::processOctaveColumn: Height of column (" - << column.size() << ") in octave " << octave - << " must be atoms-per-frame * bins-per-octave (" - << m_p.atomsPerFrame << " * " << m_binsPerOctave << " = " - << m_p.atomsPerFrame * m_binsPerOctave << ")" << endl; - throw std::invalid_argument - ("Column height must match atoms-per-frame * bins-per-octave"); - } - - ComplexSequence transformed = m_kernel->processInverse(column); - - int halfLen = m_p.fftSize/2 + 1; - - RealSequence ri(halfLen, 0); - RealSequence ii(halfLen, 0); - - for (int i = 0; i < halfLen; ++i) { - ri[i] = transformed[i].real(); - ii[i] = transformed[i].imag(); - } - - RealSequence timeDomain(m_p.fftSize, 0); - - m_fft->inverse(ri.data(), ii.data(), timeDomain.data()); - - overlapAddAndResample(octave, timeDomain); -} - -void -CQInverse::overlapAddAndResample(int octave, const RealSequence &seq) -{ - // 4. Overlap-add each octave's resynthesised blocks (unwindowed) - // - // and - // - // 5. Resample each octave's overlap-add stream to the original - // rate - - if (seq.size() != m_olaBufs[octave].size()) { - cerr << "ERROR: CQInverse::overlapAdd: input sequence length (" - << seq.size() << ") is expected to match OLA buffer size (" - << m_olaBufs[octave].size() << ")" << endl; - throw std::invalid_argument - ("Input sequence length should match OLA buffer size"); - } - - RealSequence toResample(m_olaBufs[octave].begin(), - m_olaBufs[octave].begin() + m_p.fftHop); - - RealSequence resampled = - octave > 0 ? - m_upsamplers[octave]->process(toResample.data(), toResample.size()) : - toResample; - - m_buffers[octave].insert(m_buffers[octave].end(), - resampled.begin(), - resampled.end()); - - m_olaBufs[octave] = RealSequence(m_olaBufs[octave].begin() + m_p.fftHop, - m_olaBufs[octave].end()); - - RealSequence pad(m_p.fftHop, 0); - - m_olaBufs[octave].insert(m_olaBufs[octave].end(), - pad.begin(), - pad.end()); - - for (int i = 0; i < m_p.fftSize; ++i) { - m_olaBufs[octave][i] += seq[i]; - } -} -
--- a/cpp-qm-dsp/CQInverse.h Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,91 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#ifndef CQINVERSE_H -#define CQINVERSE_H - -#include "CQBase.h" -#include "CQKernel.h" - -class Resampler; -class FFTReal; - -class CQInverse : public CQBase -{ -public: - CQInverse(double sampleRate, - double minFreq, double maxFreq, - int binsPerOctave); - virtual ~CQInverse(); - - virtual double getSampleRate() const { return m_sampleRate; } - virtual int getBinsPerOctave() const { return m_binsPerOctave; } - virtual int getOctaves() const { return m_octaves; } - virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; } - virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; } - virtual int getLatency() const { return m_outputLatency; } - virtual double getMaxFrequency() const { return m_p.maxFrequency; } - virtual double getMinFrequency() const; // actual min, not that passed to ctor - virtual double getBinFrequency(int bin) const; - - // Input is the format produced by ConstantQ class, - // i.e. uninterpolated complex, not the real-valued stuff produced - // by CQSpectrogram - - RealSequence process(const ComplexBlock &); - RealSequence getRemainingOutput(); - -private: - double m_sampleRate; - double m_maxFrequency; - double m_minFrequency; - int m_binsPerOctave; - int m_octaves; - - CQKernel *m_kernel; - CQKernel::Properties m_p; - - std::vector<Resampler *> m_upsamplers; - std::vector<RealSequence> m_buffers; - std::vector<RealSequence> m_olaBufs; // fixed-length, for overlap-add - - int m_outputLatency; - - FFTReal *m_fft; - - void initialise(); - void processOctave(int octave, const ComplexBlock &block); - void processOctaveColumn(int octave, const ComplexColumn &column); - void overlapAddAndResample(int octave, const RealSequence &); - RealSequence drawFromBuffers(); -}; - -#endif
--- a/cpp-qm-dsp/CQKernel.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,337 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#include "CQKernel.h" - -#include "maths/MathUtilities.h" -#include "dsp/transforms/FFT.h" -#include "base/Window.h" - -#include <cmath> -#include <cassert> -#include <vector> -#include <iostream> -#include <algorithm> - -using std::vector; -using std::complex; -using std::cerr; -using std::endl; - -typedef std::complex<double> C; - -CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave) : - m_fft(0) -{ - m_p.sampleRate = sampleRate; - m_p.maxFrequency = maxFreq; - m_p.binsPerOctave = binsPerOctave; - generateKernel(); -} - -CQKernel::~CQKernel() -{ - delete m_fft; -} - -void -CQKernel::generateKernel() -{ - double q = 1; - double atomHopFactor = 0.25; - double thresh = 0.0005; - - double bpo = m_p.binsPerOctave; - - m_p.minFrequency = (m_p.maxFrequency / 2) * pow(2, 1.0/bpo); - m_p.Q = q / (pow(2, 1.0/bpo) - 1.0); - - double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency); - double minNK = round - (m_p.Q * m_p.sampleRate / - (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo))); - - if (minNK == 0 || maxNK == 0) { - // most likely pathological parameters of some sort - cerr << "WARNING: CQKernel::generateKernel: minNK or maxNK is zero (minNK == " << minNK << ", maxNK == " << maxNK << "), not generating a kernel" << endl; - m_p.atomSpacing = 0; - m_p.firstCentre = 0; - m_p.fftSize = 0; - m_p.atomsPerFrame = 0; - m_p.lastCentre = 0; - m_p.fftHop = 0; - return; - } - - m_p.atomSpacing = round(minNK * atomHopFactor); - m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing); - m_p.fftSize = MathUtilities::nextPowerOfTwo - (m_p.firstCentre + ceil(maxNK / 2.0)); - - m_p.atomsPerFrame = floor - (1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing); - - cerr << "atomsPerFrame = " << m_p.atomsPerFrame << " (atomHopFactor = " << atomHopFactor << ", atomSpacing = " << m_p.atomSpacing << ", fftSize = " << m_p.fftSize << ", maxNK = " << maxNK << ", firstCentre = " << m_p.firstCentre << ")" << endl; - - m_p.lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing; - - m_p.fftHop = (m_p.lastCentre + m_p.atomSpacing) - m_p.firstCentre; - - cerr << "fftHop = " << m_p.fftHop << endl; - - m_fft = new FFT(m_p.fftSize); - - for (int k = 1; k <= m_p.binsPerOctave; ++k) { - - int nk = round(m_p.Q * m_p.sampleRate / - (m_p.minFrequency * pow(2, ((k-1.0) / bpo)))); - - // The MATLAB version uses a symmetric window, but our windows - // are periodic. A symmetric window of size N is a periodic - // one of size N-1 with the first element stuck on the end - Window<double> w(BlackmanHarrisWindow, nk-1); - vector<double> win = w.getWindowData(); - win.push_back(win[0]); - - for (int i = 0; i < (int)win.size(); ++i) { - win[i] = sqrt(win[i]) / nk; - } - - double fk = m_p.minFrequency * pow(2, ((k-1.0) / bpo)); - - vector<double> reals, imags; - - for (int i = 0; i < nk; ++i) { - double arg = (2.0 * M_PI * fk * i) / m_p.sampleRate; - reals.push_back(win[i] * cos(arg)); - imags.push_back(win[i] * sin(arg)); - } - - int atomOffset = m_p.firstCentre - int(ceil(nk/2.0)); - - for (int i = 0; i < m_p.atomsPerFrame; ++i) { - - int shift = atomOffset + (i * m_p.atomSpacing); - - vector<double> rin(m_p.fftSize, 0.0); - vector<double> iin(m_p.fftSize, 0.0); - - for (int j = 0; j < nk; ++j) { - rin[j + shift] = reals[j]; - iin[j + shift] = imags[j]; - } - - vector<double> rout(m_p.fftSize, 0.0); - vector<double> iout(m_p.fftSize, 0.0); - - m_fft->process(false, - rin.data(), iin.data(), - rout.data(), iout.data()); - - // Keep this dense for the moment (until after - // normalisation calculations) - - vector<C> row; - - for (int j = 0; j < m_p.fftSize; ++j) { - if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) { - row.push_back(C(0, 0)); - } else { - row.push_back(C(rout[j] / m_p.fftSize, - iout[j] / m_p.fftSize)); - } - } - - m_kernel.origin.push_back(0); - m_kernel.data.push_back(row); - } - } - - assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame); - - // print density as diagnostic - - int nnz = 0; - for (int i = 0; i < (int)m_kernel.data.size(); ++i) { - for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) { - if (m_kernel.data[i][j] != C(0, 0)) { - ++nnz; - } - } - } - - cerr << "size = " << m_kernel.data.size() << "*" << m_kernel.data[0].size() << " (fft size = " << m_p.fftSize << ")" << endl; - - assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame); - assert((int)m_kernel.data[0].size() == m_p.fftSize); - - cerr << "density = " << double(nnz) / double(m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize) << " (" << nnz << " of " << m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize << ")" << endl; - - finaliseKernel(); -} - -static bool ccomparator(C &c1, C &c2) -{ - return abs(c1) < abs(c2); -} - -static int maxidx(vector<C> &v) -{ - return std::max_element(v.begin(), v.end(), ccomparator) - v.begin(); -} - -void -CQKernel::finaliseKernel() -{ - // calculate weight for normalisation - - int wx1 = maxidx(m_kernel.data[0]); - int wx2 = maxidx(m_kernel.data[m_kernel.data.size()-1]); - - vector<vector<C> > subset(m_kernel.data.size()); - for (int j = wx1; j <= wx2; ++j) { - for (int i = 0; i < (int)m_kernel.data.size(); ++i) { - subset[i].push_back(m_kernel.data[i][j]); - } - } - - int nrows = subset.size(); - int ncols = subset[0].size(); - vector<vector<C> > square(ncols); // conjugate transpose of subset * subset - - for (int i = 0; i < nrows; ++i) { - assert((int)subset[i].size() == ncols); - } - - for (int j = 0; j < ncols; ++j) { - for (int i = 0; i < ncols; ++i) { - C v(0, 0); - for (int k = 0; k < nrows; ++k) { - v += subset[k][i] * conj(subset[k][j]); - } - square[i].push_back(v); - } - } - - vector<double> wK; - double q = 1; //!!! duplicated from constructor - for (int i = round(1.0/q); i < ncols - round(1.0/q) - 2; ++i) { - wK.push_back(abs(square[i][i])); - } - - double weight = double(m_p.fftHop) / m_p.fftSize; - weight /= MathUtilities::mean(wK.data(), wK.size()); - weight = sqrt(weight); - - cerr << "weight = " << weight << endl; - - // apply normalisation weight, make sparse, and store conjugate - // (we use the adjoint or conjugate transpose of the kernel matrix - // for the forward transform, the plain kernel for the inverse - // which we expect to be less common) - - KernelMatrix sk; - - for (int i = 0; i < (int)m_kernel.data.size(); ++i) { - - sk.origin.push_back(0); - sk.data.push_back(vector<C>()); - - int lastNZ = 0; - for (int j = (int)m_kernel.data[i].size()-1; j >= 0; --j) { - if (abs(m_kernel.data[i][j]) != 0.0) { - lastNZ = j; - break; - } - } - - bool haveNZ = false; - for (int j = 0; j <= lastNZ; ++j) { - if (haveNZ || abs(m_kernel.data[i][j]) != 0.0) { - if (!haveNZ) sk.origin[i] = j; - haveNZ = true; - sk.data[i].push_back(conj(m_kernel.data[i][j]) * weight); - } - } - } - - m_kernel = sk; -} - -vector<C> -CQKernel::processForward(const vector<C> &cv) -{ - // straightforward matrix multiply (taking into account m_kernel's - // slightly-sparse representation) - - if (m_kernel.data.empty()) return vector<C>(); - - int nrows = m_p.binsPerOctave * m_p.atomsPerFrame; - - vector<C> rv(nrows, C()); - - for (int i = 0; i < nrows; ++i) { - int len = m_kernel.data[i].size(); - for (int j = 0; j < len; ++j) { - rv[i] += cv[j + m_kernel.origin[i]] * m_kernel.data[i][j]; - } - } - - return rv; -} - -vector<C> -CQKernel::processInverse(const vector<C> &cv) -{ - // matrix multiply by conjugate transpose of m_kernel. This is - // actually the original kernel as calculated, we just stored the - // conjugate-transpose of the kernel because we expect to be doing - // more forward transforms than inverse ones. - - if (m_kernel.data.empty()) return vector<C>(); - - int ncols = m_p.binsPerOctave * m_p.atomsPerFrame; - int nrows = m_p.fftSize; - - vector<C> rv(nrows, C()); - - for (int j = 0; j < ncols; ++j) { - int i0 = m_kernel.origin[j]; - int i1 = i0 + m_kernel.data[j].size(); - for (int i = i0; i < i1; ++i) { - rv[i] += cv[j] * conj(m_kernel.data[j][i - i0]); - } - } - - return rv; -} - -
--- a/cpp-qm-dsp/CQKernel.h Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,82 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#ifndef CQ_KERNEL_H -#define CQ_KERNEL_H - -#include <vector> -#include <complex> - -class FFT; - -class CQKernel -{ -public: - CQKernel(double sampleRate, double maxFreq, int binsPerOctave); - ~CQKernel(); - - struct Properties { - double sampleRate; - double maxFrequency; - double minFrequency; - int binsPerOctave; - int fftSize; - int fftHop; - int atomsPerFrame; - int atomSpacing; - int firstCentre; - int lastCentre; - double Q; - }; - - Properties getProperties() const { return m_p; } - - std::vector<std::complex<double> > processForward - (const std::vector<std::complex<double> > &); - - std::vector<std::complex<double> > processInverse - (const std::vector<std::complex<double> > &); - -private: - Properties m_p; - FFT *m_fft; - - struct KernelMatrix { - std::vector<int> origin; - std::vector<std::vector<std::complex<double> > > data; - }; - KernelMatrix m_kernel; - - void generateKernel(); - void finaliseKernel(); -}; - -#endif
--- a/cpp-qm-dsp/CQSpectrogram.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,261 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#include "CQSpectrogram.h" - -#include <iostream> -#include <stdexcept> - -using std::cerr; -using std::endl; - -CQSpectrogram::CQSpectrogram(double sampleRate, - double minFreq, double maxFreq, - int binsPerOctave, - Interpolation interpolation) : - m_cq(sampleRate, minFreq, maxFreq, binsPerOctave), - m_interpolation(interpolation) -{ -} - -CQSpectrogram::~CQSpectrogram() -{ -} - -CQSpectrogram::RealBlock -CQSpectrogram::process(const RealSequence &td) -{ - return postProcess(m_cq.process(td), false); -} - -CQSpectrogram::RealBlock -CQSpectrogram::getRemainingOutput() -{ - return postProcess(m_cq.getRemainingOutput(), true); -} - -CQSpectrogram::RealBlock -CQSpectrogram::postProcess(const ComplexBlock &cq, bool insist) -{ - int width = cq.size(); - - // convert to magnitudes - RealBlock spec; - for (int i = 0; i < width; ++i) { - int height = cq[i].size(); - RealColumn col(height, 0); - for (int j = 0; j < height; ++j) { - col[j] = abs(cq[i][j]); - } - spec.push_back(col); - } - - if (m_interpolation == InterpolateZeros) { - for (int i = 0; i < width; ++i) { - int sh = spec[i].size(); - int fh = getTotalBins(); - for (int j = sh; j < fh; ++j) { - spec[i].push_back(0); - } - } - return spec; - } - - for (int i = 0; i < width; ++i) { - m_buffer.push_back(spec[i]); - } - - if (m_interpolation == InterpolateHold) { - return fetchHold(insist); - } else { - return fetchLinear(insist); - } -} - -CQSpectrogram::RealBlock -CQSpectrogram::fetchHold(bool) -{ - RealBlock out; - - int width = m_buffer.size(); - int height = getTotalBins(); - - for (int i = 0; i < width; ++i) { - - RealColumn col = m_buffer[i]; - - int thisHeight = col.size(); - int prevHeight = m_prevColumn.size(); - - for (int j = thisHeight; j < height; ++j) { - if (j < prevHeight) { - col.push_back(m_prevColumn[j]); - } else { - col.push_back(0.0); - } - } - - m_prevColumn = col; - out.push_back(col); - } - - m_buffer.clear(); - - return out; -} - -CQSpectrogram::RealBlock -CQSpectrogram::fetchLinear(bool insist) -{ - RealBlock out; - - //!!! This is surprisingly messy. I must be missing something. - - // We can only return any data when we have at least one column - // that has the full height in the buffer, that is not the first - // column. - // - // If the first col has full height, and there is another one - // later that also does, then we can interpolate between those, up - // to but not including the second full height column. Then we - // drop and return the columns we interpolated, leaving the second - // full-height col as the first col in the buffer. And repeat as - // long as enough columns are available. - // - // If the first col does not have full height, then (so long as - // we're following the logic above) we must simply have not yet - // reached the first full-height column in the CQ output, and we - // can interpolate nothing. - - int width = m_buffer.size(); - int height = getTotalBins(); - - if (width == 0) return out; - - int firstFullHeight = -1; - int secondFullHeight = -1; - - for (int i = 0; i < width; ++i) { - if ((int)m_buffer[i].size() == height) { - if (firstFullHeight == -1) { - firstFullHeight = i; - } else if (secondFullHeight == -1) { - secondFullHeight = i; - break; - } - } - } - -// cerr << "fetchLinear: firstFullHeight = " << firstFullHeight << ", secondFullHeight = " << secondFullHeight << endl; - - if (firstFullHeight < 0) { - if (insist) { - return fetchHold(true); - } else { - return out; - } - } else if (firstFullHeight > 0) { - // can interpolate nothing, stash up to first full height & recurse - out = RealBlock(m_buffer.begin(), m_buffer.begin() + firstFullHeight); - m_buffer = RealBlock(m_buffer.begin() + firstFullHeight, m_buffer.end()); - RealBlock more = fetchLinear(insist); - out.insert(out.end(), more.begin(), more.end()); - return out; - } else if (secondFullHeight < 0) { - // firstFullHeight == 0, but there is no second full height -- - // wait for it unless insist flag is set - if (insist) { - return fetchHold(true); - } else { - return out; - } - } else { - // firstFullHeight == 0 and secondFullHeight also valid. Can interpolate - out = linearInterpolated(m_buffer, 0, secondFullHeight); - m_buffer = RealBlock(m_buffer.begin() + secondFullHeight, m_buffer.end()); - RealBlock more = fetchLinear(insist); - out.insert(out.end(), more.begin(), more.end()); - return out; - } -} - -CQSpectrogram::RealBlock -CQSpectrogram::linearInterpolated(const RealBlock &g, int x0, int x1) -{ - // g must be a grid with full-height columns at x0 and x1 - - if (x0 >= x1) { - throw std::logic_error("x0 >= x1"); - } - if (x1 >= (int)g.size()) { - throw std::logic_error("x1 >= g.size()"); - } - if (g[x0].size() != g[x1].size()) { - throw std::logic_error("x0 and x1 are not the same height"); - } - - int height = g[x0].size(); - int width = x1 - x0; - - RealBlock out(g.begin() + x0, g.begin() + x1); - - for (int y = 0; y < height; ++y) { - - int spacing = width; - for (int i = 1; i < width; ++i) { - int thisHeight = g[x0 + i].size(); - if (thisHeight > height) { - throw std::logic_error("First column not full-height"); - } - if (thisHeight > y) { - spacing = i; - break; - } - } - - if (spacing < 2) continue; - - for (int i = 0; i + spacing <= width; i += spacing) { - for (int j = 1; j < spacing; ++j) { - double proportion = double(j)/double(spacing); - double interpolated = - g[x0 + i][y] * (1.0 - proportion) + - g[x0 + i + spacing][y] * proportion; - out[i + j].push_back(interpolated); - } - } - } - - return out; -} - - -
--- a/cpp-qm-dsp/CQSpectrogram.h Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,77 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#ifndef CQSPECTROGRAM_H -#define CQSPECTROGRAM_H - -#include "ConstantQ.h" - -class CQSpectrogram : public CQBase -{ -public: - enum Interpolation { - InterpolateZeros, // leave empty cells as zero - InterpolateHold, // repeat prior cell - InterpolateLinear, // linear interpolation between consecutive time cells - }; - - CQSpectrogram(double sampleRate, - double minFreq, double maxFreq, - int binsPerOctave, - Interpolation interpolation); - virtual ~CQSpectrogram(); - - virtual double getSampleRate() const { return m_cq.getSampleRate(); } - virtual int getBinsPerOctave() const { return m_cq.getBinsPerOctave(); } - virtual int getOctaves() const { return m_cq.getOctaves(); } - virtual int getTotalBins() const { return m_cq.getTotalBins(); } - virtual int getColumnHop() const { return m_cq.getColumnHop(); } - virtual int getLatency() const { return m_cq.getLatency(); } - virtual double getMaxFrequency() const { return m_cq.getMaxFrequency(); } - virtual double getMinFrequency() const { return m_cq.getMinFrequency(); } - virtual double getBinFrequency(int bin) const { return m_cq.getBinFrequency(bin); } - - RealBlock process(const RealSequence &); - RealBlock getRemainingOutput(); - -private: - ConstantQ m_cq; - Interpolation m_interpolation; - - RealBlock m_buffer; - RealBlock postProcess(const ComplexBlock &, bool insist); - RealBlock fetchHold(bool insist); - RealBlock fetchLinear(bool insist); - RealBlock linearInterpolated(const RealBlock &, int, int); - RealColumn m_prevColumn; -}; - -#endif
--- a/cpp-qm-dsp/ConstantQ.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,352 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#include "ConstantQ.h" - -#include "CQKernel.h" - -#include "dsp/rateconversion/Resampler.h" -#include "maths/MathUtilities.h" -#include "dsp/transforms/FFT.h" - -#include <algorithm> -#include <iostream> -#include <stdexcept> - -using std::vector; -using std::cerr; -using std::endl; - -//#define DEBUG_CQ 1 - -ConstantQ::ConstantQ(double sampleRate, - double minFreq, - double maxFreq, - int binsPerOctave) : - m_sampleRate(sampleRate), - m_maxFrequency(maxFreq), - m_minFrequency(minFreq), - m_binsPerOctave(binsPerOctave), - m_fft(0) -{ - if (minFreq <= 0.0 || maxFreq <= 0.0) { - throw std::invalid_argument("Frequency extents must be positive"); - } - - initialise(); -} - -ConstantQ::~ConstantQ() -{ - delete m_fft; - for (int i = 0; i < (int)m_decimators.size(); ++i) { - delete m_decimators[i]; - } - delete m_kernel; -} - -double -ConstantQ::getMinFrequency() const -{ - return m_p.minFrequency / pow(2.0, m_octaves - 1); -} - -double -ConstantQ::getBinFrequency(int bin) const -{ - return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave())); -} - -void -ConstantQ::initialise() -{ - m_octaves = int(ceil(log2(m_maxFrequency / m_minFrequency))); - m_kernel = new CQKernel(m_sampleRate, m_maxFrequency, m_binsPerOctave); - m_p = m_kernel->getProperties(); - - // Use exact powers of two for resampling rates. They don't have - // to be related to our actual samplerate: the resampler only - // cares about the ratio, but it only accepts integer source and - // target rates, and if we start from the actual samplerate we - // risk getting non-integer rates for lower octaves - - int sourceRate = pow(2, m_octaves); - vector<int> latencies; - - // top octave, no resampling - latencies.push_back(0); - m_decimators.push_back(0); - - for (int i = 1; i < m_octaves; ++i) { - - int factor = pow(2, i); - - Resampler *r = new Resampler - (sourceRate, sourceRate / factor, 50, 0.05); - -#ifdef DEBUG_CQ - cerr << "forward: octave " << i << ": resample from " << sourceRate << " to " << sourceRate / factor << endl; -#endif - - // We need to adapt the latencies so as to get the first input - // sample to be aligned, in time, at the decimator output - // across all octaves. - // - // Our decimator uses a linear phase filter, but being causal - // it is not zero phase: it has a latency that depends on the - // decimation factor. Those latencies have been calculated - // per-octave and are available to us in the latencies - // array. Left to its own devices, the first input sample will - // appear at output sample 0 in the highest octave (where no - // decimation is needed), sample number latencies[1] in the - // next octave down, latencies[2] in the next one, etc. We get - // to apply some artificial per-octave latency after the - // decimator in the processing chain, in order to compensate - // for the differing latencies associated with different - // decimation factors. How much should we insert? - // - // The outputs of the decimators are at different rates (in - // terms of the relation between clock time and samples) and - // we want them aligned in terms of time. So, for example, a - // latency of 10 samples with a decimation factor of 2 is - // equivalent to a latency of 20 with no decimation -- they - // both result in the first output sample happening at the - // same equivalent time in milliseconds. - // - // So here we record the latency added by the decimator, in - // terms of the sample rate of the undecimated signal. Then we - // use that to compensate in a moment, when we've discovered - // what the longest latency across all octaves is. - - latencies.push_back(r->getLatency() * factor); - m_decimators.push_back(r); - } - - m_bigBlockSize = m_p.fftSize * pow(2, m_octaves - 1); - - // Now add in the extra padding and compensate for hops that must - // be dropped in order to align the atom centres across - // octaves. Again this is a bit trickier because we are doing it - // at input rather than output and so must work in per-octave - // sample rates rather than output blocks - - int emptyHops = m_p.firstCentre / m_p.atomSpacing; - - vector<int> drops; - for (int i = 0; i < m_octaves; ++i) { - int factor = pow(2, i); - int dropHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops; - int drop = ((dropHops * m_p.fftHop) * factor) / m_p.atomsPerFrame; - drops.push_back(drop); - } - - int maxLatPlusDrop = 0; - for (int i = 0; i < m_octaves; ++i) { - int latPlusDrop = latencies[i] + drops[i]; - if (latPlusDrop > maxLatPlusDrop) maxLatPlusDrop = latPlusDrop; - } - - int totalLatency = maxLatPlusDrop; - - int lat0 = totalLatency - latencies[0] - drops[0]; - totalLatency = ceil(double(lat0 / m_p.fftHop) * m_p.fftHop) - + latencies[0] + drops[0]; - - // We want (totalLatency - latencies[i]) to be a multiple of 2^i - // for each octave i, so that we do not end up with fractional - // octave latencies below. In theory this is hard, in practice if - // we ensure it for the last octave we should be OK. - double finalOctLat = latencies[m_octaves-1]; - double finalOctFact = pow(2, m_octaves-1); - totalLatency = - int(round(finalOctLat + - finalOctFact * - ceil((totalLatency - finalOctLat) / finalOctFact))); - -#ifdef DEBUG_CQ - cerr << "total latency = " << totalLatency << endl; -#endif - - // Padding as in the reference (will be introduced with the - // latency compensation in the loop below) - m_outputLatency = totalLatency + m_bigBlockSize - - m_p.firstCentre * pow(2, m_octaves-1); - -#ifdef DEBUG_CQ - cerr << "m_bigBlockSize = " << m_bigBlockSize << ", firstCentre = " - << m_p.firstCentre << ", m_octaves = " << m_octaves - << ", so m_outputLatency = " << m_outputLatency << endl; -#endif - - for (int i = 0; i < m_octaves; ++i) { - - double factor = pow(2, i); - - // Calculate the difference between the total latency applied - // across all octaves, and the existing latency due to the - // decimator for this octave, and then convert it back into - // the sample rate appropriate for the output latency of this - // decimator -- including one additional big block of padding - // (as in the reference). - - double octaveLatency = - double(totalLatency - latencies[i] - drops[i] - + m_bigBlockSize) / factor; - -#ifdef DEBUG_CQ - cerr << "octave " << i << ": resampler latency = " << latencies[i] - << ", drop " << drops[i] << " (/factor = " << drops[i]/factor - << "), octaveLatency = " << octaveLatency << " -> " - << int(round(octaveLatency)) << " (diff * factor = " - << (octaveLatency - round(octaveLatency)) << " * " - << factor << " = " - << (octaveLatency - round(octaveLatency)) * factor << ")" << endl; - - cerr << "double(" << totalLatency << " - " - << latencies[i] << " - " << drops[i] << " + " - << m_bigBlockSize << ") / " << factor << " = " - << octaveLatency << endl; -#endif - - m_buffers.push_back - (RealSequence(int(round(octaveLatency)), 0.0)); - } - - m_fft = new FFTReal(m_p.fftSize); -} - -ConstantQ::ComplexBlock -ConstantQ::process(const RealSequence &td) -{ - m_buffers[0].insert(m_buffers[0].end(), td.begin(), td.end()); - - for (int i = 1; i < m_octaves; ++i) { - RealSequence dec = m_decimators[i]->process(td.data(), td.size()); - m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end()); - } - - ComplexBlock out; - - while (true) { - - // We could have quite different remaining sample counts in - // different octaves, because (apart from the predictable - // added counts for decimator output on each block) we also - // have variable additional latency per octave - bool enough = true; - for (int i = 0; i < m_octaves; ++i) { - int required = m_p.fftSize * pow(2, m_octaves - i - 1); - if ((int)m_buffers[i].size() < required) { - enough = false; - } - } - if (!enough) break; - - int base = out.size(); - int totalColumns = pow(2, m_octaves - 1) * m_p.atomsPerFrame; - for (int i = 0; i < totalColumns; ++i) { - out.push_back(ComplexColumn()); - } - - for (int octave = 0; octave < m_octaves; ++octave) { - - int blocksThisOctave = pow(2, (m_octaves - octave - 1)); - - for (int b = 0; b < blocksThisOctave; ++b) { - ComplexBlock block = processOctaveBlock(octave); - - for (int j = 0; j < m_p.atomsPerFrame; ++j) { - - int target = base + - (b * (totalColumns / blocksThisOctave) + - (j * ((totalColumns / blocksThisOctave) / - m_p.atomsPerFrame))); - - while (int(out[target].size()) < - m_p.binsPerOctave * (octave + 1)) { - out[target].push_back(Complex()); - } - - for (int i = 0; i < m_p.binsPerOctave; ++i) { - out[target][m_p.binsPerOctave * octave + i] = - block[j][m_p.binsPerOctave - i - 1]; - } - } - } - } - } - - return out; -} - -ConstantQ::ComplexBlock -ConstantQ::getRemainingOutput() -{ - // Same as padding added at start, though rounded up - int pad = ceil(double(m_outputLatency) / m_bigBlockSize) * m_bigBlockSize; - RealSequence zeros(pad, 0.0); - return process(zeros); -} - -ConstantQ::ComplexBlock -ConstantQ::processOctaveBlock(int octave) -{ - RealSequence ro(m_p.fftSize, 0.0); - RealSequence io(m_p.fftSize, 0.0); - - m_fft->forward(m_buffers[octave].data(), ro.data(), io.data()); - - RealSequence shifted; - shifted.insert(shifted.end(), - m_buffers[octave].begin() + m_p.fftHop, - m_buffers[octave].end()); - m_buffers[octave] = shifted; - - ComplexSequence cv; - for (int i = 0; i < m_p.fftSize; ++i) { - cv.push_back(Complex(ro[i], io[i])); - } - - ComplexSequence cqrowvec = m_kernel->processForward(cv); - - // Reform into a column matrix - ComplexBlock cqblock; - for (int j = 0; j < m_p.atomsPerFrame; ++j) { - cqblock.push_back(ComplexColumn()); - for (int i = 0; i < m_p.binsPerOctave; ++i) { - cqblock[j].push_back(cqrowvec[i * m_p.atomsPerFrame + j]); - } - } - - return cqblock; -} - -
--- a/cpp-qm-dsp/ConstantQ.h Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,117 +0,0 @@ -/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#ifndef CONSTANTQ_H -#define CONSTANTQ_H - -#include "CQBase.h" -#include "CQKernel.h" - -class Resampler; -class FFTReal; - -/** - * Calculate a complex sparse constant-Q representation from - * time-domain input. - * - * For a real (magnitude-only) or interpolated representation, see - * CQSpectrogram. - */ -class ConstantQ : public CQBase -{ -public: - ConstantQ(double sampleRate, - double minFreq, double maxFreq, - int binsPerOctave); - virtual ~ConstantQ(); - - virtual double getSampleRate() const { return m_sampleRate; } - virtual int getBinsPerOctave() const { return m_binsPerOctave; } - virtual int getOctaves() const { return m_octaves; } - virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; } - virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; } - virtual int getLatency() const { return m_outputLatency; } - virtual double getMaxFrequency() const { return m_p.maxFrequency; } - virtual double getMinFrequency() const; - virtual double getBinFrequency(int bin) const; - - /** - * Given a series of time-domain samples, return a series of - * constant-Q columns. Any samples left over (that did not fit - * into a constant-Q processing block) are saved for the next call - * to process or getRemainingBlocks. - * - * Each column contains a series of constant-Q bin values ordered - * from highest to lowest frequency. - * - * Columns are of variable height: each will contain at least - * getBinsPerOctave() values, because the highest-frequency octave - * is always present, but a second octave (if requested) will - * appear only in alternate columns, a third octave only in every - * fourth column, and so on. - * - * If you need a format in which all columns are of equal height - * and every bin contains a value, use CQInterpolated instead of - * ConstantQ. - */ - ComplexBlock process(const RealSequence &); - - /** - * Return the remaining constant-Q columns following the end of - * processing. Any buffered input is padded so as to ensure that - * all input provided to process() will have been returned. - */ - ComplexBlock getRemainingOutput(); - -private: - double m_sampleRate; - double m_maxFrequency; - double m_minFrequency; - int m_binsPerOctave; - int m_octaves; - - CQKernel *m_kernel; - CQKernel::Properties m_p; - int m_bigBlockSize; - - std::vector<Resampler *> m_decimators; - std::vector<RealSequence> m_buffers; - - int m_outputLatency; - - FFTReal *m_fft; - - void initialise(); - ComplexBlock processOctaveBlock(int octave); -}; - -#endif -
--- a/cpp-qm-dsp/processfile.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,294 +0,0 @@ - -#include "ConstantQ.h" -#include "CQInverse.h" - -#include <sndfile.h> - -#include <iostream> - -using std::vector; -using std::cerr; -using std::endl; - -#include <cstring> - -#include <getopt.h> -#include <unistd.h> -#include <sys/time.h> -#include <cstdlib> - -int main(int argc, char **argv) -{ - double maxFreq = 0; - double minFreq = 0; - int bpo = 0; - bool help = false; - - int c; - - while (1) { - int optionIndex = 0; - - static struct option longOpts[] = { - { "help", 0, 0, 'h', }, - { "maxfreq", 1, 0, 'x', }, - { "minfreq", 1, 0, 'n', }, - { "bpo", 1, 0, 'b' }, - { 0, 0, 0, 0 }, - }; - - c = getopt_long(argc, argv, - "hx:n:b:", - longOpts, &optionIndex); - if (c == -1) break; - - switch (c) { - case 'h': help = true; break; - case 'x': maxFreq = atof(optarg); break; - case 'n': minFreq = atof(optarg); break; - case 'b': bpo = atoi(optarg); break; - default: help = true; break; - } - } - - if (help || (optind + 2 != argc && optind + 3 != argc)) { - cerr << endl; - cerr << "Usage: " << argv[0] << " [options] infile.wav outfile.wav [differencefile.wav]" << endl; - cerr << endl; - cerr << "Options:" << endl; - cerr << " -x<X>, --maxfreq <X> Maximum frequency (default = sample rate / 3)" << endl; - cerr << " -n<X>, --minfreq <X> Minimum frequency (default = 100, actual min may vary)" << endl; - cerr << " -b<X>, --bpo <X> Bins per octave (default = 60)" << endl; - cerr << " -h, --help Print this help" << endl; - cerr << endl; - cerr << "This rather useless program simply performs a forward Constant-Q transform with" << endl; - cerr << "the requested parameters, followed by its inverse, and writes the result to the" << endl; - cerr << "output file. If a diff file name is provided, the difference between input and" << endl; - cerr << "output signals is also written to that. All this accomplishes is to produce a" << endl; - cerr << "signal that approximates the input: it's intended for test purposes only." << endl; - cerr << endl; - cerr << "(Want to calculate and obtain a Constant-Q spectrogram? Use the CQVamp plugin" << endl; - cerr << "in a Vamp plugin host.)" << endl; - cerr << endl; - return 2; - } - - char *fileName = strdup(argv[optind++]); - char *fileNameOut = strdup(argv[optind++]); - char *diffFileName = (optind < argc ? strdup(argv[optind++]) : 0); - bool doDiff = (diffFileName != 0); - - SNDFILE *sndfile; - SNDFILE *sndfileOut; - SNDFILE *sndDiffFile = 0; - SF_INFO sfinfo; - SF_INFO sfinfoOut; - SF_INFO sfinfoDiff; - 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; - } - - sfinfoOut.channels = 1; - sfinfoOut.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16; - sfinfoOut.frames = sfinfo.frames; - sfinfoOut.samplerate = sfinfo.samplerate; - 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; - } - - if (doDiff) { - sfinfoDiff = sfinfoOut; - sndDiffFile = sf_open(diffFileName, SFM_WRITE, &sfinfoDiff); - if (!sndDiffFile) { - cerr << "ERROR: Failed to open diff output file \"" << diffFileName << "\" for writing: " - << sf_strerror(sndDiffFile) << endl; - return 1; - } - } - - int ibs = 1024; - int channels = sfinfo.channels; - float *fbuf = new float[channels * ibs]; - - if (maxFreq == 0.0) maxFreq = sfinfo.samplerate / 3; - if (minFreq == 0.0) minFreq = 100; - if (bpo == 0) bpo = 60; - - ConstantQ cq(sfinfo.samplerate, minFreq, maxFreq, bpo); - CQInverse cqi(sfinfo.samplerate, minFreq, maxFreq, bpo); - - cerr << "max freq = " << cq.getMaxFrequency() << ", min freq = " - << cq.getMinFrequency() << ", octaves = " << cq.getOctaves() << endl; - - cerr << "octave boundaries: "; - for (int i = 0; i < cq.getOctaves(); ++i) { - cerr << cq.getMaxFrequency() / pow(2, i) << " "; - } - cerr << endl; - - int inframe = 0; - int outframe = 0; - int latency = cq.getLatency() + cqi.getLatency(); - - vector<double> buffer; - - double maxdiff = 0.0; - int maxdiffidx = 0; - - cerr << "forward latency = " << cq.getLatency() << ", inverse latency = " - << cqi.getLatency() << ", total = " << latency << endl; - - timeval tv; - (void)gettimeofday(&tv, 0); - - while (inframe < sfinfo.frames) { - - int count = -1; - - if ((count = sf_readf_float(sndfile, fbuf, ibs)) < 0) { - break; - } - - vector<double> cqin; - for (int i = 0; i < count; ++i) { - double v = fbuf[i * channels]; - if (channels > 1) { - for (int c = 1; c < channels; ++c) { - v += fbuf[i * channels + c]; - } - v /= channels; - } - cqin.push_back(v); - } - - if (doDiff) { - buffer.insert(buffer.end(), cqin.begin(), cqin.end()); - } - - vector<double> cqout = cqi.process(cq.process(cqin)); - - for (int i = 0; i < int(cqout.size()); ++i) { - if (cqout[i] > 1.0) cqout[i] = 1.0; - if (cqout[i] < -1.0) cqout[i] = -1.0; - } - - if (outframe >= latency) { - - sf_writef_double(sndfileOut, - cqout.data(), - cqout.size()); - - } else if (outframe + (int)cqout.size() >= latency) { - - int offset = latency - outframe; - sf_writef_double(sndfileOut, - cqout.data() + offset, - cqout.size() - offset); - } - - if (doDiff) { - for (int i = 0; i < (int)cqout.size(); ++i) { - if (outframe + i >= latency) { - int dframe = outframe + i - latency; - if (dframe >= (int)buffer.size()) cqout[i] = 0; - else cqout[i] -= buffer[dframe]; - if (fabs(cqout[i]) > maxdiff && - dframe > sfinfo.samplerate && // ignore first/last sec - dframe + sfinfo.samplerate < sfinfo.frames) { - maxdiff = fabs(cqout[i]); - maxdiffidx = dframe; - } - } - } - - if (outframe >= latency) { - - sf_writef_double(sndDiffFile, - cqout.data(), - cqout.size()); - - } else if (outframe + (int)cqout.size() >= latency) { - - int offset = latency - outframe; - sf_writef_double(sndDiffFile, - cqout.data() + offset, - cqout.size() - offset); - } - } - - inframe += count; - outframe += cqout.size(); - } - - vector<double> r = cqi.process(cq.getRemainingOutput()); - vector<double> r2 = cqi.getRemainingOutput(); - - r.insert(r.end(), r2.begin(), r2.end()); - - for (int i = 0; i < int(r.size()); ++i) { - if (r[i] > 1.0) r[i] = 1.0; - if (r[i] < -1.0) r[i] = -1.0; - } - - sf_writef_double(sndfileOut, r.data(), r.size()); - if (doDiff) { - for (int i = 0; i < (int)r.size(); ++i) { - if (outframe + i >= latency) { - int dframe = outframe + i - latency; - if (dframe >= (int)buffer.size()) r[i] = 0; - else r[i] -= buffer[dframe]; - if (fabs(r[i]) > maxdiff && - dframe > sfinfo.samplerate && // ignore first/last sec - dframe + sfinfo.samplerate < sfinfo.frames) { - maxdiff = fabs(r[i]); - maxdiffidx = dframe; - } - } - } - sf_writef_double(sndDiffFile, r.data(), r.size()); - } - outframe += r.size(); - - sf_close(sndfile); - sf_close(sndfileOut); - - if (doDiff) { - sf_close(sndDiffFile); - } - - cerr << "in: " << inframe << ", out: " << outframe - latency << endl; - - if (doDiff) { - double db = 10 * log10(maxdiff); - cerr << "max diff [excluding first and last second of audio] is " - << maxdiff << " (" << db << " dBFS)" - << " at sample index " << maxdiffidx << endl; - } - - 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 << "elapsed time (not counting init): " << sec << " sec, frames/sec at input: " << inframe/sec << endl; - - return 0; -} -
--- a/cpp-qm-dsp/test.cpp Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -#include "CQSpectrogram.h" - -#include <iostream> -#include <vector> - -using std::vector; -using std::cerr; -using std::cout; -using std::endl; - -#include <cstdio> - -int main(int argc, char **argv) -{ - vector<double> in; - - for (int i = 0; i < 64; ++i) { -// if (i == 0) in.push_back(1); -// else in.push_back(0); - in.push_back(sin(i * M_PI / 2.0)); - } - - CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros); - - vector<vector<double> > out = k.process(in); - vector<vector<double> > rest = k.getRemainingOutput(); - - out.insert(out.end(), rest.begin(), rest.end()); - - cerr << "got " << out.size() << " back (" << out[0].size() << " in each?)" << endl; - - for (int b = 0; b < (int)out.size() / 8; ++b) { - printf("\nColumns %d to %d:\n\n", b * 8, b * 8 + 7); - for (int j = int(out[0].size()) - 1; j >= 0; --j) { - for (int i = 0; i < 8; ++i) { - if (i + b * 8 < (int)out.size()) { - double v = out[i + b * 8][j]; - if (v < 0.0001) printf(" 0 "); - else printf(" %.4f ", out[i + b * 8][j]); - } - } - printf("\n"); - } - } -} -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cq/CQBase.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,60 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef CQBASE_H +#define CQBASE_H + +#include <vector> +#include <complex> + +class CQBase // interface class +{ +public: + typedef std::complex<double> Complex; + typedef std::vector<double> RealSequence; + typedef std::vector<double> RealColumn; + typedef std::vector<Complex> ComplexSequence; + typedef std::vector<Complex> ComplexColumn; + typedef std::vector<RealColumn> RealBlock; + typedef std::vector<ComplexColumn> ComplexBlock; + + virtual double getSampleRate() const = 0; + virtual int getBinsPerOctave() const = 0; + virtual int getOctaves() const = 0; + virtual int getTotalBins() const = 0; + virtual int getColumnHop() const = 0; + virtual int getLatency() const = 0; + virtual double getMaxFrequency() const = 0; + virtual double getMinFrequency() const = 0; // actual min, not that passed to ctor + virtual double getBinFrequency(int bin) const = 0; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cq/CQInverse.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,91 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef CQINVERSE_H +#define CQINVERSE_H + +#include "CQBase.h" +#include "CQKernel.h" + +class Resampler; +class FFTReal; + +class CQInverse : public CQBase +{ +public: + CQInverse(double sampleRate, + double minFreq, double maxFreq, + int binsPerOctave); + virtual ~CQInverse(); + + virtual double getSampleRate() const { return m_sampleRate; } + virtual int getBinsPerOctave() const { return m_binsPerOctave; } + virtual int getOctaves() const { return m_octaves; } + virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; } + virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; } + virtual int getLatency() const { return m_outputLatency; } + virtual double getMaxFrequency() const { return m_p.maxFrequency; } + virtual double getMinFrequency() const; // actual min, not that passed to ctor + virtual double getBinFrequency(int bin) const; + + // Input is the format produced by ConstantQ class, + // i.e. uninterpolated complex, not the real-valued stuff produced + // by CQSpectrogram + + RealSequence process(const ComplexBlock &); + RealSequence getRemainingOutput(); + +private: + double m_sampleRate; + double m_maxFrequency; + double m_minFrequency; + int m_binsPerOctave; + int m_octaves; + + CQKernel *m_kernel; + CQKernel::Properties m_p; + + std::vector<Resampler *> m_upsamplers; + std::vector<RealSequence> m_buffers; + std::vector<RealSequence> m_olaBufs; // fixed-length, for overlap-add + + int m_outputLatency; + + FFTReal *m_fft; + + void initialise(); + void processOctave(int octave, const ComplexBlock &block); + void processOctaveColumn(int octave, const ComplexColumn &column); + void overlapAddAndResample(int octave, const RealSequence &); + RealSequence drawFromBuffers(); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cq/CQSpectrogram.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,77 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef CQSPECTROGRAM_H +#define CQSPECTROGRAM_H + +#include "ConstantQ.h" + +class CQSpectrogram : public CQBase +{ +public: + enum Interpolation { + InterpolateZeros, // leave empty cells as zero + InterpolateHold, // repeat prior cell + InterpolateLinear, // linear interpolation between consecutive time cells + }; + + CQSpectrogram(double sampleRate, + double minFreq, double maxFreq, + int binsPerOctave, + Interpolation interpolation); + virtual ~CQSpectrogram(); + + virtual double getSampleRate() const { return m_cq.getSampleRate(); } + virtual int getBinsPerOctave() const { return m_cq.getBinsPerOctave(); } + virtual int getOctaves() const { return m_cq.getOctaves(); } + virtual int getTotalBins() const { return m_cq.getTotalBins(); } + virtual int getColumnHop() const { return m_cq.getColumnHop(); } + virtual int getLatency() const { return m_cq.getLatency(); } + virtual double getMaxFrequency() const { return m_cq.getMaxFrequency(); } + virtual double getMinFrequency() const { return m_cq.getMinFrequency(); } + virtual double getBinFrequency(int bin) const { return m_cq.getBinFrequency(bin); } + + RealBlock process(const RealSequence &); + RealBlock getRemainingOutput(); + +private: + ConstantQ m_cq; + Interpolation m_interpolation; + + RealBlock m_buffer; + RealBlock postProcess(const ComplexBlock &, bool insist); + RealBlock fetchHold(bool insist); + RealBlock fetchLinear(bool insist); + RealBlock linearInterpolated(const RealBlock &, int, int); + RealColumn m_prevColumn; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cq/ConstantQ.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,117 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef CONSTANTQ_H +#define CONSTANTQ_H + +#include "CQBase.h" +#include "CQKernel.h" + +class Resampler; +class FFTReal; + +/** + * Calculate a complex sparse constant-Q representation from + * time-domain input. + * + * For a real (magnitude-only) or interpolated representation, see + * CQSpectrogram. + */ +class ConstantQ : public CQBase +{ +public: + ConstantQ(double sampleRate, + double minFreq, double maxFreq, + int binsPerOctave); + virtual ~ConstantQ(); + + virtual double getSampleRate() const { return m_sampleRate; } + virtual int getBinsPerOctave() const { return m_binsPerOctave; } + virtual int getOctaves() const { return m_octaves; } + virtual int getTotalBins() const { return m_octaves * m_binsPerOctave; } + virtual int getColumnHop() const { return m_p.fftHop / m_p.atomsPerFrame; } + virtual int getLatency() const { return m_outputLatency; } + virtual double getMaxFrequency() const { return m_p.maxFrequency; } + virtual double getMinFrequency() const; + virtual double getBinFrequency(int bin) const; + + /** + * Given a series of time-domain samples, return a series of + * constant-Q columns. Any samples left over (that did not fit + * into a constant-Q processing block) are saved for the next call + * to process or getRemainingBlocks. + * + * Each column contains a series of constant-Q bin values ordered + * from highest to lowest frequency. + * + * Columns are of variable height: each will contain at least + * getBinsPerOctave() values, because the highest-frequency octave + * is always present, but a second octave (if requested) will + * appear only in alternate columns, a third octave only in every + * fourth column, and so on. + * + * If you need a format in which all columns are of equal height + * and every bin contains a value, use CQInterpolated instead of + * ConstantQ. + */ + ComplexBlock process(const RealSequence &); + + /** + * Return the remaining constant-Q columns following the end of + * processing. Any buffered input is padded so as to ensure that + * all input provided to process() will have been returned. + */ + ComplexBlock getRemainingOutput(); + +private: + double m_sampleRate; + double m_maxFrequency; + double m_minFrequency; + int m_binsPerOctave; + int m_octaves; + + CQKernel *m_kernel; + CQKernel::Properties m_p; + int m_bigBlockSize; + + std::vector<Resampler *> m_decimators; + std::vector<RealSequence> m_buffers; + + int m_outputLatency; + + FFTReal *m_fft; + + void initialise(); + ComplexBlock processOctaveBlock(int octave); +}; + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/Makefile Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,6 @@ + +MAY_DIR := ../../may +MAY_YC := $(MAY_DIR)/bin/yc + +test: + $(MAY_YC) ./test.yeti
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/build.xml Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,74 @@ +<project name="cqt" default="test" basedir="."> + + <property name="maydir" value="${basedir}/../../may"/> + <property name="yetidir" value="${basedir}/../../yeti"/> + + <property name="jardir" value="${maydir}/ext/jar"/> + + <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/jzy3d-swt-0.9.1.jar:${jardir}/jzy3d-api-0.9.1.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/> + + <condition property="archtag" value="linux32"> + <os family="unix" arch="i386"/> + </condition> + <condition property="archtag" value="linux64"> + <os family="unix" arch="amd64"/> + </condition> + <condition property="archtag" value="win32"> + <os family="windows" arch="x86"/> + </condition> + <condition property="archtag" value="win64"> + <os family="windows" arch="amd64"/> + </condition> + <condition property="archtag" value="osx"> + <os family="mac"/> + </condition> + + <target name="taskdef"> + <taskdef name="yetic" classname="yeti.lang.compiler.YetiTask" + classpath="${yetidir}/yeti.jar:${maydir}/may.jar:${extjars}" /> + </target> + + <target name="prepare"> + <mkdir dir="${basedir}/classes"/> + </target> + + <target name="yeticlasses" depends="taskdef"> + <yetic srcdir="${basedir}" + destdir="${basedir}/classes" + includes="**/*.yeti" + preload="yeti/lang/std:yeti/lang/io"/> + </target> + + <target name="classes" depends="prepare,yeticlasses"/> + + <target name="jar" depends="classes,taskdef"> + <jar jarfile="${basedir}/cqt.jar"> + <fileset dir="${basedir}/classes" + includes="**/*.class" + excludes="**/test*.class"/> + </jar> + </target> + + <target name="testjar" depends="classes,taskdef"> + <jar jarfile="${basedir}/test.jar"> + <fileset dir="${basedir}/classes" + includes="**/test*.class"/> + </jar> + </target> + + <target name="test" depends="jar,testjar,taskdef"> + <java classpath="${basedir}/test.jar:${basedir}/cqt.jar:${maydir}/may.jar:${yetidir}/yeti.jar:${extjars}" + classname="test" + fork="true" failonerror="true"> + <sysproperty key="java.library.path" path="${maydir}/ext/native/${archtag}"/> + </java> + </target> + + <target name="clean"> + <delete dir="${basedir}/classes"/> + </target> + + <target name="rebuild" depends="clean,jar"/> + +</project> +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/cqt.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,260 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +module cqt; + +cqtkernel = load cqtkernel; +resample = load may.stream.resample; +manipulate = load may.stream.manipulate; +mat = load may.matrix; +cm = load may.matrix.complex; +framer = load may.stream.framer; +cplx = load may.complex; +fft = load may.transform.fft; +vec = load may.vector; +ch = load may.stream.channels; + +{ pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc; + +cqt { maxFreq, minFreq, binsPerOctave } str = + (sampleRate = str.sampleRate; + octaves = ceil (log2 (maxFreq / minFreq)); +// actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave)); + + kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave }; + +// eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave), fftSize = \(kdata.fftSize), hop = \(kdata.fftHop)"; + +// eprintln "atomsPerFrame = \(kdata.atomsPerFrame)"; + + padding = (kdata.fftSize * (pow 2 (octaves-1))); + +// eprintln "padding = \(padding)"; + + str = manipulate.paddedBy padding str; + + streams = manipulate.duplicated octaves str; + + // forward transform uses the conjugate-transposed kernel, inverse + // uses the original + kernel = cm.transposed (cm.conjugateTransposed kdata.kernel); + +// eprintln "have kernel"; + + fftFunc = fft.forward kdata.fftSize; + + cqblocks = + map do octave: + frames = map ch.mixedDown //!!! mono for now + (framer.frames kdata.fftSize [ Hop kdata.fftHop, Padded false ] + (resample.decimated (pow 2 octave) streams[octave])); + map do frame: + freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize)); +// eprintln "octave = \(octave), frame = \(vec.list frame)"; +// eprintln "octave = \(octave), freq = \(freq)"; + cm.product kernel (cm.newComplexColumnVector freq); + done frames; + done [0..octaves-1]; + + // The cqblocks list is a list<list<matrix>>. Each top-level list + // corresponds to an octave, from highest to lowest, each having + // twice as many elements in its list as the next octave. The + // sub-lists are sampled in time with an effective spacing of + // fftSize * 2^(octave-1) audio frames, and the matrices are row + // vectors with atomsPerFrame * binsPerOctave complex elements. + // + // *** + // + // In a typical constant-Q structure, each (2^(octaves-1) * + // fftHop) input frames gives us an output structure conceptually + // like this: + // + // [][][][][][][][] <- fftHop frames per highest-octave output value + // [][][][][][][][] layered as many times as binsPerOctave (here 2) + // [--][--][--][--] <- fftHop*2 frames for the next lower octave + // [--][--][--][--] etc + // [------][------] + // [------][------] + // [--------------] + // [--------------] + // + // *** + // + // But the kernel we're using here has more than one temporally + // spaced atom; each individual cell is a row vector with + // atomsPerFrame * binsPerOctave elements, but that actually + // represents a rectangular matrix of result cells with width + // atomsPerFrame and height binsPerOctave. The columns of this + // matrix (the atoms) then need to be spaced by 2^(octave-1) + // relative to those from the highest octave. + + // Reshape each row vector into the appropriate rectangular matrix + // and split into single-atom columns + + emptyHops = kdata.firstCentre / kdata.atomSpacing; //!!! int? round? +// maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops; +// eprintln "maxDrop = \(maxDrop)"; + + cqblocks = + map do octlist: + concatMap do rv: + cm.asColumns + (cm.generate do row col: + cm.at rv ((row * kdata.atomsPerFrame) + col) 0 + done { + rows = kdata.binsPerOctave, + columns = kdata.atomsPerFrame + }) + done octlist + done cqblocks; + + cqblocks = array (map2 do octlist octave: + d = emptyHops * (pow 2 (octaves-octave)) - emptyHops; +// eprintln "dropping \(d)"; + drop d octlist; + done cqblocks [1..octaves]); + + assembleBlock bits = + (//eprintln "assembleBlock: structure of bits is:"; + //eprintln (map length bits); + + rows = octaves * kdata.binsPerOctave; + columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame; + + cm.generate do row col: + + // bits structure: [1,2,4,8,...] + + // each elt of bits is a list of the chunks that should + // make up this block in that octave (lowest octave first) + + // each chunk has atomsPerFrame * binsPerOctave elts in it + + // row is disposed with 0 at the top, highest octave (in + // both pitch and index into bits structure) + + oct = int (row / binsPerOctave); + binNo = row % kdata.binsPerOctave; + + chunks = pow 2 oct; + colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame)); + atomNo = int (col / colsPerAtom); + atomOffset = col % colsPerAtom; + + if atomOffset == 0 and atomNo < length bits[oct] then + bits[oct][atomNo][binNo]; + else + cplx.zero + fi; + + done { rows, columns }; + ); + + assembleBlockSpectrogram bits = + (// As assembleBlock, but producing a dense magnitude + // spectrogram (rather than a complex output). (todo: + // interpolation, smoothing) + + //eprintln "assembleBlockSpectrogram: structure of bits is:"; + //eprintln (map length bits); + + rows = octaves * kdata.binsPerOctave; + columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame; + + mat.generate do row col: + + oct = int (row / binsPerOctave); + binNo = row % kdata.binsPerOctave; + + chunks = pow 2 oct; + colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame)); + atomNo = int (col / colsPerAtom); + + if atomNo < length bits[oct] then + cplx.magnitude bits[oct][atomNo][binNo]; + else + 0 + fi; + + done { rows, columns }; + ); + + processOctaveLists assembler octs = + case octs[0] of + block::rest: + (toAssemble = array + (map do oct: + n = kdata.atomsPerFrame * pow 2 oct; + if not empty? octs[oct] then + forBlock = array (take n octs[oct]); + octs[oct] := drop n octs[oct]; + forBlock + else + array [] + fi + done (keys octs)); + assembler toAssemble :. \(processOctaveLists assembler octs)); + _: [] + esac; + +//eprintln "cqblocks has \(length cqblocks) entries"; + + octaveLists = [:]; + + cqblocks = array cqblocks; + for [1..octaves] do oct: + octaveLists[octaves - oct] := cqblocks[oct-1]; + done; +/* + \() (map2 do octlist octave: +println "oct \(octaves) - \(octave) = \(octaves - octave)"; + octaveLists[octaves - octave] := octlist + done cqblocks [1..octaves]); +*/ +//eprintln "octaveLists keys are: \(keys octaveLists)"; + + { + kernel = kdata with { + binFrequencies = array + (concatMap do octave: + map do freq: + freq / (pow 2 octave); + done (reverse (list kdata.binFrequencies)) + done [0..octaves-1]) + }, + sampleRate, + octaves, + get cqComplex () = processOctaveLists assembleBlock octaveLists, + get cqSpectrogram () = processOctaveLists assembleBlockSpectrogram octaveLists, + } + ); + +{ cqt } +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/cqtkernel.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,158 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +module cqtkernel; + +vec = load may.vector; +complex = load may.complex; +window = load may.signal.window; +fft = load may.transform.fft; +cm = load may.matrix.complex; + +{ pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc; + +makeKernel { sampleRate, maxFreq, binsPerOctave } = + (q = 1; + atomHopFactor = 0.25; + thresh = 0.0005; + minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave)); + bigQ = q / ((pow 2 (1/binsPerOctave)) - 1); + + maxNK = round(bigQ * sampleRate / minFreq); + minNK = round(bigQ * sampleRate / + (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave)))); + + atomHop = round(minNK * atomHopFactor); + + firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop)); + + fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2)); + +// eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)"; +// eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)"; + + winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1; + + lastCentre = firstCentre + (winNr - 1) * atomHop; + + fftHop = (lastCentre + atomHop) - firstCentre; + +// eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)"; + + fftFunc = fft.forward fftSize; + + // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating + // function. We can't do that here; we need to generate real and imag + // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x). + + binFrequencies = array []; + + kernels = map do k: + + nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave)))); + + // the cq MATLAB toolbox uses a symmetric window for + // blackmanharris -- which is odd because it uses a periodic one + // for other types. Oh well + win = vec.divideBy nk + (vec.sqrt + (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk)); + + fk = minFreq * (pow 2 ((k-1)/binsPerOctave)); + + push binFrequencies fk; + + genKernel f = vec.multiply + [win, + vec.fromList + (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])]; + + reals = genKernel cos; + imags = genKernel sin; + + atomOffset = firstCentre - ceil(nk/2); + + map do i: + + shift = vec.zeros (atomOffset + ((i-1) * atomHop)); + + specKernel = fftFunc + (complex.complexArray + (vec.concat [shift, reals]) + (vec.concat [shift, imags])); + + map do c: + if complex.magnitude c <= thresh then complex.zero else c fi + done specKernel; + + done [1..winNr]; + + done [1..binsPerOctave]; + + kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels))); + +// eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))"; + + // Normalisation + + wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat)); + wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat)); + + subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1)); + square = cm.product (cm.conjugateTransposed subset) subset; + + diag = complex.magnitudes (cm.getDiagonal 0 square); + wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2); + + weight = (fftHop / fftSize) / (vec.mean (vec.abs wK)); + weight = sqrt(weight); + + kernel = cm.scaled weight kmat; + +// eprintln "weight = \(weight)"; + + { + kernel, + fftSize, + fftHop, + binsPerOctave, + atomsPerFrame = winNr, + atomSpacing = atomHop, + firstCentre, + maxFrequency = maxFreq, + minFrequency = minFreq, + binFrequencies, + bigQ + }); + +{ + makeKernel +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/experiment.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,50 @@ + +program experiment; + +af = load may.stream.audiofile; +plot = load may.plot; +cm = load may.matrix.complex; +mat = load may.matrix; +vec = load may.vector; +manipulate = load may.stream.manipulate; +syn = load may.stream.syntheticstream; + +{ cqt } = load cqt; +{ icqt } = load icqt; + +//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500); +//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4); +testStream = af.open "sweep-48000.wav"; +//testStream = af.open "sweep.wav"; + +eprintln "have test stream"; + +cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream; + +eprintln "bin frequencies: \(cq.kernel.binFrequencies)"; + +icq = icqt cq; + +eprintln "calculated icq..."; + +outfile = "test-output.wav"; +written = af.write icq outfile; + +eprintln "wrote \(written) to \(outfile)"; + +//for cq.cqSpectrogram do mm: +// for (mat.asColumns mm) (println . strJoin "," . vec.list); +//done; + +/* +bigM = mat.concatHorizontal (map cm.magnitudes cq.output); + +eprintln "overall output size = \(mat.size bigM)"; + +mat.print bigM; +*/ +//\() (plot.plot [Contour bigM]); +//\() (plot.plot [Grid bigM]); + +() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/icqt.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,145 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +module icqt; + +cqt = load cqt; +cm = load may.matrix.complex; +framer = load may.stream.framer; +win = load may.signal.window; +mm = load may.mathmisc; +vec = load may.vector; +resamp = load may.stream.resample; +manip = load may.stream.manipulate; +mat = load may.matrix; + +icqt cq = + (kdata = cq.kernel; + + // kdata.kernel is the kernel matrix for a single octave. It has + // width kdata.fftSize and height kdata.binsPerOctave * + // kdata.atomsPerFrame. + // + // kdata.fftHop is the overlap between kernel matrices in a single + // octave. + // + // cq.sampleRate is the output stream sample rate and cq.octaves + // the number of octaves. + // + // cq.cqComplex is the list of complex matrices containing the CQ + // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1) + // and height kdata.binsPerOctave * cq.octaves. + + bpo = kdata.binsPerOctave; + atomsPerFrame = kdata.atomsPerFrame; + octaves = cq.octaves; + + // transform a single block, all octaves tall, into an array + // (indexed by octave) of lists of individual columns (valid + // values for that octave only) + decomposeOctaves mat = array + (map do oct: + inverted = cm.fromRows (reverse (cm.asRows mat)); + octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1)); + gap = (mm.pow 2 oct) - 1; + pickFrom cols = + case cols of + c::cs: c::(pickFrom (drop gap cs)); + _: []; + esac; + pickFrom (cm.asColumns octMat); + done [0..octaves-1]); + + // transform a list of the arrays produced by decomposeOctaves + // into a single array (indexed by octave) of lists of the + // individual columns + flattenOctaves decomposed = + (flattenAux acc decomposed = + case decomposed of + chunk::rest: + flattenAux + (array + (map do oct: + acc[oct] ++ chunk[oct] + done [0..octaves-1])) + rest; + _: acc; + esac; + flattenAux (array (map \[] [0..octaves-1])) decomposed); + + // given a list of columns, deinterleave and pile up each sequence + // of atomsPerFrame columns into a single tall column and return + // the resulting list of tall columns + pile cols = + (pileAux acc cols = + if (length cols) >= atomsPerFrame then + atoms = take atomsPerFrame cols; + juggled = concat (reverse (cm.asRows (cm.fromColumns atoms))); + pileAux (juggled :: acc) (drop atomsPerFrame cols); + else + reverse acc + fi; + pileAux [] cols); + + octaveColumnLists = + map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex))); + + for octaveColumnLists do l: println "octave column list length: \(length l)" done; + + kernel = cm.transposed kdata.kernel; // right way around for the multiply + + spectra = + map do l: + map do col: + res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col)); + cm.columnSlice res 0 (kdata.fftSize / 2 + 1) + done l; + done octaveColumnLists; + + eprintln "calculated spectra, now to ifft, overlap-add..."; + + rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1]; + + resynthesised = + map2 do frames rate: + framer.complexStreamed rate kdata.fftSize + [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ] + frames; + done spectra rates; + + eprintln "have streams:"; + for resynthesised eprintln; + + resampled = map (resamp.resampledTo cq.sampleRate) resynthesised; + + manip.sum resampled; +); + +{icqt}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/plotfile.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,41 @@ + +program plotfile; + +af = load may.stream.audiofile; +plot = load may.plot; +cm = load may.matrix.complex; +mat = load may.matrix; + +{ cqt } = load cqt; + +minFreq = 27.5; +maxRateDivisor = 3; +binsPerOctave = 60; + +plotfile f = + (testStream = af.open f; + eprintln "Opened file stream..."; + start = System#currentTimeMillis(); + cq = cqt { + maxFreq = testStream.sampleRate/maxRateDivisor, + minFreq, + binsPerOctave + } testStream; + middle = System#currentTimeMillis(); + eprintln "Generated kernel \(cm.size cq.kernel.kernel) and primed transform (\(cq.octaves) octaves), took \(middle-start)ms, now calculating..."; + bigM = mat.concatHorizontal cq.cqSpectrogram; + finish = System#currentTimeMillis(); + eprintln "Done, that part took \(finish-middle)ms, all CQ stuff took \(finish-start)ms"; + eprintln "Plotting..."; + \() (plot.plot [Contour bigM])); + +usage () = + (eprintln "\nUsage: plotfile file.wav"; + eprintln "\n Loads audio from file.wav and plots a \(binsPerOctave)bpo Constant-Q spectrogram"; + eprintln " from \(minFreq)Hz up to 1/\(maxRateDivisor) of the file samplerate"); + +case (list _argv) of +file::[]: plotfile file; +_: usage (); +esac; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/test.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,96 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +program test; + +{ runTests } = load may.test; + +// We want to test: +// +// Kernel design -- check size (number of bins, number of atoms); +// check an example kernel against known data +// +// Time alignment -- feed a dirac train, check that peaks in all bins +// align +// +// Frequency discrimination -- feed a sinusoid, check peaks +// +// Latency compensation -- for dirac at 0, check peak can be found at +// 0 plus the declared latency +// +// Signal-noise ratio +// +// Specimen output for simple test case + +tests = [ +"kernel" : load test_cqtkernel, +"frequency" : load test_frequency, +]; + +bad = sum (mapHash do name testHash: runTests name testHash done tests); + +if (bad > 0) then + println "\n** \(bad) test(s) failed!"; + threadExit 1 +else + () +fi + + + + +/* +//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500); +//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4); +testStream = af.open "sweep-48000.wav"; +//testStream = af.open "sweep.wav"; + +// So the stream is [ 0, 1, 0, -1, 0, 1, 0, -1, ... ] : +//testStream = manipulate.withDuration 64 (syn.sinusoid 8 2); + +testStream = manipulate.withDuration 32 (syn.pulseTrain 8 0.001); + +eprintln "have test stream"; + +cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream; + +eprintln "bin frequencies: \(cq.kernel.binFrequencies)"; + +bigM = mat.concatHorizontal (map cm.magnitudes cq.output); + +eprintln "overall output size = \(mat.size bigM)"; + +mat.print bigM; + +//\() (plot.plot [Contour bigM]); +\() (plot.plot [Grid bigM]); +*/ +//() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/test_cqtkernel.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,63 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +module test_cqtkernel; + +cm = load may.matrix.complex; +mm = load may.mathmisc; + +{ compare, compareUsing } = load may.test; + +{ makeKernel } = load cqtkernel; + +eps = 1e-7; + +compareClose = compareUsing do a b: abs (a - b) < eps done; + +[ + +"minimal": \( + k = makeKernel { sampleRate = 16, maxFreq = 8, binsPerOctave = 4 }; + compare k.binsPerOctave 4 and + compare (cm.size k.kernel) { + rows = k.binsPerOctave * k.atomsPerFrame, + columns = k.fftSize + } and + compareClose k.maxFrequency 8 and + compareClose k.minFrequency (4 * (mm.pow 2 (1/4))) and + compare k.atomsPerFrame 5 and + compare k.fftSize 32 and + compare (length k.binFrequencies) k.binsPerOctave and + compareClose (head k.binFrequencies) k.minFrequency and + compareClose (head (reverse k.binFrequencies)) k.maxFrequency +), + +] is hash<string, () -> boolean> +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/yeti/test_frequency.yeti Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,132 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +module test_frequency; + +mat = load may.matrix; +vec = load may.vector; +win = load may.signal.window; +mm = load may.mathmisc; +cm = load may.matrix.complex; +syn = load may.stream.syntheticstream; +plot = load may.plot; + +{ compare } = load may.test; + +{ cqt } = load cqt; + +// Test with a single windowed sinusoid, repeating at various frequencies + +sinTestStream sampleRate duration signalFreq = // duration is in samples + (sin = syn.sinusoid sampleRate signalFreq; + chunk = mat.getRow 0 (sin.read duration); + syn.precalculatedMono sampleRate (win.windowed win.hann chunk)); + +// We want to make a CQ transform spanning more than one octave, but +// not going all the way to fs/2 so we can test it also with +// frequencies above and below its extents + +sampleRate = 100; + +// fs/2 = 50 so 10->40 gives us 2 octaves +cqmin = 10; +cqmax = 40; +bpo = 4; // fairly arbitrary + +testFreqs = map (* 5) [ 0..10 ]; +duration = sampleRate * 2; + +threshold = 0.08; + +streamBuilder = sinTestStream sampleRate duration; + +binForFreq f = + mm.round (bpo * mm.log2 (f / cqmin)) - 1; + +report message matrix = + (eprintln message; + eprintln "matrix is:"; + mat.eprint matrix); +// chart = plot.plot [Grid matrix]; +// sleep 100; +// chart#dispose()); + +tests = mapIntoHash + do f: "freq_\(f)" done + do f: \( + str = streamBuilder f; + cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str; + spec = cq.cqSpectrogram; + rightSize = all id + (map do s: + compare (mat.size s) { + rows = cq.kernel.binsPerOctave * cq.octaves, + columns = cq.kernel.atomsPerFrame * mm.pow 2 (cq.octaves - 1) + } + done spec); + m = mat.concatHorizontal spec; +// println "binFrequencies = \(cq.kernel.binFrequencies)"; +// println "binForFreq \(f) = \(binForFreq f)"; + var colno = 0; + success = all id + (rightSize :: map do c: + // The test passes for this column if: + // + // * the max bin is the expected one, or + // + // * the expected max is out of range entirely (but + // we need to test _something_ in this case -- + // what?), or + // + // * all bins are below a threshold, or + // + // * this is an odd column and the expected max is in + // the lower octave + // + // We should also check that all values in the lower + // octave are zero for odd columns. + // + expected = binForFreq f; + good = + (expected < 0 or expected >= vec.length c) or + ((colno % 2 == 1) and expected < (vec.length c / 2)) or + (vec.max c < threshold) or + (vec.maxindex c == binForFreq f); + if (not good) then + report " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)" m; + fi; + colno := colno + 1; + good; + done (mat.asColumns m)); + success; + ) done + testFreqs; + +tests is hash<string, () -> boolean>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/CQInverse.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,411 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "CQInverse.h" + +#include "dsp/Resampler.h" +#include "dsp/MathUtilities.h" +#include "dsp/FFT.h" + +#include <algorithm> +#include <iostream> +#include <stdexcept> + +using std::vector; +using std::cerr; +using std::endl; + +//#define DEBUG_CQ 1 + +CQInverse::CQInverse(double sampleRate, + double minFreq, + double maxFreq, + int binsPerOctave) : + m_sampleRate(sampleRate), + m_maxFrequency(maxFreq), + m_minFrequency(minFreq), + m_binsPerOctave(binsPerOctave), + m_fft(0) +{ + if (minFreq <= 0.0 || maxFreq <= 0.0) { + throw std::invalid_argument("Frequency extents must be positive"); + } + + initialise(); +} + +CQInverse::~CQInverse() +{ + delete m_fft; + for (int i = 0; i < (int)m_upsamplers.size(); ++i) { + delete m_upsamplers[i]; + } + delete m_kernel; +} + +double +CQInverse::getMinFrequency() const +{ + return m_p.minFrequency / pow(2.0, m_octaves - 1); +} + +double +CQInverse::getBinFrequency(int bin) const +{ + return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave())); +} + +void +CQInverse::initialise() +{ + m_octaves = int(ceil(log2(m_maxFrequency / m_minFrequency))); + m_kernel = new CQKernel(m_sampleRate, m_maxFrequency, m_binsPerOctave); + m_p = m_kernel->getProperties(); + + // Use exact powers of two for resampling rates. They don't have + // to be related to our actual samplerate: the resampler only + // cares about the ratio, but it only accepts integer source and + // target rates, and if we start from the actual samplerate we + // risk getting non-integer rates for lower octaves + + int sourceRate = pow(2, m_octaves); + vector<int> latencies; + + // top octave, no resampling + latencies.push_back(0); + m_upsamplers.push_back(0); + + for (int i = 1; i < m_octaves; ++i) { + + int factor = pow(2, i); + + Resampler *r = new Resampler + (sourceRate / factor, sourceRate, 50, 0.05); + +#ifdef DEBUG_CQ + cerr << "inverse: octave " << i << ": resample from " << sourceRate/factor << " to " << sourceRate << endl; +#endif + + // See ConstantQ.cpp for discussion on latency -- output + // latency here is at target rate which, this way around, is + // what we want + + latencies.push_back(r->getLatency()); + m_upsamplers.push_back(r); + } + + // additionally we will have fftHop latency at individual octave + // rate (before upsampling) for the overlap-add in each octave + for (int i = 0; i < m_octaves; ++i) { + latencies[i] += m_p.fftHop * pow(2, i); + } + + // Now reverse the drop adjustment made in ConstantQ to align the + // atom centres across different octaves (but this time at output + // sample rate) + + int emptyHops = m_p.firstCentre / m_p.atomSpacing; + + vector<int> pushes; + for (int i = 0; i < m_octaves; ++i) { + int factor = pow(2, i); + int pushHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops; + int push = ((pushHops * m_p.fftHop) * factor) / m_p.atomsPerFrame; + pushes.push_back(push); + } + + int maxLatLessPush = 0; + for (int i = 0; i < m_octaves; ++i) { + int latLessPush = latencies[i] - pushes[i]; + if (latLessPush > maxLatLessPush) maxLatLessPush = latLessPush; + } + + int totalLatency = maxLatLessPush + 10; + if (totalLatency < 0) totalLatency = 0; + + m_outputLatency = totalLatency + m_p.firstCentre * pow(2, m_octaves-1); + +#ifdef DEBUG_CQ + cerr << "totalLatency = " << totalLatency << ", m_outputLatency = " << m_outputLatency << endl; +#endif + + for (int i = 0; i < m_octaves; ++i) { + + // Calculate the difference between the total latency applied + // across all octaves, and the existing latency due to the + // upsampler for this octave. + + int latencyPadding = totalLatency - latencies[i] + pushes[i]; + +#ifdef DEBUG_CQ + cerr << "octave " << i << ": push " << pushes[i] << ", resampler latency inc overlap space " << latencies[i] << ", latencyPadding = " << latencyPadding << " (/factor = " << latencyPadding / pow(2, i) << ")" << endl; +#endif + + m_buffers.push_back(RealSequence(latencyPadding, 0.0)); + } + + for (int i = 0; i < m_octaves; ++i) { + // Fixed-size buffer for IFFT overlap-add + m_olaBufs.push_back(RealSequence(m_p.fftSize, 0.0)); + } + + m_fft = new FFTReal(m_p.fftSize); +} + +CQInverse::RealSequence +CQInverse::process(const ComplexBlock &block) +{ + // The input data is of the form produced by ConstantQ::process -- + // an unknown number N of columns of varying height. We assert + // that N is a multiple of atomsPerFrame * 2^(octaves-1), as must + // be the case for data that came directly from our ConstantQ + // implementation. + + int widthProvided = block.size(); + + if (widthProvided == 0) { + return drawFromBuffers(); + } + + int blockWidth = m_p.atomsPerFrame * int(pow(2, m_octaves - 1)); + + if (widthProvided % blockWidth != 0) { + cerr << "ERROR: CQInverse::process: Input block size (" + << widthProvided + << ") must be a multiple of processing block width " + << "(atoms-per-frame * 2^(octaves-1) = " + << m_p.atomsPerFrame << " * 2^(" << m_octaves << "-1) = " + << blockWidth << ")" << endl; + throw std::invalid_argument + ("Input block size must be a multiple of processing block width"); + } + + // Procedure: + // + // 1. Slice the list of columns into a set of lists of columns, + // one per octave, each of width N / (2^octave-1) and height + // binsPerOctave, containing the values present in that octave + // + // 2. Group each octave list by atomsPerFrame columns at a time, + // and stack these so as to achieve a list, for each octave, of + // taller columns of height binsPerOctave * atomsPerFrame + // + // 3. For each taller column, take the product with the inverse CQ + // kernel (which is the conjugate of the forward kernel) and + // perform an inverse FFT + // + // 4. Overlap-add each octave's resynthesised blocks (unwindowed) + // + // 5. Resample each octave's overlap-add stream to the original + // rate + // + // 6. Sum the resampled streams and return + + for (int i = 0; i < m_octaves; ++i) { + + // Step 1 + + ComplexBlock oct; + + for (int j = 0; j < widthProvided; ++j) { + int h = block[j].size(); + if (h < m_binsPerOctave * (i+1)) { + continue; + } + ComplexColumn col(block[j].begin() + m_binsPerOctave * i, + block[j].begin() + m_binsPerOctave * (i+1)); + oct.push_back(col); + } + + // Steps 2, 3, 4, 5 + processOctave(i, oct); + } + + // Step 6 + return drawFromBuffers(); +} + +CQInverse::RealSequence +CQInverse::drawFromBuffers() +{ + // 6. Sum the resampled streams and return + + int available = 0; + + for (int i = 0; i < m_octaves; ++i) { + if (i == 0 || int(m_buffers[i].size()) < available) { + available = m_buffers[i].size(); + } + } + + RealSequence result(available, 0); + + if (available == 0) { + return result; + } + + for (int i = 0; i < m_octaves; ++i) { + for (int j = 0; j < available; ++j) { + result[j] += m_buffers[i][j]; + } + m_buffers[i] = RealSequence(m_buffers[i].begin() + available, + m_buffers[i].end()); + } + + return result; +} + +CQInverse::RealSequence +CQInverse::getRemainingOutput() +{ + for (int j = 0; j < m_octaves; ++j) { + int factor = pow(2, j); + int latency = (j > 0 ? m_upsamplers[j]->getLatency() : 0) / factor; + for (int i = 0; i < (latency + m_p.fftSize) / m_p.fftHop; ++i) { + overlapAddAndResample(j, RealSequence(m_olaBufs[j].size(), 0)); + } + } + + return drawFromBuffers(); +} + +void +CQInverse::processOctave(int octave, const ComplexBlock &columns) +{ + // 2. Group each octave list by atomsPerFrame columns at a time, + // and stack these so as to achieve a list, for each octave, of + // taller columns of height binsPerOctave * atomsPerFrame + + int ncols = columns.size(); + + if (ncols % m_p.atomsPerFrame != 0) { + cerr << "ERROR: CQInverse::process: Number of columns (" + << ncols + << ") in octave " << octave + << " must be a multiple of atoms-per-frame (" + << m_p.atomsPerFrame << ")" << endl; + throw std::invalid_argument + ("Columns in octave must be a multiple of atoms per frame"); + } + + for (int i = 0; i < ncols; i += m_p.atomsPerFrame) { + + ComplexColumn tallcol; + for (int b = 0; b < m_binsPerOctave; ++b) { + for (int a = 0; a < m_p.atomsPerFrame; ++a) { + tallcol.push_back(columns[i + a][m_binsPerOctave - b - 1]); + } + } + + processOctaveColumn(octave, tallcol); + } +} + +void +CQInverse::processOctaveColumn(int octave, const ComplexColumn &column) +{ + // 3. For each taller column, take the product with the inverse CQ + // kernel (which is the conjugate of the forward kernel) and + // perform an inverse FFT + + if ((int)column.size() != m_p.atomsPerFrame * m_binsPerOctave) { + cerr << "ERROR: CQInverse::processOctaveColumn: Height of column (" + << column.size() << ") in octave " << octave + << " must be atoms-per-frame * bins-per-octave (" + << m_p.atomsPerFrame << " * " << m_binsPerOctave << " = " + << m_p.atomsPerFrame * m_binsPerOctave << ")" << endl; + throw std::invalid_argument + ("Column height must match atoms-per-frame * bins-per-octave"); + } + + ComplexSequence transformed = m_kernel->processInverse(column); + + int halfLen = m_p.fftSize/2 + 1; + + RealSequence ri(halfLen, 0); + RealSequence ii(halfLen, 0); + + for (int i = 0; i < halfLen; ++i) { + ri[i] = transformed[i].real(); + ii[i] = transformed[i].imag(); + } + + RealSequence timeDomain(m_p.fftSize, 0); + + m_fft->inverse(ri.data(), ii.data(), timeDomain.data()); + + overlapAddAndResample(octave, timeDomain); +} + +void +CQInverse::overlapAddAndResample(int octave, const RealSequence &seq) +{ + // 4. Overlap-add each octave's resynthesised blocks (unwindowed) + // + // and + // + // 5. Resample each octave's overlap-add stream to the original + // rate + + if (seq.size() != m_olaBufs[octave].size()) { + cerr << "ERROR: CQInverse::overlapAdd: input sequence length (" + << seq.size() << ") is expected to match OLA buffer size (" + << m_olaBufs[octave].size() << ")" << endl; + throw std::invalid_argument + ("Input sequence length should match OLA buffer size"); + } + + RealSequence toResample(m_olaBufs[octave].begin(), + m_olaBufs[octave].begin() + m_p.fftHop); + + RealSequence resampled = + octave > 0 ? + m_upsamplers[octave]->process(toResample.data(), toResample.size()) : + toResample; + + m_buffers[octave].insert(m_buffers[octave].end(), + resampled.begin(), + resampled.end()); + + m_olaBufs[octave] = RealSequence(m_olaBufs[octave].begin() + m_p.fftHop, + m_olaBufs[octave].end()); + + RealSequence pad(m_p.fftHop, 0); + + m_olaBufs[octave].insert(m_olaBufs[octave].end(), + pad.begin(), + pad.end()); + + for (int i = 0; i < m_p.fftSize; ++i) { + m_olaBufs[octave][i] += seq[i]; + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/CQKernel.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,337 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "CQKernel.h" + +#include "dsp/MathUtilities.h" +#include "dsp/FFT.h" +#include "dsp/Window.h" + +#include <cmath> +#include <cassert> +#include <vector> +#include <iostream> +#include <algorithm> + +using std::vector; +using std::complex; +using std::cerr; +using std::endl; + +typedef std::complex<double> C; + +CQKernel::CQKernel(double sampleRate, double maxFreq, int binsPerOctave) : + m_fft(0) +{ + m_p.sampleRate = sampleRate; + m_p.maxFrequency = maxFreq; + m_p.binsPerOctave = binsPerOctave; + generateKernel(); +} + +CQKernel::~CQKernel() +{ + delete m_fft; +} + +void +CQKernel::generateKernel() +{ + double q = 1; + double atomHopFactor = 0.25; + double thresh = 0.0005; + + double bpo = m_p.binsPerOctave; + + m_p.minFrequency = (m_p.maxFrequency / 2) * pow(2, 1.0/bpo); + m_p.Q = q / (pow(2, 1.0/bpo) - 1.0); + + double maxNK = round(m_p.Q * m_p.sampleRate / m_p.minFrequency); + double minNK = round + (m_p.Q * m_p.sampleRate / + (m_p.minFrequency * pow(2, (bpo - 1.0) / bpo))); + + if (minNK == 0 || maxNK == 0) { + // most likely pathological parameters of some sort + cerr << "WARNING: CQKernel::generateKernel: minNK or maxNK is zero (minNK == " << minNK << ", maxNK == " << maxNK << "), not generating a kernel" << endl; + m_p.atomSpacing = 0; + m_p.firstCentre = 0; + m_p.fftSize = 0; + m_p.atomsPerFrame = 0; + m_p.lastCentre = 0; + m_p.fftHop = 0; + return; + } + + m_p.atomSpacing = round(minNK * atomHopFactor); + m_p.firstCentre = m_p.atomSpacing * ceil(ceil(maxNK / 2.0) / m_p.atomSpacing); + m_p.fftSize = MathUtilities::nextPowerOfTwo + (m_p.firstCentre + ceil(maxNK / 2.0)); + + m_p.atomsPerFrame = floor + (1.0 + (m_p.fftSize - ceil(maxNK / 2.0) - m_p.firstCentre) / m_p.atomSpacing); + + cerr << "atomsPerFrame = " << m_p.atomsPerFrame << " (atomHopFactor = " << atomHopFactor << ", atomSpacing = " << m_p.atomSpacing << ", fftSize = " << m_p.fftSize << ", maxNK = " << maxNK << ", firstCentre = " << m_p.firstCentre << ")" << endl; + + m_p.lastCentre = m_p.firstCentre + (m_p.atomsPerFrame - 1) * m_p.atomSpacing; + + m_p.fftHop = (m_p.lastCentre + m_p.atomSpacing) - m_p.firstCentre; + + cerr << "fftHop = " << m_p.fftHop << endl; + + m_fft = new FFT(m_p.fftSize); + + for (int k = 1; k <= m_p.binsPerOctave; ++k) { + + int nk = round(m_p.Q * m_p.sampleRate / + (m_p.minFrequency * pow(2, ((k-1.0) / bpo)))); + + // The MATLAB version uses a symmetric window, but our windows + // are periodic. A symmetric window of size N is a periodic + // one of size N-1 with the first element stuck on the end + Window<double> w(BlackmanHarrisWindow, nk-1); + vector<double> win = w.getWindowData(); + win.push_back(win[0]); + + for (int i = 0; i < (int)win.size(); ++i) { + win[i] = sqrt(win[i]) / nk; + } + + double fk = m_p.minFrequency * pow(2, ((k-1.0) / bpo)); + + vector<double> reals, imags; + + for (int i = 0; i < nk; ++i) { + double arg = (2.0 * M_PI * fk * i) / m_p.sampleRate; + reals.push_back(win[i] * cos(arg)); + imags.push_back(win[i] * sin(arg)); + } + + int atomOffset = m_p.firstCentre - int(ceil(nk/2.0)); + + for (int i = 0; i < m_p.atomsPerFrame; ++i) { + + int shift = atomOffset + (i * m_p.atomSpacing); + + vector<double> rin(m_p.fftSize, 0.0); + vector<double> iin(m_p.fftSize, 0.0); + + for (int j = 0; j < nk; ++j) { + rin[j + shift] = reals[j]; + iin[j + shift] = imags[j]; + } + + vector<double> rout(m_p.fftSize, 0.0); + vector<double> iout(m_p.fftSize, 0.0); + + m_fft->process(false, + rin.data(), iin.data(), + rout.data(), iout.data()); + + // Keep this dense for the moment (until after + // normalisation calculations) + + vector<C> row; + + for (int j = 0; j < m_p.fftSize; ++j) { + if (sqrt(rout[j] * rout[j] + iout[j] * iout[j]) < thresh) { + row.push_back(C(0, 0)); + } else { + row.push_back(C(rout[j] / m_p.fftSize, + iout[j] / m_p.fftSize)); + } + } + + m_kernel.origin.push_back(0); + m_kernel.data.push_back(row); + } + } + + assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame); + + // print density as diagnostic + + int nnz = 0; + for (int i = 0; i < (int)m_kernel.data.size(); ++i) { + for (int j = 0; j < (int)m_kernel.data[i].size(); ++j) { + if (m_kernel.data[i][j] != C(0, 0)) { + ++nnz; + } + } + } + + cerr << "size = " << m_kernel.data.size() << "*" << m_kernel.data[0].size() << " (fft size = " << m_p.fftSize << ")" << endl; + + assert((int)m_kernel.data.size() == m_p.binsPerOctave * m_p.atomsPerFrame); + assert((int)m_kernel.data[0].size() == m_p.fftSize); + + cerr << "density = " << double(nnz) / double(m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize) << " (" << nnz << " of " << m_p.binsPerOctave * m_p.atomsPerFrame * m_p.fftSize << ")" << endl; + + finaliseKernel(); +} + +static bool ccomparator(C &c1, C &c2) +{ + return abs(c1) < abs(c2); +} + +static int maxidx(vector<C> &v) +{ + return std::max_element(v.begin(), v.end(), ccomparator) - v.begin(); +} + +void +CQKernel::finaliseKernel() +{ + // calculate weight for normalisation + + int wx1 = maxidx(m_kernel.data[0]); + int wx2 = maxidx(m_kernel.data[m_kernel.data.size()-1]); + + vector<vector<C> > subset(m_kernel.data.size()); + for (int j = wx1; j <= wx2; ++j) { + for (int i = 0; i < (int)m_kernel.data.size(); ++i) { + subset[i].push_back(m_kernel.data[i][j]); + } + } + + int nrows = subset.size(); + int ncols = subset[0].size(); + vector<vector<C> > square(ncols); // conjugate transpose of subset * subset + + for (int i = 0; i < nrows; ++i) { + assert((int)subset[i].size() == ncols); + } + + for (int j = 0; j < ncols; ++j) { + for (int i = 0; i < ncols; ++i) { + C v(0, 0); + for (int k = 0; k < nrows; ++k) { + v += subset[k][i] * conj(subset[k][j]); + } + square[i].push_back(v); + } + } + + vector<double> wK; + double q = 1; //!!! duplicated from constructor + for (int i = round(1.0/q); i < ncols - round(1.0/q) - 2; ++i) { + wK.push_back(abs(square[i][i])); + } + + double weight = double(m_p.fftHop) / m_p.fftSize; + weight /= MathUtilities::mean(wK.data(), wK.size()); + weight = sqrt(weight); + + cerr << "weight = " << weight << endl; + + // apply normalisation weight, make sparse, and store conjugate + // (we use the adjoint or conjugate transpose of the kernel matrix + // for the forward transform, the plain kernel for the inverse + // which we expect to be less common) + + KernelMatrix sk; + + for (int i = 0; i < (int)m_kernel.data.size(); ++i) { + + sk.origin.push_back(0); + sk.data.push_back(vector<C>()); + + int lastNZ = 0; + for (int j = (int)m_kernel.data[i].size()-1; j >= 0; --j) { + if (abs(m_kernel.data[i][j]) != 0.0) { + lastNZ = j; + break; + } + } + + bool haveNZ = false; + for (int j = 0; j <= lastNZ; ++j) { + if (haveNZ || abs(m_kernel.data[i][j]) != 0.0) { + if (!haveNZ) sk.origin[i] = j; + haveNZ = true; + sk.data[i].push_back(conj(m_kernel.data[i][j]) * weight); + } + } + } + + m_kernel = sk; +} + +vector<C> +CQKernel::processForward(const vector<C> &cv) +{ + // straightforward matrix multiply (taking into account m_kernel's + // slightly-sparse representation) + + if (m_kernel.data.empty()) return vector<C>(); + + int nrows = m_p.binsPerOctave * m_p.atomsPerFrame; + + vector<C> rv(nrows, C()); + + for (int i = 0; i < nrows; ++i) { + int len = m_kernel.data[i].size(); + for (int j = 0; j < len; ++j) { + rv[i] += cv[j + m_kernel.origin[i]] * m_kernel.data[i][j]; + } + } + + return rv; +} + +vector<C> +CQKernel::processInverse(const vector<C> &cv) +{ + // matrix multiply by conjugate transpose of m_kernel. This is + // actually the original kernel as calculated, we just stored the + // conjugate-transpose of the kernel because we expect to be doing + // more forward transforms than inverse ones. + + if (m_kernel.data.empty()) return vector<C>(); + + int ncols = m_p.binsPerOctave * m_p.atomsPerFrame; + int nrows = m_p.fftSize; + + vector<C> rv(nrows, C()); + + for (int j = 0; j < ncols; ++j) { + int i0 = m_kernel.origin[j]; + int i1 = i0 + m_kernel.data[j].size(); + for (int i = i0; i < i1; ++i) { + rv[i] += cv[j] * conj(m_kernel.data[j][i - i0]); + } + } + + return rv; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/CQKernel.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,82 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef CQ_KERNEL_H +#define CQ_KERNEL_H + +#include <vector> +#include <complex> + +class FFT; + +class CQKernel +{ +public: + CQKernel(double sampleRate, double maxFreq, int binsPerOctave); + ~CQKernel(); + + struct Properties { + double sampleRate; + double maxFrequency; + double minFrequency; + int binsPerOctave; + int fftSize; + int fftHop; + int atomsPerFrame; + int atomSpacing; + int firstCentre; + int lastCentre; + double Q; + }; + + Properties getProperties() const { return m_p; } + + std::vector<std::complex<double> > processForward + (const std::vector<std::complex<double> > &); + + std::vector<std::complex<double> > processInverse + (const std::vector<std::complex<double> > &); + +private: + Properties m_p; + FFT *m_fft; + + struct KernelMatrix { + std::vector<int> origin; + std::vector<std::vector<std::complex<double> > > data; + }; + KernelMatrix m_kernel; + + void generateKernel(); + void finaliseKernel(); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/CQSpectrogram.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,261 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "CQSpectrogram.h" + +#include <iostream> +#include <stdexcept> + +using std::cerr; +using std::endl; + +CQSpectrogram::CQSpectrogram(double sampleRate, + double minFreq, double maxFreq, + int binsPerOctave, + Interpolation interpolation) : + m_cq(sampleRate, minFreq, maxFreq, binsPerOctave), + m_interpolation(interpolation) +{ +} + +CQSpectrogram::~CQSpectrogram() +{ +} + +CQSpectrogram::RealBlock +CQSpectrogram::process(const RealSequence &td) +{ + return postProcess(m_cq.process(td), false); +} + +CQSpectrogram::RealBlock +CQSpectrogram::getRemainingOutput() +{ + return postProcess(m_cq.getRemainingOutput(), true); +} + +CQSpectrogram::RealBlock +CQSpectrogram::postProcess(const ComplexBlock &cq, bool insist) +{ + int width = cq.size(); + + // convert to magnitudes + RealBlock spec; + for (int i = 0; i < width; ++i) { + int height = cq[i].size(); + RealColumn col(height, 0); + for (int j = 0; j < height; ++j) { + col[j] = abs(cq[i][j]); + } + spec.push_back(col); + } + + if (m_interpolation == InterpolateZeros) { + for (int i = 0; i < width; ++i) { + int sh = spec[i].size(); + int fh = getTotalBins(); + for (int j = sh; j < fh; ++j) { + spec[i].push_back(0); + } + } + return spec; + } + + for (int i = 0; i < width; ++i) { + m_buffer.push_back(spec[i]); + } + + if (m_interpolation == InterpolateHold) { + return fetchHold(insist); + } else { + return fetchLinear(insist); + } +} + +CQSpectrogram::RealBlock +CQSpectrogram::fetchHold(bool) +{ + RealBlock out; + + int width = m_buffer.size(); + int height = getTotalBins(); + + for (int i = 0; i < width; ++i) { + + RealColumn col = m_buffer[i]; + + int thisHeight = col.size(); + int prevHeight = m_prevColumn.size(); + + for (int j = thisHeight; j < height; ++j) { + if (j < prevHeight) { + col.push_back(m_prevColumn[j]); + } else { + col.push_back(0.0); + } + } + + m_prevColumn = col; + out.push_back(col); + } + + m_buffer.clear(); + + return out; +} + +CQSpectrogram::RealBlock +CQSpectrogram::fetchLinear(bool insist) +{ + RealBlock out; + + //!!! This is surprisingly messy. I must be missing something. + + // We can only return any data when we have at least one column + // that has the full height in the buffer, that is not the first + // column. + // + // If the first col has full height, and there is another one + // later that also does, then we can interpolate between those, up + // to but not including the second full height column. Then we + // drop and return the columns we interpolated, leaving the second + // full-height col as the first col in the buffer. And repeat as + // long as enough columns are available. + // + // If the first col does not have full height, then (so long as + // we're following the logic above) we must simply have not yet + // reached the first full-height column in the CQ output, and we + // can interpolate nothing. + + int width = m_buffer.size(); + int height = getTotalBins(); + + if (width == 0) return out; + + int firstFullHeight = -1; + int secondFullHeight = -1; + + for (int i = 0; i < width; ++i) { + if ((int)m_buffer[i].size() == height) { + if (firstFullHeight == -1) { + firstFullHeight = i; + } else if (secondFullHeight == -1) { + secondFullHeight = i; + break; + } + } + } + +// cerr << "fetchLinear: firstFullHeight = " << firstFullHeight << ", secondFullHeight = " << secondFullHeight << endl; + + if (firstFullHeight < 0) { + if (insist) { + return fetchHold(true); + } else { + return out; + } + } else if (firstFullHeight > 0) { + // can interpolate nothing, stash up to first full height & recurse + out = RealBlock(m_buffer.begin(), m_buffer.begin() + firstFullHeight); + m_buffer = RealBlock(m_buffer.begin() + firstFullHeight, m_buffer.end()); + RealBlock more = fetchLinear(insist); + out.insert(out.end(), more.begin(), more.end()); + return out; + } else if (secondFullHeight < 0) { + // firstFullHeight == 0, but there is no second full height -- + // wait for it unless insist flag is set + if (insist) { + return fetchHold(true); + } else { + return out; + } + } else { + // firstFullHeight == 0 and secondFullHeight also valid. Can interpolate + out = linearInterpolated(m_buffer, 0, secondFullHeight); + m_buffer = RealBlock(m_buffer.begin() + secondFullHeight, m_buffer.end()); + RealBlock more = fetchLinear(insist); + out.insert(out.end(), more.begin(), more.end()); + return out; + } +} + +CQSpectrogram::RealBlock +CQSpectrogram::linearInterpolated(const RealBlock &g, int x0, int x1) +{ + // g must be a grid with full-height columns at x0 and x1 + + if (x0 >= x1) { + throw std::logic_error("x0 >= x1"); + } + if (x1 >= (int)g.size()) { + throw std::logic_error("x1 >= g.size()"); + } + if (g[x0].size() != g[x1].size()) { + throw std::logic_error("x0 and x1 are not the same height"); + } + + int height = g[x0].size(); + int width = x1 - x0; + + RealBlock out(g.begin() + x0, g.begin() + x1); + + for (int y = 0; y < height; ++y) { + + int spacing = width; + for (int i = 1; i < width; ++i) { + int thisHeight = g[x0 + i].size(); + if (thisHeight > height) { + throw std::logic_error("First column not full-height"); + } + if (thisHeight > y) { + spacing = i; + break; + } + } + + if (spacing < 2) continue; + + for (int i = 0; i + spacing <= width; i += spacing) { + for (int j = 1; j < spacing; ++j) { + double proportion = double(j)/double(spacing); + double interpolated = + g[x0 + i][y] * (1.0 - proportion) + + g[x0 + i + spacing][y] * proportion; + out[i + j].push_back(interpolated); + } + } + } + + return out; +} + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/ConstantQ.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,352 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "ConstantQ.h" + +#include "CQKernel.h" + +#include "dsp/Resampler.h" +#include "dsp/MathUtilities.h" +#include "dsp/FFT.h" + +#include <algorithm> +#include <iostream> +#include <stdexcept> + +using std::vector; +using std::cerr; +using std::endl; + +//#define DEBUG_CQ 1 + +ConstantQ::ConstantQ(double sampleRate, + double minFreq, + double maxFreq, + int binsPerOctave) : + m_sampleRate(sampleRate), + m_maxFrequency(maxFreq), + m_minFrequency(minFreq), + m_binsPerOctave(binsPerOctave), + m_fft(0) +{ + if (minFreq <= 0.0 || maxFreq <= 0.0) { + throw std::invalid_argument("Frequency extents must be positive"); + } + + initialise(); +} + +ConstantQ::~ConstantQ() +{ + delete m_fft; + for (int i = 0; i < (int)m_decimators.size(); ++i) { + delete m_decimators[i]; + } + delete m_kernel; +} + +double +ConstantQ::getMinFrequency() const +{ + return m_p.minFrequency / pow(2.0, m_octaves - 1); +} + +double +ConstantQ::getBinFrequency(int bin) const +{ + return getMinFrequency() * pow(2, (double(bin) / getBinsPerOctave())); +} + +void +ConstantQ::initialise() +{ + m_octaves = int(ceil(log2(m_maxFrequency / m_minFrequency))); + m_kernel = new CQKernel(m_sampleRate, m_maxFrequency, m_binsPerOctave); + m_p = m_kernel->getProperties(); + + // Use exact powers of two for resampling rates. They don't have + // to be related to our actual samplerate: the resampler only + // cares about the ratio, but it only accepts integer source and + // target rates, and if we start from the actual samplerate we + // risk getting non-integer rates for lower octaves + + int sourceRate = pow(2, m_octaves); + vector<int> latencies; + + // top octave, no resampling + latencies.push_back(0); + m_decimators.push_back(0); + + for (int i = 1; i < m_octaves; ++i) { + + int factor = pow(2, i); + + Resampler *r = new Resampler + (sourceRate, sourceRate / factor, 50, 0.05); + +#ifdef DEBUG_CQ + cerr << "forward: octave " << i << ": resample from " << sourceRate << " to " << sourceRate / factor << endl; +#endif + + // We need to adapt the latencies so as to get the first input + // sample to be aligned, in time, at the decimator output + // across all octaves. + // + // Our decimator uses a linear phase filter, but being causal + // it is not zero phase: it has a latency that depends on the + // decimation factor. Those latencies have been calculated + // per-octave and are available to us in the latencies + // array. Left to its own devices, the first input sample will + // appear at output sample 0 in the highest octave (where no + // decimation is needed), sample number latencies[1] in the + // next octave down, latencies[2] in the next one, etc. We get + // to apply some artificial per-octave latency after the + // decimator in the processing chain, in order to compensate + // for the differing latencies associated with different + // decimation factors. How much should we insert? + // + // The outputs of the decimators are at different rates (in + // terms of the relation between clock time and samples) and + // we want them aligned in terms of time. So, for example, a + // latency of 10 samples with a decimation factor of 2 is + // equivalent to a latency of 20 with no decimation -- they + // both result in the first output sample happening at the + // same equivalent time in milliseconds. + // + // So here we record the latency added by the decimator, in + // terms of the sample rate of the undecimated signal. Then we + // use that to compensate in a moment, when we've discovered + // what the longest latency across all octaves is. + + latencies.push_back(r->getLatency() * factor); + m_decimators.push_back(r); + } + + m_bigBlockSize = m_p.fftSize * pow(2, m_octaves - 1); + + // Now add in the extra padding and compensate for hops that must + // be dropped in order to align the atom centres across + // octaves. Again this is a bit trickier because we are doing it + // at input rather than output and so must work in per-octave + // sample rates rather than output blocks + + int emptyHops = m_p.firstCentre / m_p.atomSpacing; + + vector<int> drops; + for (int i = 0; i < m_octaves; ++i) { + int factor = pow(2, i); + int dropHops = emptyHops * pow(2, m_octaves - i - 1) - emptyHops; + int drop = ((dropHops * m_p.fftHop) * factor) / m_p.atomsPerFrame; + drops.push_back(drop); + } + + int maxLatPlusDrop = 0; + for (int i = 0; i < m_octaves; ++i) { + int latPlusDrop = latencies[i] + drops[i]; + if (latPlusDrop > maxLatPlusDrop) maxLatPlusDrop = latPlusDrop; + } + + int totalLatency = maxLatPlusDrop; + + int lat0 = totalLatency - latencies[0] - drops[0]; + totalLatency = ceil(double(lat0 / m_p.fftHop) * m_p.fftHop) + + latencies[0] + drops[0]; + + // We want (totalLatency - latencies[i]) to be a multiple of 2^i + // for each octave i, so that we do not end up with fractional + // octave latencies below. In theory this is hard, in practice if + // we ensure it for the last octave we should be OK. + double finalOctLat = latencies[m_octaves-1]; + double finalOctFact = pow(2, m_octaves-1); + totalLatency = + int(round(finalOctLat + + finalOctFact * + ceil((totalLatency - finalOctLat) / finalOctFact))); + +#ifdef DEBUG_CQ + cerr << "total latency = " << totalLatency << endl; +#endif + + // Padding as in the reference (will be introduced with the + // latency compensation in the loop below) + m_outputLatency = totalLatency + m_bigBlockSize + - m_p.firstCentre * pow(2, m_octaves-1); + +#ifdef DEBUG_CQ + cerr << "m_bigBlockSize = " << m_bigBlockSize << ", firstCentre = " + << m_p.firstCentre << ", m_octaves = " << m_octaves + << ", so m_outputLatency = " << m_outputLatency << endl; +#endif + + for (int i = 0; i < m_octaves; ++i) { + + double factor = pow(2, i); + + // Calculate the difference between the total latency applied + // across all octaves, and the existing latency due to the + // decimator for this octave, and then convert it back into + // the sample rate appropriate for the output latency of this + // decimator -- including one additional big block of padding + // (as in the reference). + + double octaveLatency = + double(totalLatency - latencies[i] - drops[i] + + m_bigBlockSize) / factor; + +#ifdef DEBUG_CQ + cerr << "octave " << i << ": resampler latency = " << latencies[i] + << ", drop " << drops[i] << " (/factor = " << drops[i]/factor + << "), octaveLatency = " << octaveLatency << " -> " + << int(round(octaveLatency)) << " (diff * factor = " + << (octaveLatency - round(octaveLatency)) << " * " + << factor << " = " + << (octaveLatency - round(octaveLatency)) * factor << ")" << endl; + + cerr << "double(" << totalLatency << " - " + << latencies[i] << " - " << drops[i] << " + " + << m_bigBlockSize << ") / " << factor << " = " + << octaveLatency << endl; +#endif + + m_buffers.push_back + (RealSequence(int(round(octaveLatency)), 0.0)); + } + + m_fft = new FFTReal(m_p.fftSize); +} + +ConstantQ::ComplexBlock +ConstantQ::process(const RealSequence &td) +{ + m_buffers[0].insert(m_buffers[0].end(), td.begin(), td.end()); + + for (int i = 1; i < m_octaves; ++i) { + RealSequence dec = m_decimators[i]->process(td.data(), td.size()); + m_buffers[i].insert(m_buffers[i].end(), dec.begin(), dec.end()); + } + + ComplexBlock out; + + while (true) { + + // We could have quite different remaining sample counts in + // different octaves, because (apart from the predictable + // added counts for decimator output on each block) we also + // have variable additional latency per octave + bool enough = true; + for (int i = 0; i < m_octaves; ++i) { + int required = m_p.fftSize * pow(2, m_octaves - i - 1); + if ((int)m_buffers[i].size() < required) { + enough = false; + } + } + if (!enough) break; + + int base = out.size(); + int totalColumns = pow(2, m_octaves - 1) * m_p.atomsPerFrame; + for (int i = 0; i < totalColumns; ++i) { + out.push_back(ComplexColumn()); + } + + for (int octave = 0; octave < m_octaves; ++octave) { + + int blocksThisOctave = pow(2, (m_octaves - octave - 1)); + + for (int b = 0; b < blocksThisOctave; ++b) { + ComplexBlock block = processOctaveBlock(octave); + + for (int j = 0; j < m_p.atomsPerFrame; ++j) { + + int target = base + + (b * (totalColumns / blocksThisOctave) + + (j * ((totalColumns / blocksThisOctave) / + m_p.atomsPerFrame))); + + while (int(out[target].size()) < + m_p.binsPerOctave * (octave + 1)) { + out[target].push_back(Complex()); + } + + for (int i = 0; i < m_p.binsPerOctave; ++i) { + out[target][m_p.binsPerOctave * octave + i] = + block[j][m_p.binsPerOctave - i - 1]; + } + } + } + } + } + + return out; +} + +ConstantQ::ComplexBlock +ConstantQ::getRemainingOutput() +{ + // Same as padding added at start, though rounded up + int pad = ceil(double(m_outputLatency) / m_bigBlockSize) * m_bigBlockSize; + RealSequence zeros(pad, 0.0); + return process(zeros); +} + +ConstantQ::ComplexBlock +ConstantQ::processOctaveBlock(int octave) +{ + RealSequence ro(m_p.fftSize, 0.0); + RealSequence io(m_p.fftSize, 0.0); + + m_fft->forward(m_buffers[octave].data(), ro.data(), io.data()); + + RealSequence shifted; + shifted.insert(shifted.end(), + m_buffers[octave].begin() + m_p.fftHop, + m_buffers[octave].end()); + m_buffers[octave] = shifted; + + ComplexSequence cv; + for (int i = 0; i < m_p.fftSize; ++i) { + cv.push_back(Complex(ro[i], io[i])); + } + + ComplexSequence cqrowvec = m_kernel->processForward(cv); + + // Reform into a column matrix + ComplexBlock cqblock; + for (int j = 0; j < m_p.atomsPerFrame; ++j) { + cqblock.push_back(ComplexColumn()); + for (int i = 0; i < m_p.binsPerOctave; ++i) { + cqblock[j].push_back(cqrowvec[i * m_p.atomsPerFrame + j]); + } + } + + return cqblock; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/FFT.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,225 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "FFT.h" + +#include "MathUtilities.h" + +#include "kiss_fft.h" +#include "kiss_fftr.h" + +#include <cmath> + +#include <iostream> + +#include <stdexcept> + +class FFT::D +{ +public: + D(int n) : m_n(n) { + m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL); + m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL); + m_kin = new kiss_fft_cpx[m_n]; + m_kout = new kiss_fft_cpx[m_n]; + } + + ~D() { + kiss_fft_free(m_planf); + kiss_fft_free(m_plani); + delete[] m_kin; + delete[] m_kout; + } + + void process(bool inverse, + const double *ri, + const double *ii, + double *ro, + double *io) { + + for (int i = 0; i < m_n; ++i) { + m_kin[i].r = ri[i]; + m_kin[i].i = (ii ? ii[i] : 0.0); + } + + if (!inverse) { + + kiss_fft(m_planf, m_kin, m_kout); + + for (int i = 0; i < m_n; ++i) { + ro[i] = m_kout[i].r; + io[i] = m_kout[i].i; + } + + } else { + + kiss_fft(m_plani, m_kin, m_kout); + + double scale = 1.0 / m_n; + + for (int i = 0; i < m_n; ++i) { + ro[i] = m_kout[i].r * scale; + io[i] = m_kout[i].i * scale; + } + } + } + +private: + int m_n; + kiss_fft_cfg m_planf; + kiss_fft_cfg m_plani; + kiss_fft_cpx *m_kin; + kiss_fft_cpx *m_kout; +}; + +FFT::FFT(int n) : + m_d(new D(n)) +{ +} + +FFT::~FFT() +{ + delete m_d; +} + +void +FFT::process(bool inverse, + const double *p_lpRealIn, const double *p_lpImagIn, + double *p_lpRealOut, double *p_lpImagOut) +{ + m_d->process(inverse, + p_lpRealIn, p_lpImagIn, + p_lpRealOut, p_lpImagOut); +} + +class FFTReal::D +{ +public: + D(int n) : m_n(n) { + if (n % 2) { + throw std::invalid_argument + ("nsamples must be even in FFTReal constructor"); + } + m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL); + m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL); + m_c = new kiss_fft_cpx[m_n]; + } + + ~D() { + kiss_fftr_free(m_planf); + kiss_fftr_free(m_plani); + delete[] m_c; + } + + void forward(const double *ri, double *ro, double *io) { + + kiss_fftr(m_planf, ri, m_c); + + for (int i = 0; i <= m_n/2; ++i) { + ro[i] = m_c[i].r; + io[i] = m_c[i].i; + } + + for (int i = 0; i + 1 < m_n/2; ++i) { + ro[m_n - i - 1] = ro[i + 1]; + io[m_n - i - 1] = -io[i + 1]; + } + } + + void forwardMagnitude(const double *ri, double *mo) { + + double *io = new double[m_n]; + + forward(ri, mo, io); + + for (int i = 0; i < m_n; ++i) { + mo[i] = sqrt(mo[i] * mo[i] + io[i] * io[i]); + } + + delete[] io; + } + + void inverse(const double *ri, const double *ii, double *ro) { + + // kiss_fftr.h says + // "input freqdata has nfft/2+1 complex points" + + for (int i = 0; i < m_n/2 + 1; ++i) { + m_c[i].r = ri[i]; + m_c[i].i = ii[i]; + } + + kiss_fftri(m_plani, m_c, ro); + + double scale = 1.0 / m_n; + + for (int i = 0; i < m_n; ++i) { + ro[i] *= scale; + } + } + +private: + int m_n; + kiss_fftr_cfg m_planf; + kiss_fftr_cfg m_plani; + kiss_fft_cpx *m_c; +}; + +FFTReal::FFTReal(int n) : + m_d(new D(n)) +{ +} + +FFTReal::~FFTReal() +{ + delete m_d; +} + +void +FFTReal::forward(const double *ri, double *ro, double *io) +{ + m_d->forward(ri, ro, io); +} + +void +FFTReal::forwardMagnitude(const double *ri, double *mo) +{ + m_d->forwardMagnitude(ri, mo); +} + +void +FFTReal::inverse(const double *ri, const double *ii, double *ro) +{ + m_d->inverse(ri, ii, ro); +} + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/FFT.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,128 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef FFT_H +#define FFT_H + +class FFT +{ +public: + /** + * Construct an FFT object to carry out complex-to-complex + * transforms of size nsamples. nsamples does not have to be a + * power of two. + */ + FFT(int nsamples); + ~FFT(); + + /** + * Carry out a forward or inverse transform (depending on the + * value of inverse) of size nsamples, where nsamples is the value + * provided to the constructor above. + * + * realIn and (where present) imagIn should contain nsamples each, + * and realOut and imagOut should point to enough space to receive + * nsamples each. + * + * imagIn may be NULL if the signal is real, but the other + * pointers must be valid. + * + * The inverse transform is scaled by 1/nsamples. + */ + void process(bool inverse, + const double *realIn, const double *imagIn, + double *realOut, double *imagOut); + +private: + class D; + D *m_d; +}; + +class FFTReal +{ +public: + /** + * Construct an FFT object to carry out real-to-complex transforms + * of size nsamples. nsamples does not have to be a power of two, + * but it does have to be even. (Use the complex-complex FFT above + * if you need an odd FFT size. This constructor will throw + * std::invalid_argument if nsamples is odd.) + */ + FFTReal(int nsamples); + ~FFTReal(); + + /** + * Carry out a forward real-to-complex transform of size nsamples, + * where nsamples is the value provided to the constructor above. + * + * realIn, realOut, and imagOut must point to (enough space for) + * nsamples values. For consistency with the FFT class above, and + * compatibility with existing code, the conjugate half of the + * output is returned even though it is redundant. + */ + void forward(const double *realIn, + double *realOut, double *imagOut); + + /** + * Carry out a forward real-to-complex transform of size nsamples, + * where nsamples is the value provided to the constructor + * above. Return only the magnitudes of the complex output values. + * + * realIn and magOut must point to (enough space for) nsamples + * values. For consistency with the FFT class above, and + * compatibility with existing code, the conjugate half of the + * output is returned even though it is redundant. + */ + void forwardMagnitude(const double *realIn, double *magOut); + + /** + * Carry out an inverse real transform (i.e. complex-to-real) of + * size nsamples, where nsamples is the value provided to the + * constructor above. + * + * realIn and imagIn should point to at least nsamples/2+1 values; + * if more are provided, only the first nsamples/2+1 values of + * each will be used (the conjugate half will always be deduced + * from the first nsamples/2+1 rather than being read from the + * input data). realOut should point to enough space to receive + * nsamples values. + * + * The inverse transform is scaled by 1/nsamples. + */ + void inverse(const double *realIn, const double *imagIn, + double *realOut); + +private: + class D; + D *m_d; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/KaiserWindow.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,83 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "KaiserWindow.h" + +#include "MathUtilities.h" + +KaiserWindow::Parameters +KaiserWindow::parametersForTransitionWidth(double attenuation, + double transition) +{ + Parameters p; + p.length = 1 + (attenuation > 21.0 ? + ceil((attenuation - 7.95) / (2.285 * transition)) : + ceil(5.79 / transition)); + p.beta = (attenuation > 50.0 ? + 0.1102 * (attenuation - 8.7) : + attenuation > 21.0 ? + 0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) : + 0); + return p; +} + +static double besselTerm(double x, int i) +{ + if (i == 0) { + return 1; + } else { + double f = MathUtilities::factorial(i); + return pow(x/2, i*2) / (f*f); + } +} + +static double bessel0(double x) +{ + double b = 0.0; + for (int i = 0; i < 20; ++i) { + b += besselTerm(x, i); + } + return b; +} + +void +KaiserWindow::init() +{ + double denominator = bessel0(m_beta); + bool even = (m_length % 2 == 0); + for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) { + double k = double(2*i) / double(m_length-1) - 1.0; + m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator); + } + for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) { + m_window.push_back(m_window[int(m_length/2) - i - 1]); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/KaiserWindow.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,123 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef KAISER_WINDOW_H +#define KAISER_WINDOW_H + +#include <vector> +#include <cmath> + +/** + * Kaiser window: A windower whose bandwidth and sidelobe height + * (signal-noise ratio) can be specified. These parameters are traded + * off against the window length. + */ +class KaiserWindow +{ +public: + struct Parameters { + int length; + double beta; + }; + + /** + * Construct a Kaiser windower with the given length and beta + * parameter. + */ + KaiserWindow(Parameters p) : m_length(p.length), m_beta(p.beta) { init(); } + + /** + * Construct a Kaiser windower with the given attenuation in dB + * and transition width in samples. + */ + static KaiserWindow byTransitionWidth(double attenuation, + double transition) { + return KaiserWindow + (parametersForTransitionWidth(attenuation, transition)); + } + + /** + * Construct a Kaiser windower with the given attenuation in dB + * and transition bandwidth in Hz for the given samplerate. + */ + static KaiserWindow byBandwidth(double attenuation, + double bandwidth, + double samplerate) { + return KaiserWindow + (parametersForBandwidth(attenuation, bandwidth, samplerate)); + } + + /** + * Obtain the parameters necessary for a Kaiser window of the + * given attenuation in dB and transition width in samples. + */ + static Parameters parametersForTransitionWidth(double attenuation, + double transition); + + /** + * Obtain the parameters necessary for a Kaiser window of the + * given attenuation in dB and transition bandwidth in Hz for the + * given samplerate. + */ + static Parameters parametersForBandwidth(double attenuation, + double bandwidth, + double samplerate) { + return parametersForTransitionWidth + (attenuation, (bandwidth * 2 * M_PI) / samplerate); + } + + int getLength() const { + return m_length; + } + + const double *getWindow() const { + return m_window.data(); + } + + void cut(double *src) const { + cut(src, src); + } + + void cut(const double *src, double *dst) const { + for (int i = 0; i < m_length; ++i) { + dst[i] = src[i] * m_window[i]; + } + } + +private: + int m_length; + double m_beta; + std::vector<double> m_window; + + void init(); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/MathUtilities.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,419 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "MathUtilities.h" + +#include <iostream> +#include <algorithm> +#include <vector> +#include <cmath> + + +double MathUtilities::mod(double x, double y) +{ + double a = floor( x / y ); + + double b = x - ( y * a ); + return b; +} + +double MathUtilities::princarg(double ang) +{ + double ValOut; + + ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI; + + return ValOut; +} + +void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm) +{ + unsigned int i; + double temp = 0.0; + double a=0.0; + + for( i = 0; i < len; i++) + { + temp = data[ i ]; + + a += ::pow( fabs(temp), double(alpha) ); + } + a /= ( double )len; + a = ::pow( a, ( 1.0 / (double) alpha ) ); + + *ANorm = a; +} + +double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha ) +{ + unsigned int i; + unsigned int len = data.size(); + double temp = 0.0; + double a=0.0; + + for( i = 0; i < len; i++) + { + temp = data[ i ]; + + a += ::pow( fabs(temp), double(alpha) ); + } + a /= ( double )len; + a = ::pow( a, ( 1.0 / (double) alpha ) ); + + return a; +} + +double MathUtilities::round(double x) +{ + if (x < 0) { + return -floor(-x + 0.5); + } else { + return floor(x + 0.5); + } +} + +double MathUtilities::median(const double *src, unsigned int len) +{ + if (len == 0) return 0; + + std::vector<double> scratch; + for (int i = 0; i < len; ++i) scratch.push_back(src[i]); + std::sort(scratch.begin(), scratch.end()); + + int middle = len/2; + if (len % 2 == 0) { + return (scratch[middle] + scratch[middle - 1]) / 2; + } else { + return scratch[middle]; + } +} + +double MathUtilities::sum(const double *src, unsigned int len) +{ + unsigned int i ; + double retVal =0.0; + + for( i = 0; i < len; i++) + { + retVal += src[ i ]; + } + + return retVal; +} + +double MathUtilities::mean(const double *src, unsigned int len) +{ + double retVal =0.0; + + if (len == 0) return 0; + + double s = sum( src, len ); + + retVal = s / (double)len; + + return retVal; +} + +double MathUtilities::mean(const std::vector<double> &src, + unsigned int start, + unsigned int count) +{ + double sum = 0.; + + if (count == 0) return 0; + + for (int i = 0; i < (int)count; ++i) + { + sum += src[start + i]; + } + + return sum / count; +} + +void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max) +{ + unsigned int i; + double temp = 0.0; + + if (len == 0) { + *min = *max = 0; + return; + } + + *min = data[0]; + *max = data[0]; + + for( i = 0; i < len; i++) + { + temp = data[ i ]; + + if( temp < *min ) + { + *min = temp ; + } + if( temp > *max ) + { + *max = temp ; + } + + } +} + +int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) +{ + unsigned int index = 0; + unsigned int i; + double temp = 0.0; + + double max = pData[0]; + + for( i = 0; i < Length; i++) + { + temp = pData[ i ]; + + if( temp > max ) + { + max = temp ; + index = i; + } + + } + + if (pMax) *pMax = max; + + + return index; +} + +int MathUtilities::getMax( const std::vector<double> & data, double* pMax ) +{ + unsigned int index = 0; + unsigned int i; + double temp = 0.0; + + double max = data[0]; + + for( i = 0; i < data.size(); i++) + { + temp = data[ i ]; + + if( temp > max ) + { + max = temp ; + index = i; + } + + } + + if (pMax) *pMax = max; + + + return index; +} + +void MathUtilities::circShift( double* pData, int length, int shift) +{ + shift = shift % length; + double temp; + int i,n; + + for( i = 0; i < shift; i++) + { + temp=*(pData + length - 1); + + for( n = length-2; n >= 0; n--) + { + *(pData+n+1)=*(pData+n); + } + + *pData = temp; + } +} + +int MathUtilities::compareInt (const void * a, const void * b) +{ + return ( *(int*)a - *(int*)b ); +} + +void MathUtilities::normalise(double *data, int length, NormaliseType type) +{ + switch (type) { + + case NormaliseNone: return; + + case NormaliseUnitSum: + { + double sum = 0.0; + for (int i = 0; i < length; ++i) { + sum += data[i]; + } + if (sum != 0.0) { + for (int i = 0; i < length; ++i) { + data[i] /= sum; + } + } + } + break; + + case NormaliseUnitMax: + { + double max = 0.0; + for (int i = 0; i < length; ++i) { + if (fabs(data[i]) > max) { + max = fabs(data[i]); + } + } + if (max != 0.0) { + for (int i = 0; i < length; ++i) { + data[i] /= max; + } + } + } + break; + + } +} + +void MathUtilities::normalise(std::vector<double> &data, NormaliseType type) +{ + switch (type) { + + case NormaliseNone: return; + + case NormaliseUnitSum: + { + double sum = 0.0; + for (int i = 0; i < (int)data.size(); ++i) sum += data[i]; + if (sum != 0.0) { + for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum; + } + } + break; + + case NormaliseUnitMax: + { + double max = 0.0; + for (int i = 0; i < (int)data.size(); ++i) { + if (fabs(data[i]) > max) max = fabs(data[i]); + } + if (max != 0.0) { + for (int i = 0; i < (int)data.size(); ++i) data[i] /= max; + } + } + break; + + } +} + +void MathUtilities::adaptiveThreshold(std::vector<double> &data) +{ + int sz = int(data.size()); + if (sz == 0) return; + + std::vector<double> smoothed(sz); + + int p_pre = 8; + int p_post = 7; + + for (int i = 0; i < sz; ++i) { + + int first = std::max(0, i - p_pre); + int last = std::min(sz - 1, i + p_post); + + smoothed[i] = mean(data, first, last - first + 1); + } + + for (int i = 0; i < sz; i++) { + data[i] -= smoothed[i]; + if (data[i] < 0.0) data[i] = 0.0; + } +} + +bool +MathUtilities::isPowerOfTwo(int x) +{ + if (x < 1) return false; + if (x & (x-1)) return false; + return true; +} + +int +MathUtilities::nextPowerOfTwo(int x) +{ + if (isPowerOfTwo(x)) return x; + if (x < 1) return 1; + int n = 1; + while (x) { x >>= 1; n <<= 1; } + return n; +} + +int +MathUtilities::previousPowerOfTwo(int x) +{ + if (isPowerOfTwo(x)) return x; + if (x < 1) return 1; + int n = 1; + x >>= 1; + while (x) { x >>= 1; n <<= 1; } + return n; +} + +int +MathUtilities::nearestPowerOfTwo(int x) +{ + if (isPowerOfTwo(x)) return x; + int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x); + if (x - n0 < n1 - x) return n0; + else return n1; +} + +double +MathUtilities::factorial(int x) +{ + if (x < 0) return 0; + double f = 1; + for (int i = 1; i <= x; ++i) { + f = f * i; + } + return f; +} + +int +MathUtilities::gcd(int a, int b) +{ + int c = a % b; + if (c == 0) { + return b; + } else { + return gcd(b, c); + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/MathUtilities.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,153 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef MATHUTILITIES_H +#define MATHUTILITIES_H + +#include <vector> + +#include "nan-inf.h" + +/** + * Static helper functions for simple mathematical calculations. + */ +class MathUtilities +{ +public: + /** + * Round x to the nearest integer. + */ + static double round( double x ); + + /** + * Return through min and max pointers the highest and lowest + * values in the given array of the given length. + */ + static void getFrameMinMax( const double* data, unsigned int len, double* min, double* max ); + + /** + * Return the mean of the given array of the given length. + */ + static double mean( const double* src, unsigned int len ); + + /** + * Return the mean of the subset of the given vector identified by + * start and count. + */ + static double mean( const std::vector<double> &data, + unsigned int start, unsigned int count ); + + /** + * Return the sum of the values in the given array of the given + * length. + */ + static double sum( const double* src, unsigned int len ); + + /** + * Return the median of the values in the given array of the given + * length. If the array is even in length, the returned value will + * be half-way between the two values adjacent to median. + */ + static double median( const double* src, unsigned int len ); + + /** + * The principle argument function. Map the phase angle ang into + * the range [-pi,pi). + */ + static double princarg( double ang ); + + /** + * Floating-point division modulus: return x % y. + */ + static double mod( double x, double y); + + static void getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm); + static double getAlphaNorm(const std::vector <double> &data, unsigned int alpha ); + + static void circShift( double* data, int length, int shift); + + static int getMax( double* data, unsigned int length, double* max = 0 ); + static int getMax( const std::vector<double> &data, double* max = 0 ); + static int compareInt(const void * a, const void * b); + + enum NormaliseType { + NormaliseNone, + NormaliseUnitSum, + NormaliseUnitMax + }; + + static void normalise(double *data, int length, + NormaliseType n = NormaliseUnitMax); + + static void normalise(std::vector<double> &data, + NormaliseType n = NormaliseUnitMax); + + /** + * Threshold the input/output vector data against a moving-mean + * average filter. + */ + static void adaptiveThreshold(std::vector<double> &data); + + /** + * Return true if x is 2^n for some integer n >= 0. + */ + static bool isPowerOfTwo(int x); + + /** + * Return the next higher integer power of two from x, e.g. 1300 + * -> 2048, 2048 -> 2048. + */ + static int nextPowerOfTwo(int x); + + /** + * Return the next lower integer power of two from x, e.g. 1300 -> + * 1024, 2048 -> 2048. + */ + static int previousPowerOfTwo(int x); + + /** + * Return the nearest integer power of two to x, e.g. 1300 -> 1024, + * 12 -> 16 (not 8; if two are equidistant, the higher is returned). + */ + static int nearestPowerOfTwo(int x); + + /** + * Return x! + */ + static double factorial(int x); // returns double in case it is large + + /** + * Return the greatest common divisor of natural numbers a and b. + */ + static int gcd(int a, int b); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/Resampler.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,436 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "Resampler.h" + +#include "MathUtilities.h" +#include "KaiserWindow.h" +#include "SincWindow.h" + +#include <iostream> +#include <vector> +#include <map> +#include <cassert> + +#include <pthread.h> + +using std::vector; +using std::map; +using std::cerr; +using std::endl; + +//#define DEBUG_RESAMPLER 1 +//#define DEBUG_RESAMPLER_VERBOSE 1 + +Resampler::Resampler(int sourceRate, int targetRate) : + m_sourceRate(sourceRate), + m_targetRate(targetRate) +{ + initialise(100, 0.02); +} + +Resampler::Resampler(int sourceRate, int targetRate, + double snr, double bandwidth) : + m_sourceRate(sourceRate), + m_targetRate(targetRate) +{ + initialise(snr, bandwidth); +} + +Resampler::~Resampler() +{ + delete[] m_phaseData; +} + +// peakToPole -> length -> beta -> window +static map<double, map<int, map<double, vector<double> > > > +knownFilters; + +static pthread_mutex_t +knownFilterMutex = PTHREAD_MUTEX_INITIALIZER; + +void +Resampler::initialise(double snr, double bandwidth) +{ + int higher = std::max(m_sourceRate, m_targetRate); + int lower = std::min(m_sourceRate, m_targetRate); + + m_gcd = MathUtilities::gcd(lower, higher); + m_peakToPole = higher / m_gcd; + + if (m_targetRate < m_sourceRate) { + // antialiasing filter, should be slightly below nyquist + m_peakToPole = m_peakToPole / (1.0 - bandwidth/2.0); + } + + KaiserWindow::Parameters params = + KaiserWindow::parametersForBandwidth(snr, bandwidth, higher / m_gcd); + + params.length = + (params.length % 2 == 0 ? params.length + 1 : params.length); + + params.length = + (params.length > 200001 ? 200001 : params.length); + + m_filterLength = params.length; + + vector<double> filter; + + pthread_mutex_lock(&knownFilterMutex); + + if (knownFilters[m_peakToPole][m_filterLength].find(params.beta) == + knownFilters[m_peakToPole][m_filterLength].end()) { + + KaiserWindow kw(params); + SincWindow sw(m_filterLength, m_peakToPole * 2); + + filter = vector<double>(m_filterLength, 0.0); + for (int i = 0; i < m_filterLength; ++i) filter[i] = 1.0; + sw.cut(filter.data()); + kw.cut(filter.data()); + + knownFilters[m_peakToPole][m_filterLength][params.beta] = filter; + } + + filter = knownFilters[m_peakToPole][m_filterLength][params.beta]; + + pthread_mutex_unlock(&knownFilterMutex); + + int inputSpacing = m_targetRate / m_gcd; + int outputSpacing = m_sourceRate / m_gcd; + +#ifdef DEBUG_RESAMPLER + cerr << "resample " << m_sourceRate << " -> " << m_targetRate + << ": inputSpacing " << inputSpacing << ", outputSpacing " + << outputSpacing << ": filter length " << m_filterLength + << endl; +#endif + + // Now we have a filter of (odd) length flen in which the lower + // sample rate corresponds to every n'th point and the higher rate + // to every m'th where n and m are higher and lower rates divided + // by their gcd respectively. So if x coordinates are on the same + // scale as our filter resolution, then source sample i is at i * + // (targetRate / gcd) and target sample j is at j * (sourceRate / + // gcd). + + // To reconstruct a single target sample, we want a buffer (real + // or virtual) of flen values formed of source samples spaced at + // intervals of (targetRate / gcd), in our example case 3. This + // is initially formed with the first sample at the filter peak. + // + // 0 0 0 0 a 0 0 b 0 + // + // and of course we have our filter + // + // f1 f2 f3 f4 f5 f6 f7 f8 f9 + // + // We take the sum of products of non-zero values from this buffer + // with corresponding values in the filter + // + // a * f5 + b * f8 + // + // Then we drop (sourceRate / gcd) values, in our example case 4, + // from the start of the buffer and fill until it has flen values + // again + // + // a 0 0 b 0 0 c 0 0 + // + // repeat to reconstruct the next target sample + // + // a * f1 + b * f4 + c * f7 + // + // and so on. + // + // Above I said the buffer could be "real or virtual" -- ours is + // virtual. We don't actually store all the zero spacing values, + // except for padding at the start; normally we store only the + // values that actually came from the source stream, along with a + // phase value that tells us how many virtual zeroes there are at + // the start of the virtual buffer. So the two examples above are + // + // 0 a b [ with phase 1 ] + // a b c [ with phase 0 ] + // + // Having thus broken down the buffer so that only the elements we + // need to multiply are present, we can also unzip the filter into + // every-nth-element subsets at each phase, allowing us to do the + // filter multiplication as a simply vector multiply. That is, rather + // than store + // + // f1 f2 f3 f4 f5 f6 f7 f8 f9 + // + // we store separately + // + // f1 f4 f7 + // f2 f5 f8 + // f3 f6 f9 + // + // Each time we complete a multiply-and-sum, we need to work out + // how many (real) samples to drop from the start of our buffer, + // and how many to add at the end of it for the next multiply. We + // know we want to drop enough real samples to move along by one + // computed output sample, which is our outputSpacing number of + // virtual buffer samples. Depending on the relationship between + // input and output spacings, this may mean dropping several real + // samples, one real sample, or none at all (and simply moving to + // a different "phase"). + + m_phaseData = new Phase[inputSpacing]; + + for (int phase = 0; phase < inputSpacing; ++phase) { + + Phase p; + + p.nextPhase = phase - outputSpacing; + while (p.nextPhase < 0) p.nextPhase += inputSpacing; + p.nextPhase %= inputSpacing; + + p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase)) + / inputSpacing)); + + int filtZipLength = int(ceil(double(m_filterLength - phase) + / inputSpacing)); + + for (int i = 0; i < filtZipLength; ++i) { + p.filter.push_back(filter[i * inputSpacing + phase]); + } + + m_phaseData[phase] = p; + } + +#ifdef DEBUG_RESAMPLER + int cp = 0; + int totDrop = 0; + for (int i = 0; i < inputSpacing; ++i) { + cerr << "phase = " << cp << ", drop = " << m_phaseData[cp].drop + << ", filter length = " << m_phaseData[cp].filter.size() + << ", next phase = " << m_phaseData[cp].nextPhase << endl; + totDrop += m_phaseData[cp].drop; + cp = m_phaseData[cp].nextPhase; + } + cerr << "total drop = " << totDrop << endl; +#endif + + // The May implementation of this uses a pull model -- we ask the + // resampler for a certain number of output samples, and it asks + // its source stream for as many as it needs to calculate + // those. This means (among other things) that the source stream + // can be asked for enough samples up-front to fill the buffer + // before the first output sample is generated. + // + // In this implementation we're using a push model in which a + // certain number of source samples is provided and we're asked + // for as many output samples as that makes available. But we + // can't return any samples from the beginning until half the + // filter length has been provided as input. This means we must + // either return a very variable number of samples (none at all + // until the filter fills, then half the filter length at once) or + // else have a lengthy declared latency on the output. We do the + // latter. (What do other implementations do?) + // + // We want to make sure the first "real" sample will eventually be + // aligned with the centre sample in the filter (it's tidier, and + // easier to do diagnostic calculations that way). So we need to + // pick the initial phase and buffer fill accordingly. + // + // Example: if the inputSpacing is 2, outputSpacing is 3, and + // filter length is 7, + // + // x x x x a b c ... input samples + // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ... + // i j k l ... output samples + // [--------|--------] <- filter with centre mark + // + // Let h be the index of the centre mark, here 3 (generally + // int(filterLength/2) for odd-length filters). + // + // The smallest n such that h + n * outputSpacing > filterLength + // is 2 (that is, ceil((filterLength - h) / outputSpacing)), and + // (h + 2 * outputSpacing) % inputSpacing == 1, so the initial + // phase is 1. + // + // To achieve our n, we need to pre-fill the "virtual" buffer with + // 4 zero samples: the x's above. This is int((h + n * + // outputSpacing) / inputSpacing). It's the phase that makes this + // buffer get dealt with in such a way as to give us an effective + // index for sample a of 9 rather than 8 or 10 or whatever. + // + // This gives us output latency of 2 (== n), i.e. output samples i + // and j will appear before the one in which input sample a is at + // the centre of the filter. + + int h = int(m_filterLength / 2); + int n = ceil(double(m_filterLength - h) / outputSpacing); + + m_phase = (h + n * outputSpacing) % inputSpacing; + + int fill = (h + n * outputSpacing) / inputSpacing; + + m_latency = n; + + m_buffer = vector<double>(fill, 0); + m_bufferOrigin = 0; + +#ifdef DEBUG_RESAMPLER + cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")" + << ", latency " << m_latency << endl; +#endif +} + +double +Resampler::reconstructOne() +{ + Phase &pd = m_phaseData[m_phase]; + double v = 0.0; + int n = pd.filter.size(); + + assert(n + m_bufferOrigin <= (int)m_buffer.size()); + + const double *const __restrict__ buf = m_buffer.data() + m_bufferOrigin; + const double *const __restrict__ filt = pd.filter.data(); + + for (int i = 0; i < n; ++i) { + // NB gcc can only vectorize this with -ffast-math + v += buf[i] * filt[i]; + } + + m_bufferOrigin += pd.drop; + m_phase = pd.nextPhase; + return v; +} + +int +Resampler::process(const double *src, double *dst, int n) +{ + for (int i = 0; i < n; ++i) { + m_buffer.push_back(src[i]); + } + + int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate)); + int outidx = 0; + +#ifdef DEBUG_RESAMPLER + cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << endl; +#endif + + double scaleFactor = (double(m_targetRate) / m_gcd) / m_peakToPole; + + while (outidx < maxout && + m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) { + dst[outidx] = scaleFactor * reconstructOne(); + outidx++; + } + + m_buffer = vector<double>(m_buffer.begin() + m_bufferOrigin, m_buffer.end()); + m_bufferOrigin = 0; + + return outidx; +} + +vector<double> +Resampler::process(const double *src, int n) +{ + int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate)); + vector<double> out(maxout, 0.0); + int got = process(src, out.data(), n); + assert(got <= maxout); + if (got < maxout) out.resize(got); + return out; +} + +vector<double> +Resampler::resample(int sourceRate, int targetRate, const double *data, int n) +{ + Resampler r(sourceRate, targetRate); + + int latency = r.getLatency(); + + // latency is the output latency. We need to provide enough + // padding input samples at the end of input to guarantee at + // *least* the latency's worth of output samples. that is, + + int inputPad = int(ceil((double(latency) * sourceRate) / targetRate)); + + // that means we are providing this much input in total: + + int n1 = n + inputPad; + + // and obtaining this much output in total: + + int m1 = int(ceil((double(n1) * targetRate) / sourceRate)); + + // in order to return this much output to the user: + + int m = int(ceil((double(n) * targetRate) / sourceRate)); + +#ifdef DEBUG_RESAMPLER + cerr << "n = " << n << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", m = " << m << ", latency = " << latency << ", inputPad = " << inputPad << ", m1 = " << m1 << ", n1 = " << n1 << ", n1 - n = " << n1 - n << endl; +#endif + + vector<double> pad(n1 - n, 0.0); + vector<double> out(m1 + 1, 0.0); + + int gotData = r.process(data, out.data(), n); + int gotPad = r.process(pad.data(), out.data() + gotData, pad.size()); + int got = gotData + gotPad; + +#ifdef DEBUG_RESAMPLER + cerr << "resample: " << n << " in, " << pad.size() << " padding, " << got << " out (" << gotData << " data, " << gotPad << " padding, latency = " << latency << ")" << endl; +#endif +#ifdef DEBUG_RESAMPLER_VERBOSE + int printN = 50; + cerr << "first " << printN << " in:" << endl; + for (int i = 0; i < printN && i < n; ++i) { + if (i % 5 == 0) cerr << endl << i << "... "; + cerr << data[i] << " "; + } + cerr << endl; +#endif + + int toReturn = got - latency; + if (toReturn > m) toReturn = m; + + vector<double> sliced(out.begin() + latency, + out.begin() + latency + toReturn); + +#ifdef DEBUG_RESAMPLER_VERBOSE + cerr << "first " << printN << " out (after latency compensation), length " << sliced.size() << ":"; + for (int i = 0; i < printN && i < sliced.size(); ++i) { + if (i % 5 == 0) cerr << endl << i << "... "; + cerr << sliced[i] << " "; + } + cerr << endl; +#endif + + return sliced; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/Resampler.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,119 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef RESAMPLER_H +#define RESAMPLER_H + +#include <vector> + +/** + * Resampler resamples a stream from one integer sample rate to + * another (arbitrary) rate, using a kaiser-windowed sinc filter. The + * results and performance are pretty similar to libraries such as + * libsamplerate, though this implementation does not support + * time-varying ratios (the ratio is fixed on construction). + * + * See also Decimator, which is faster and rougher but supports only + * power-of-two downsampling factors. + */ +class Resampler +{ +public: + /** + * Construct a Resampler to resample from sourceRate to + * targetRate. + */ + Resampler(int sourceRate, int targetRate); + + /** + * Construct a Resampler to resample from sourceRate to + * targetRate, using the given filter parameters. + */ + Resampler(int sourceRate, int targetRate, + double snr, double bandwidth); + + virtual ~Resampler(); + + /** + * Read n input samples from src and write resampled data to + * dst. The return value is the number of samples written, which + * will be no more than ceil((n * targetRate) / sourceRate). The + * caller must ensure the dst buffer has enough space for the + * samples returned. + */ + int process(const double *src, double *dst, int n); + + /** + * Read n input samples from src and return resampled data by + * value. + */ + std::vector<double> process(const double *src, int n); + + /** + * Return the number of samples of latency at the output due by + * the filter. (That is, the output will be delayed by this number + * of samples relative to the input.) + */ + int getLatency() const { return m_latency; } + + /** + * Carry out a one-off resample of a single block of n + * samples. The output is latency-compensated. + */ + static std::vector<double> resample + (int sourceRate, int targetRate, const double *data, int n); + +private: + int m_sourceRate; + int m_targetRate; + int m_gcd; + int m_filterLength; + int m_bufferLength; + int m_latency; + double m_peakToPole; + + struct Phase { + int nextPhase; + std::vector<double> filter; + int drop; + }; + + Phase *m_phaseData; + int m_phase; + std::vector<double> m_buffer; + int m_bufferOrigin; + + void initialise(double, double); + double reconstructOne(); +}; + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/SincWindow.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,63 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "SincWindow.h" + +#include <cmath> + +void +SincWindow::init() +{ + if (m_length < 1) { + return; + } else if (m_length < 2) { + m_window.push_back(1); + return; + } else { + + int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2); + int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2); + double m = 2 * M_PI / m_p; + + for (int i = 0; i < n0; ++i) { + double x = ((m_length / 2) - i) * m; + m_window.push_back(sin(x) / x); + } + + m_window.push_back(1.0); + + for (int i = 1; i < n1; ++i) { + double x = i * m; + m_window.push_back(sin(x) / x); + } + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/SincWindow.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,79 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef SINC_WINDOW_H +#define SINC_WINDOW_H + +#include <vector> + +/** + * A window containing values of the sinc function, i.e. sin(x)/x with + * sinc(0) == 1, with x == 0 at the centre. + */ +class SincWindow +{ +public: + /** + * Construct a windower of the given length, containing the values + * of sinc(x) with x=0 in the middle, i.e. at sample (length-1)/2 + * for odd or (length/2)+1 for even length, such that the distance + * from -pi to pi (the nearest zero crossings either side of the + * peak) is p samples. + */ + SincWindow(int length, double p) : m_length(length), m_p(p) { init(); } + + int getLength() const { + return m_length; + } + + const double *getWindow() const { + return m_window.data(); + } + + void cut(double *src) const { + cut(src, src); + } + + void cut(const double *src, double *dst) const { + for (int i = 0; i < m_length; ++i) { + dst[i] = src[i] * m_window[i]; + } + } + +private: + int m_length; + double m_p; + std::vector<double> m_window; + + void init(); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/Window.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,172 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef _WINDOW_H_ +#define _WINDOW_H_ + +#include <cmath> +#include <iostream> +#include <map> +#include <vector> + +enum WindowType { + RectangularWindow, + BartlettWindow, + HammingWindow, + HanningWindow, + BlackmanWindow, + BlackmanHarrisWindow, + + FirstWindow = RectangularWindow, + LastWindow = BlackmanHarrisWindow +}; + +/** + * Various shaped windows for sample frame conditioning, including + * cosine windows (Hann etc) and triangular and rectangular windows. + */ +template <typename T> +class Window +{ +public: + /** + * Construct a windower of the given type and size. + * + * Note that the cosine windows are periodic by design, rather + * than symmetrical. (A window of size N is equivalent to a + * symmetrical window of size N+1 with the final element missing.) + */ + Window(WindowType type, int size) : m_type(type), m_size(size) { encache(); } + Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); } + Window &operator=(const Window &w) { + if (&w == this) return *this; + m_type = w.m_type; + m_size = w.m_size; + encache(); + return *this; + } + virtual ~Window() { delete[] m_cache; } + + void cut(T *src) const { cut(src, src); } + void cut(const T *src, T *dst) const { + for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i]; + } + + WindowType getType() const { return m_type; } + int getSize() const { return m_size; } + + std::vector<T> getWindowData() const { + std::vector<T> d; + for (int i = 0; i < m_size; ++i) { + d.push_back(m_cache[i]); + } + return d; + } + +protected: + WindowType m_type; + int m_size; + T *m_cache; + + void encache(); +}; + +template <typename T> +void Window<T>::encache() +{ + int n = m_size; + T *mult = new T[n]; + int i; + for (i = 0; i < n; ++i) mult[i] = 1.0; + + switch (m_type) { + + case RectangularWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * 0.5; + } + break; + + case BartlettWindow: + if (n == 2) { + mult[0] = mult[1] = 0; // "matlab compatible" + } else if (n == 3) { + mult[0] = 0; + mult[1] = mult[2] = 2./3.; + } else if (n > 3) { + for (i = 0; i < n/2; ++i) { + mult[i] = mult[i] * (i / T(n/2)); + mult[i + n - n/2] = mult[i + n - n/2] * (1.0 - (i / T(n/2))); + } + } + break; + + case HammingWindow: + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n)); + } + } + break; + + case HanningWindow: + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n)); + } + } + break; + + case BlackmanWindow: + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n) + + 0.08 * cos(4 * M_PI * i / n)); + } + } + break; + + case BlackmanHarrisWindow: + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.35875 + - 0.48829 * cos(2 * M_PI * i / n) + + 0.14128 * cos(4 * M_PI * i / n) + - 0.01168 * cos(6 * M_PI * i / n)); + } + } + break; + } + + m_cache = mult; +} + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/dsp/nan-inf.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,13 @@ + +#ifndef NAN_INF_H +#define NAN_INF_H + +#define ISNAN(x) (sizeof(x) == sizeof(double) ? ISNANd(x) : ISNANf(x)) +static inline int ISNANf(float x) { return x != x; } +static inline int ISNANd(double x) { return x != x; } + +#define ISINF(x) (sizeof(x) == sizeof(double) ? ISINFd(x) : ISINFf(x)) +static inline int ISINFf(float x) { return !ISNANf(x) && ISNANf(x - x); } +static inline int ISINFd(double x) { return !ISNANd(x) && ISNANd(x - x); } + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/processfile.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,294 @@ + +#include "ConstantQ.h" +#include "CQInverse.h" + +#include <sndfile.h> + +#include <iostream> + +using std::vector; +using std::cerr; +using std::endl; + +#include <cstring> + +#include <getopt.h> +#include <unistd.h> +#include <sys/time.h> +#include <cstdlib> + +int main(int argc, char **argv) +{ + double maxFreq = 0; + double minFreq = 0; + int bpo = 0; + bool help = false; + + int c; + + while (1) { + int optionIndex = 0; + + static struct option longOpts[] = { + { "help", 0, 0, 'h', }, + { "maxfreq", 1, 0, 'x', }, + { "minfreq", 1, 0, 'n', }, + { "bpo", 1, 0, 'b' }, + { 0, 0, 0, 0 }, + }; + + c = getopt_long(argc, argv, + "hx:n:b:", + longOpts, &optionIndex); + if (c == -1) break; + + switch (c) { + case 'h': help = true; break; + case 'x': maxFreq = atof(optarg); break; + case 'n': minFreq = atof(optarg); break; + case 'b': bpo = atoi(optarg); break; + default: help = true; break; + } + } + + if (help || (optind + 2 != argc && optind + 3 != argc)) { + cerr << endl; + cerr << "Usage: " << argv[0] << " [options] infile.wav outfile.wav [differencefile.wav]" << endl; + cerr << endl; + cerr << "Options:" << endl; + cerr << " -x<X>, --maxfreq <X> Maximum frequency (default = sample rate / 3)" << endl; + cerr << " -n<X>, --minfreq <X> Minimum frequency (default = 100, actual min may vary)" << endl; + cerr << " -b<X>, --bpo <X> Bins per octave (default = 60)" << endl; + cerr << " -h, --help Print this help" << endl; + cerr << endl; + cerr << "This rather useless program simply performs a forward Constant-Q transform with" << endl; + cerr << "the requested parameters, followed by its inverse, and writes the result to the" << endl; + cerr << "output file. If a diff file name is provided, the difference between input and" << endl; + cerr << "output signals is also written to that. All this accomplishes is to produce a" << endl; + cerr << "signal that approximates the input: it's intended for test purposes only." << endl; + cerr << endl; + cerr << "(Want to calculate and obtain a Constant-Q spectrogram? Use the CQVamp plugin" << endl; + cerr << "in a Vamp plugin host.)" << endl; + cerr << endl; + return 2; + } + + char *fileName = strdup(argv[optind++]); + char *fileNameOut = strdup(argv[optind++]); + char *diffFileName = (optind < argc ? strdup(argv[optind++]) : 0); + bool doDiff = (diffFileName != 0); + + SNDFILE *sndfile; + SNDFILE *sndfileOut; + SNDFILE *sndDiffFile = 0; + SF_INFO sfinfo; + SF_INFO sfinfoOut; + SF_INFO sfinfoDiff; + 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; + } + + sfinfoOut.channels = 1; + sfinfoOut.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16; + sfinfoOut.frames = sfinfo.frames; + sfinfoOut.samplerate = sfinfo.samplerate; + 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; + } + + if (doDiff) { + sfinfoDiff = sfinfoOut; + sndDiffFile = sf_open(diffFileName, SFM_WRITE, &sfinfoDiff); + if (!sndDiffFile) { + cerr << "ERROR: Failed to open diff output file \"" << diffFileName << "\" for writing: " + << sf_strerror(sndDiffFile) << endl; + return 1; + } + } + + int ibs = 1024; + int channels = sfinfo.channels; + float *fbuf = new float[channels * ibs]; + + if (maxFreq == 0.0) maxFreq = sfinfo.samplerate / 3; + if (minFreq == 0.0) minFreq = 100; + if (bpo == 0) bpo = 60; + + ConstantQ cq(sfinfo.samplerate, minFreq, maxFreq, bpo); + CQInverse cqi(sfinfo.samplerate, minFreq, maxFreq, bpo); + + cerr << "max freq = " << cq.getMaxFrequency() << ", min freq = " + << cq.getMinFrequency() << ", octaves = " << cq.getOctaves() << endl; + + cerr << "octave boundaries: "; + for (int i = 0; i < cq.getOctaves(); ++i) { + cerr << cq.getMaxFrequency() / pow(2, i) << " "; + } + cerr << endl; + + int inframe = 0; + int outframe = 0; + int latency = cq.getLatency() + cqi.getLatency(); + + vector<double> buffer; + + double maxdiff = 0.0; + int maxdiffidx = 0; + + cerr << "forward latency = " << cq.getLatency() << ", inverse latency = " + << cqi.getLatency() << ", total = " << latency << endl; + + timeval tv; + (void)gettimeofday(&tv, 0); + + while (inframe < sfinfo.frames) { + + int count = -1; + + if ((count = sf_readf_float(sndfile, fbuf, ibs)) < 0) { + break; + } + + vector<double> cqin; + for (int i = 0; i < count; ++i) { + double v = fbuf[i * channels]; + if (channels > 1) { + for (int c = 1; c < channels; ++c) { + v += fbuf[i * channels + c]; + } + v /= channels; + } + cqin.push_back(v); + } + + if (doDiff) { + buffer.insert(buffer.end(), cqin.begin(), cqin.end()); + } + + vector<double> cqout = cqi.process(cq.process(cqin)); + + for (int i = 0; i < int(cqout.size()); ++i) { + if (cqout[i] > 1.0) cqout[i] = 1.0; + if (cqout[i] < -1.0) cqout[i] = -1.0; + } + + if (outframe >= latency) { + + sf_writef_double(sndfileOut, + cqout.data(), + cqout.size()); + + } else if (outframe + (int)cqout.size() >= latency) { + + int offset = latency - outframe; + sf_writef_double(sndfileOut, + cqout.data() + offset, + cqout.size() - offset); + } + + if (doDiff) { + for (int i = 0; i < (int)cqout.size(); ++i) { + if (outframe + i >= latency) { + int dframe = outframe + i - latency; + if (dframe >= (int)buffer.size()) cqout[i] = 0; + else cqout[i] -= buffer[dframe]; + if (fabs(cqout[i]) > maxdiff && + dframe > sfinfo.samplerate && // ignore first/last sec + dframe + sfinfo.samplerate < sfinfo.frames) { + maxdiff = fabs(cqout[i]); + maxdiffidx = dframe; + } + } + } + + if (outframe >= latency) { + + sf_writef_double(sndDiffFile, + cqout.data(), + cqout.size()); + + } else if (outframe + (int)cqout.size() >= latency) { + + int offset = latency - outframe; + sf_writef_double(sndDiffFile, + cqout.data() + offset, + cqout.size() - offset); + } + } + + inframe += count; + outframe += cqout.size(); + } + + vector<double> r = cqi.process(cq.getRemainingOutput()); + vector<double> r2 = cqi.getRemainingOutput(); + + r.insert(r.end(), r2.begin(), r2.end()); + + for (int i = 0; i < int(r.size()); ++i) { + if (r[i] > 1.0) r[i] = 1.0; + if (r[i] < -1.0) r[i] = -1.0; + } + + sf_writef_double(sndfileOut, r.data(), r.size()); + if (doDiff) { + for (int i = 0; i < (int)r.size(); ++i) { + if (outframe + i >= latency) { + int dframe = outframe + i - latency; + if (dframe >= (int)buffer.size()) r[i] = 0; + else r[i] -= buffer[dframe]; + if (fabs(r[i]) > maxdiff && + dframe > sfinfo.samplerate && // ignore first/last sec + dframe + sfinfo.samplerate < sfinfo.frames) { + maxdiff = fabs(r[i]); + maxdiffidx = dframe; + } + } + } + sf_writef_double(sndDiffFile, r.data(), r.size()); + } + outframe += r.size(); + + sf_close(sndfile); + sf_close(sndfileOut); + + if (doDiff) { + sf_close(sndDiffFile); + } + + cerr << "in: " << inframe << ", out: " << outframe - latency << endl; + + if (doDiff) { + double db = 10 * log10(maxdiff); + cerr << "max diff [excluding first and last second of audio] is " + << maxdiff << " (" << db << " dBFS)" + << " at sample index " << maxdiffidx << endl; + } + + 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 << "elapsed time (not counting init): " << sec << " sec, frames/sec at input: " << inframe/sec << endl; + + return 0; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/test.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,76 @@ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "CQSpectrogram.h" + +#include <iostream> +#include <vector> + +using std::vector; +using std::cerr; +using std::cout; +using std::endl; + +#include <cstdio> + +int main(int argc, char **argv) +{ + vector<double> in; + + for (int i = 0; i < 64; ++i) { +// if (i == 0) in.push_back(1); +// else in.push_back(0); + in.push_back(sin(i * M_PI / 2.0)); + } + + CQSpectrogram k(8, 1, 4, 4, CQSpectrogram::InterpolateZeros); + + vector<vector<double> > out = k.process(in); + vector<vector<double> > rest = k.getRemainingOutput(); + + out.insert(out.end(), rest.begin(), rest.end()); + + cerr << "got " << out.size() << " back (" << out[0].size() << " in each?)" << endl; + + for (int b = 0; b < (int)out.size() / 8; ++b) { + printf("\nColumns %d to %d:\n\n", b * 8, b * 8 + 7); + for (int j = int(out[0].size()) - 1; j >= 0; --j) { + for (int i = 0; i < 8; ++i) { + if (i + b * 8 < (int)out.size()) { + double v = out[i + b * 8][j]; + if (v < 0.0001) printf(" 0 "); + else printf(" %.4f ", out[i + b * 8][j]); + } + } + printf("\n"); + } + } +} +
--- a/vamp/CQChromaVamp.cpp Thu May 15 11:59:11 2014 +0100 +++ b/vamp/CQChromaVamp.cpp Thu May 15 14:42:11 2014 +0100 @@ -31,9 +31,9 @@ #include "CQChromaVamp.h" -#include "cpp-qm-dsp/CQSpectrogram.h" +#include "cq/CQSpectrogram.h" -#include "base/Pitch.h" +#include "Pitch.h" #include <algorithm> #include <cstdio>
--- a/vamp/CQVamp.cpp Thu May 15 11:59:11 2014 +0100 +++ b/vamp/CQVamp.cpp Thu May 15 14:42:11 2014 +0100 @@ -31,7 +31,7 @@ #include "CQVamp.h" -#include "base/Pitch.h" +#include "Pitch.h" #include <algorithm> #include <cstdio>
--- a/vamp/CQVamp.h Thu May 15 11:59:11 2014 +0100 +++ b/vamp/CQVamp.h Thu May 15 14:42:11 2014 +0100 @@ -34,7 +34,7 @@ #include <vamp-sdk/Plugin.h> -#include "cpp-qm-dsp/CQSpectrogram.h" +#include "cq/CQSpectrogram.h" class ConstantQ;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/vamp/Pitch.cpp Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,63 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#include "Pitch.h" + +#include <math.h> + +float +Pitch::getFrequencyForPitch(int midiPitch, + float centsOffset, + float concertA) +{ + float p = float(midiPitch) + (centsOffset / 100); + return concertA * powf(2.0, (p - 69.0) / 12.0); +} + +int +Pitch::getPitchForFrequency(float frequency, + float *centsOffsetReturn, + float concertA) +{ + float p = 12.0 * (log(frequency / (concertA / 2.0)) / log(2.0)) + 57.0; + + int midiPitch = int(p + 0.00001); + float centsOffset = (p - midiPitch) * 100.0; + + if (centsOffset >= 50.0) { + midiPitch = midiPitch + 1; + centsOffset = -(100.0 - centsOffset); + } + + if (centsOffsetReturn) *centsOffsetReturn = centsOffset; + return midiPitch; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/vamp/Pitch.h Thu May 15 14:42:11 2014 +0100 @@ -0,0 +1,52 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + Constant-Q library + Copyright (c) 2013-2014 Queen Mary, University of London + + Permission is hereby granted, free of charge, to any person + obtaining a copy of this software and associated documentation + files (the "Software"), to deal in the Software without + restriction, including without limitation the rights to use, copy, + modify, merge, publish, distribute, sublicense, and/or sell copies + of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY + CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + + Except as contained in this notice, the names of the Centre for + Digital Music; Queen Mary, University of London; and Chris Cannam + shall not be used in advertising or otherwise to promote the sale, + use or other dealings in this Software without prior written + authorization. +*/ + +#ifndef _PITCH_H_ +#define _PITCH_H_ + +/** + * Convert between musical pitch (i.e. MIDI pitch number) and + * fundamental frequency. + */ +class Pitch +{ +public: + static float getFrequencyForPitch(int midiPitch, + float centsOffset = 0, + float concertA = 440.0); + + static int getPitchForFrequency(float frequency, + float *centsOffsetReturn = 0, + float concertA = 440.0); +}; + + +#endif
--- a/yeti/Makefile Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - -MAY_DIR := ../../may -MAY_YC := $(MAY_DIR)/bin/yc - -test: - $(MAY_YC) ./test.yeti
--- a/yeti/build.xml Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -<project name="cqt" default="test" basedir="."> - - <property name="maydir" value="${basedir}/../../may"/> - <property name="yetidir" value="${basedir}/../../yeti"/> - - <property name="jardir" value="${maydir}/ext/jar"/> - - <property name="extjars" value="${jardir}/jvamp.jar:${jardir}/yertle.jar:${jardir}/jtransforms-2.4.jar:${jardir}/jzy3d-swt-0.9.1.jar:${jardir}/jzy3d-api-0.9.1.jar:${jardir}/jogl-all.jar:${jardir}/gluegen.jar:${jardir}/gluegen-rt.jar:${jardir}/opencsv-2.1.jar:${jardir}/org.convexhull.jar"/> - - <condition property="archtag" value="linux32"> - <os family="unix" arch="i386"/> - </condition> - <condition property="archtag" value="linux64"> - <os family="unix" arch="amd64"/> - </condition> - <condition property="archtag" value="win32"> - <os family="windows" arch="x86"/> - </condition> - <condition property="archtag" value="win64"> - <os family="windows" arch="amd64"/> - </condition> - <condition property="archtag" value="osx"> - <os family="mac"/> - </condition> - - <target name="taskdef"> - <taskdef name="yetic" classname="yeti.lang.compiler.YetiTask" - classpath="${yetidir}/yeti.jar:${maydir}/may.jar:${extjars}" /> - </target> - - <target name="prepare"> - <mkdir dir="${basedir}/classes"/> - </target> - - <target name="yeticlasses" depends="taskdef"> - <yetic srcdir="${basedir}" - destdir="${basedir}/classes" - includes="**/*.yeti" - preload="yeti/lang/std:yeti/lang/io"/> - </target> - - <target name="classes" depends="prepare,yeticlasses"/> - - <target name="jar" depends="classes,taskdef"> - <jar jarfile="${basedir}/cqt.jar"> - <fileset dir="${basedir}/classes" - includes="**/*.class" - excludes="**/test*.class"/> - </jar> - </target> - - <target name="testjar" depends="classes,taskdef"> - <jar jarfile="${basedir}/test.jar"> - <fileset dir="${basedir}/classes" - includes="**/test*.class"/> - </jar> - </target> - - <target name="test" depends="jar,testjar,taskdef"> - <java classpath="${basedir}/test.jar:${basedir}/cqt.jar:${maydir}/may.jar:${yetidir}/yeti.jar:${extjars}" - classname="test" - fork="true" failonerror="true"> - <sysproperty key="java.library.path" path="${maydir}/ext/native/${archtag}"/> - </java> - </target> - - <target name="clean"> - <delete dir="${basedir}/classes"/> - </target> - - <target name="rebuild" depends="clean,jar"/> - -</project> -
--- a/yeti/cqt.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,260 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -module cqt; - -cqtkernel = load cqtkernel; -resample = load may.stream.resample; -manipulate = load may.stream.manipulate; -mat = load may.matrix; -cm = load may.matrix.complex; -framer = load may.stream.framer; -cplx = load may.complex; -fft = load may.transform.fft; -vec = load may.vector; -ch = load may.stream.channels; - -{ pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc; - -cqt { maxFreq, minFreq, binsPerOctave } str = - (sampleRate = str.sampleRate; - octaves = ceil (log2 (maxFreq / minFreq)); -// actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave)); - - kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave }; - -// eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave), fftSize = \(kdata.fftSize), hop = \(kdata.fftHop)"; - -// eprintln "atomsPerFrame = \(kdata.atomsPerFrame)"; - - padding = (kdata.fftSize * (pow 2 (octaves-1))); - -// eprintln "padding = \(padding)"; - - str = manipulate.paddedBy padding str; - - streams = manipulate.duplicated octaves str; - - // forward transform uses the conjugate-transposed kernel, inverse - // uses the original - kernel = cm.transposed (cm.conjugateTransposed kdata.kernel); - -// eprintln "have kernel"; - - fftFunc = fft.forward kdata.fftSize; - - cqblocks = - map do octave: - frames = map ch.mixedDown //!!! mono for now - (framer.frames kdata.fftSize [ Hop kdata.fftHop, Padded false ] - (resample.decimated (pow 2 octave) streams[octave])); - map do frame: - freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize)); -// eprintln "octave = \(octave), frame = \(vec.list frame)"; -// eprintln "octave = \(octave), freq = \(freq)"; - cm.product kernel (cm.newComplexColumnVector freq); - done frames; - done [0..octaves-1]; - - // The cqblocks list is a list<list<matrix>>. Each top-level list - // corresponds to an octave, from highest to lowest, each having - // twice as many elements in its list as the next octave. The - // sub-lists are sampled in time with an effective spacing of - // fftSize * 2^(octave-1) audio frames, and the matrices are row - // vectors with atomsPerFrame * binsPerOctave complex elements. - // - // *** - // - // In a typical constant-Q structure, each (2^(octaves-1) * - // fftHop) input frames gives us an output structure conceptually - // like this: - // - // [][][][][][][][] <- fftHop frames per highest-octave output value - // [][][][][][][][] layered as many times as binsPerOctave (here 2) - // [--][--][--][--] <- fftHop*2 frames for the next lower octave - // [--][--][--][--] etc - // [------][------] - // [------][------] - // [--------------] - // [--------------] - // - // *** - // - // But the kernel we're using here has more than one temporally - // spaced atom; each individual cell is a row vector with - // atomsPerFrame * binsPerOctave elements, but that actually - // represents a rectangular matrix of result cells with width - // atomsPerFrame and height binsPerOctave. The columns of this - // matrix (the atoms) then need to be spaced by 2^(octave-1) - // relative to those from the highest octave. - - // Reshape each row vector into the appropriate rectangular matrix - // and split into single-atom columns - - emptyHops = kdata.firstCentre / kdata.atomSpacing; //!!! int? round? -// maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops; -// eprintln "maxDrop = \(maxDrop)"; - - cqblocks = - map do octlist: - concatMap do rv: - cm.asColumns - (cm.generate do row col: - cm.at rv ((row * kdata.atomsPerFrame) + col) 0 - done { - rows = kdata.binsPerOctave, - columns = kdata.atomsPerFrame - }) - done octlist - done cqblocks; - - cqblocks = array (map2 do octlist octave: - d = emptyHops * (pow 2 (octaves-octave)) - emptyHops; -// eprintln "dropping \(d)"; - drop d octlist; - done cqblocks [1..octaves]); - - assembleBlock bits = - (//eprintln "assembleBlock: structure of bits is:"; - //eprintln (map length bits); - - rows = octaves * kdata.binsPerOctave; - columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame; - - cm.generate do row col: - - // bits structure: [1,2,4,8,...] - - // each elt of bits is a list of the chunks that should - // make up this block in that octave (lowest octave first) - - // each chunk has atomsPerFrame * binsPerOctave elts in it - - // row is disposed with 0 at the top, highest octave (in - // both pitch and index into bits structure) - - oct = int (row / binsPerOctave); - binNo = row % kdata.binsPerOctave; - - chunks = pow 2 oct; - colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame)); - atomNo = int (col / colsPerAtom); - atomOffset = col % colsPerAtom; - - if atomOffset == 0 and atomNo < length bits[oct] then - bits[oct][atomNo][binNo]; - else - cplx.zero - fi; - - done { rows, columns }; - ); - - assembleBlockSpectrogram bits = - (// As assembleBlock, but producing a dense magnitude - // spectrogram (rather than a complex output). (todo: - // interpolation, smoothing) - - //eprintln "assembleBlockSpectrogram: structure of bits is:"; - //eprintln (map length bits); - - rows = octaves * kdata.binsPerOctave; - columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame; - - mat.generate do row col: - - oct = int (row / binsPerOctave); - binNo = row % kdata.binsPerOctave; - - chunks = pow 2 oct; - colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame)); - atomNo = int (col / colsPerAtom); - - if atomNo < length bits[oct] then - cplx.magnitude bits[oct][atomNo][binNo]; - else - 0 - fi; - - done { rows, columns }; - ); - - processOctaveLists assembler octs = - case octs[0] of - block::rest: - (toAssemble = array - (map do oct: - n = kdata.atomsPerFrame * pow 2 oct; - if not empty? octs[oct] then - forBlock = array (take n octs[oct]); - octs[oct] := drop n octs[oct]; - forBlock - else - array [] - fi - done (keys octs)); - assembler toAssemble :. \(processOctaveLists assembler octs)); - _: [] - esac; - -//eprintln "cqblocks has \(length cqblocks) entries"; - - octaveLists = [:]; - - cqblocks = array cqblocks; - for [1..octaves] do oct: - octaveLists[octaves - oct] := cqblocks[oct-1]; - done; -/* - \() (map2 do octlist octave: -println "oct \(octaves) - \(octave) = \(octaves - octave)"; - octaveLists[octaves - octave] := octlist - done cqblocks [1..octaves]); -*/ -//eprintln "octaveLists keys are: \(keys octaveLists)"; - - { - kernel = kdata with { - binFrequencies = array - (concatMap do octave: - map do freq: - freq / (pow 2 octave); - done (reverse (list kdata.binFrequencies)) - done [0..octaves-1]) - }, - sampleRate, - octaves, - get cqComplex () = processOctaveLists assembleBlock octaveLists, - get cqSpectrogram () = processOctaveLists assembleBlockSpectrogram octaveLists, - } - ); - -{ cqt } -
--- a/yeti/cqtkernel.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,158 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -module cqtkernel; - -vec = load may.vector; -complex = load may.complex; -window = load may.signal.window; -fft = load may.transform.fft; -cm = load may.matrix.complex; - -{ pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc; - -makeKernel { sampleRate, maxFreq, binsPerOctave } = - (q = 1; - atomHopFactor = 0.25; - thresh = 0.0005; - minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave)); - bigQ = q / ((pow 2 (1/binsPerOctave)) - 1); - - maxNK = round(bigQ * sampleRate / minFreq); - minNK = round(bigQ * sampleRate / - (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave)))); - - atomHop = round(minNK * atomHopFactor); - - firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop)); - - fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2)); - -// eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)"; -// eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)"; - - winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1; - - lastCentre = firstCentre + (winNr - 1) * atomHop; - - fftHop = (lastCentre + atomHop) - firstCentre; - -// eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)"; - - fftFunc = fft.forward fftSize; - - // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating - // function. We can't do that here; we need to generate real and imag - // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x). - - binFrequencies = array []; - - kernels = map do k: - - nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave)))); - - // the cq MATLAB toolbox uses a symmetric window for - // blackmanharris -- which is odd because it uses a periodic one - // for other types. Oh well - win = vec.divideBy nk - (vec.sqrt - (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk)); - - fk = minFreq * (pow 2 ((k-1)/binsPerOctave)); - - push binFrequencies fk; - - genKernel f = vec.multiply - [win, - vec.fromList - (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])]; - - reals = genKernel cos; - imags = genKernel sin; - - atomOffset = firstCentre - ceil(nk/2); - - map do i: - - shift = vec.zeros (atomOffset + ((i-1) * atomHop)); - - specKernel = fftFunc - (complex.complexArray - (vec.concat [shift, reals]) - (vec.concat [shift, imags])); - - map do c: - if complex.magnitude c <= thresh then complex.zero else c fi - done specKernel; - - done [1..winNr]; - - done [1..binsPerOctave]; - - kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels))); - -// eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))"; - - // Normalisation - - wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat)); - wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat)); - - subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1)); - square = cm.product (cm.conjugateTransposed subset) subset; - - diag = complex.magnitudes (cm.getDiagonal 0 square); - wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2); - - weight = (fftHop / fftSize) / (vec.mean (vec.abs wK)); - weight = sqrt(weight); - - kernel = cm.scaled weight kmat; - -// eprintln "weight = \(weight)"; - - { - kernel, - fftSize, - fftHop, - binsPerOctave, - atomsPerFrame = winNr, - atomSpacing = atomHop, - firstCentre, - maxFrequency = maxFreq, - minFrequency = minFreq, - binFrequencies, - bigQ - }); - -{ - makeKernel -} -
--- a/yeti/experiment.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,50 +0,0 @@ - -program experiment; - -af = load may.stream.audiofile; -plot = load may.plot; -cm = load may.matrix.complex; -mat = load may.matrix; -vec = load may.vector; -manipulate = load may.stream.manipulate; -syn = load may.stream.syntheticstream; - -{ cqt } = load cqt; -{ icqt } = load icqt; - -//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500); -//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4); -testStream = af.open "sweep-48000.wav"; -//testStream = af.open "sweep.wav"; - -eprintln "have test stream"; - -cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream; - -eprintln "bin frequencies: \(cq.kernel.binFrequencies)"; - -icq = icqt cq; - -eprintln "calculated icq..."; - -outfile = "test-output.wav"; -written = af.write icq outfile; - -eprintln "wrote \(written) to \(outfile)"; - -//for cq.cqSpectrogram do mm: -// for (mat.asColumns mm) (println . strJoin "," . vec.list); -//done; - -/* -bigM = mat.concatHorizontal (map cm.magnitudes cq.output); - -eprintln "overall output size = \(mat.size bigM)"; - -mat.print bigM; -*/ -//\() (plot.plot [Contour bigM]); -//\() (plot.plot [Grid bigM]); - -() -
--- a/yeti/icqt.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,145 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -module icqt; - -cqt = load cqt; -cm = load may.matrix.complex; -framer = load may.stream.framer; -win = load may.signal.window; -mm = load may.mathmisc; -vec = load may.vector; -resamp = load may.stream.resample; -manip = load may.stream.manipulate; -mat = load may.matrix; - -icqt cq = - (kdata = cq.kernel; - - // kdata.kernel is the kernel matrix for a single octave. It has - // width kdata.fftSize and height kdata.binsPerOctave * - // kdata.atomsPerFrame. - // - // kdata.fftHop is the overlap between kernel matrices in a single - // octave. - // - // cq.sampleRate is the output stream sample rate and cq.octaves - // the number of octaves. - // - // cq.cqComplex is the list of complex matrices containing the CQ - // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1) - // and height kdata.binsPerOctave * cq.octaves. - - bpo = kdata.binsPerOctave; - atomsPerFrame = kdata.atomsPerFrame; - octaves = cq.octaves; - - // transform a single block, all octaves tall, into an array - // (indexed by octave) of lists of individual columns (valid - // values for that octave only) - decomposeOctaves mat = array - (map do oct: - inverted = cm.fromRows (reverse (cm.asRows mat)); - octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1)); - gap = (mm.pow 2 oct) - 1; - pickFrom cols = - case cols of - c::cs: c::(pickFrom (drop gap cs)); - _: []; - esac; - pickFrom (cm.asColumns octMat); - done [0..octaves-1]); - - // transform a list of the arrays produced by decomposeOctaves - // into a single array (indexed by octave) of lists of the - // individual columns - flattenOctaves decomposed = - (flattenAux acc decomposed = - case decomposed of - chunk::rest: - flattenAux - (array - (map do oct: - acc[oct] ++ chunk[oct] - done [0..octaves-1])) - rest; - _: acc; - esac; - flattenAux (array (map \[] [0..octaves-1])) decomposed); - - // given a list of columns, deinterleave and pile up each sequence - // of atomsPerFrame columns into a single tall column and return - // the resulting list of tall columns - pile cols = - (pileAux acc cols = - if (length cols) >= atomsPerFrame then - atoms = take atomsPerFrame cols; - juggled = concat (reverse (cm.asRows (cm.fromColumns atoms))); - pileAux (juggled :: acc) (drop atomsPerFrame cols); - else - reverse acc - fi; - pileAux [] cols); - - octaveColumnLists = - map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex))); - - for octaveColumnLists do l: println "octave column list length: \(length l)" done; - - kernel = cm.transposed kdata.kernel; // right way around for the multiply - - spectra = - map do l: - map do col: - res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col)); - cm.columnSlice res 0 (kdata.fftSize / 2 + 1) - done l; - done octaveColumnLists; - - eprintln "calculated spectra, now to ifft, overlap-add..."; - - rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1]; - - resynthesised = - map2 do frames rate: - framer.complexStreamed rate kdata.fftSize - [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ] - frames; - done spectra rates; - - eprintln "have streams:"; - for resynthesised eprintln; - - resampled = map (resamp.resampledTo cq.sampleRate) resynthesised; - - manip.sum resampled; -); - -{icqt}
--- a/yeti/nbproject/ide-file-targets.xml Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<project basedir=".." name="cqt-IDE"> - <!-- TODO: edit the following target according to your needs --> - <!-- (more info: http://www.netbeans.org/kb/articles/freeform-config.html#debugj2se) --> - <property name="maydir" value="${basedir}/../../may"/> - <property name="yetidir" value="${basedir}/../../yeti"/> - <property name="jardir" value="${maydir}/ext/jar"/> - <path id="cp"> - <fileset dir="${basedir}"> - <include name="cqt.jar"/> - </fileset> - <fileset dir="${jardir}"> - <include name="*.jar"/> - </fileset> - <fileset dir="${yetidir}"> - <include name="yeti.jar"/> - </fileset> - <fileset dir="${maydir}"> - <include name="may.jar"/> - </fileset> - </path> - <target name="debug-nb"> - <nbjpdastart addressproperty="jpda.address" name="cqt" transport="dt_socket"> - <classpath refid="cp"/> - </nbjpdastart> - <!-- TODO configure the main class for your project here: --> - <java classname="plotfile" fork="true"> - <classpath refid="cp"/> - <jvmarg value="-Xdebug"/> - <jvmarg value="-Xrunjdwp:transport=dt_socket,address=${jpda.address}"/> - <arg value="${basedir}/test.wav"/> - </java> - </target> - <target name="-profile-check"> - <startprofiler freeform="true"/> - </target> - <!-- TODO: edit the following target according to your needs --> - <!-- (more info: http://netbeans.org/kb/articles/freeform-config.html#profilej2se) --> - <target depends="-profile-check" if="profiler.configured" name="profile-nb"> - <!-- TODO configure the main class for your project here: --> - <java classname="plotfile" fork="true"> - <classpath refid="cp"/> - <jvmarg line="${agent.jvmargs}"/> - <arg value="${basedir}/test.wav"/> - </java> - </target> -</project>
--- a/yeti/nbproject/project.xml Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<project xmlns="http://www.netbeans.org/ns/project/1"> - <type>org.netbeans.modules.ant.freeform</type> - <configuration> - <general-data xmlns="http://www.netbeans.org/ns/freeform-project/1"> - <!-- Do not use Project Properties customizer when editing this file manually. - To prevent the customizer from showing, create nbproject/project.properties file and enter -auxiliary.show.customizer=false -property there. Adding -auxiliary.show.customizer.message=<message> - will show your customized message when someone attempts to open the customizer. --> - <name>cqt</name> - <properties/> - <folders/> - <ide-actions> - <action name="build"> - <target>jar</target> - </action> - <action name="clean"> - <target>clean</target> - </action> - <action name="rebuild"> - <target>clean</target> - <target>jar</target> - </action> - <action name="debug"> - <script>nbproject/ide-file-targets.xml</script> - <target>debug-nb</target> - </action> - <action name="profile"> - <script>nbproject/ide-file-targets.xml</script> - <target>profile-nb</target> - </action> - </ide-actions> - <view> - <items> - <source-file> - <location>build.xml</location> - </source-file> - </items> - <context-menu> - <ide-action name="build"/> - <ide-action name="rebuild"/> - <ide-action name="clean"/> - <ide-action name="debug"/> - <ide-action name="profile"/> - </context-menu> - </view> - <subprojects/> - </general-data> - <java-data xmlns="http://www.netbeans.org/ns/freeform-project-java/1"/> - </configuration> -</project>
--- a/yeti/plotfile.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ - -program plotfile; - -af = load may.stream.audiofile; -plot = load may.plot; -cm = load may.matrix.complex; -mat = load may.matrix; - -{ cqt } = load cqt; - -minFreq = 27.5; -maxRateDivisor = 3; -binsPerOctave = 60; - -plotfile f = - (testStream = af.open f; - eprintln "Opened file stream..."; - start = System#currentTimeMillis(); - cq = cqt { - maxFreq = testStream.sampleRate/maxRateDivisor, - minFreq, - binsPerOctave - } testStream; - middle = System#currentTimeMillis(); - eprintln "Generated kernel \(cm.size cq.kernel.kernel) and primed transform (\(cq.octaves) octaves), took \(middle-start)ms, now calculating..."; - bigM = mat.concatHorizontal cq.cqSpectrogram; - finish = System#currentTimeMillis(); - eprintln "Done, that part took \(finish-middle)ms, all CQ stuff took \(finish-start)ms"; - eprintln "Plotting..."; - \() (plot.plot [Contour bigM])); - -usage () = - (eprintln "\nUsage: plotfile file.wav"; - eprintln "\n Loads audio from file.wav and plots a \(binsPerOctave)bpo Constant-Q spectrogram"; - eprintln " from \(minFreq)Hz up to 1/\(maxRateDivisor) of the file samplerate"); - -case (list _argv) of -file::[]: plotfile file; -_: usage (); -esac; -
--- a/yeti/test.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,96 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -program test; - -{ runTests } = load may.test; - -// We want to test: -// -// Kernel design -- check size (number of bins, number of atoms); -// check an example kernel against known data -// -// Time alignment -- feed a dirac train, check that peaks in all bins -// align -// -// Frequency discrimination -- feed a sinusoid, check peaks -// -// Latency compensation -- for dirac at 0, check peak can be found at -// 0 plus the declared latency -// -// Signal-noise ratio -// -// Specimen output for simple test case - -tests = [ -"kernel" : load test_cqtkernel, -"frequency" : load test_frequency, -]; - -bad = sum (mapHash do name testHash: runTests name testHash done tests); - -if (bad > 0) then - println "\n** \(bad) test(s) failed!"; - threadExit 1 -else - () -fi - - - - -/* -//testStream = manipulate.withDuration 96000 (syn.sinusoid 48000 500); -//testStream = manipulate.withDuration 96000 (syn.pulseTrain 48000 4); -testStream = af.open "sweep-48000.wav"; -//testStream = af.open "sweep.wav"; - -// So the stream is [ 0, 1, 0, -1, 0, 1, 0, -1, ... ] : -//testStream = manipulate.withDuration 64 (syn.sinusoid 8 2); - -testStream = manipulate.withDuration 32 (syn.pulseTrain 8 0.001); - -eprintln "have test stream"; - -cq = cqt { maxFreq = testStream.sampleRate/2, minFreq = 50, binsPerOctave = 24 } testStream; - -eprintln "bin frequencies: \(cq.kernel.binFrequencies)"; - -bigM = mat.concatHorizontal (map cm.magnitudes cq.output); - -eprintln "overall output size = \(mat.size bigM)"; - -mat.print bigM; - -//\() (plot.plot [Contour bigM]); -\() (plot.plot [Grid bigM]); -*/ -//() -
--- a/yeti/test_cqtkernel.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,63 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -module test_cqtkernel; - -cm = load may.matrix.complex; -mm = load may.mathmisc; - -{ compare, compareUsing } = load may.test; - -{ makeKernel } = load cqtkernel; - -eps = 1e-7; - -compareClose = compareUsing do a b: abs (a - b) < eps done; - -[ - -"minimal": \( - k = makeKernel { sampleRate = 16, maxFreq = 8, binsPerOctave = 4 }; - compare k.binsPerOctave 4 and - compare (cm.size k.kernel) { - rows = k.binsPerOctave * k.atomsPerFrame, - columns = k.fftSize - } and - compareClose k.maxFrequency 8 and - compareClose k.minFrequency (4 * (mm.pow 2 (1/4))) and - compare k.atomsPerFrame 5 and - compare k.fftSize 32 and - compare (length k.binFrequencies) k.binsPerOctave and - compareClose (head k.binFrequencies) k.minFrequency and - compareClose (head (reverse k.binFrequencies)) k.maxFrequency -), - -] is hash<string, () -> boolean> -
--- a/yeti/test_frequency.yeti Thu May 15 11:59:11 2014 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,132 +0,0 @@ -/* - Constant-Q library - Copyright (c) 2013-2014 Queen Mary, University of London - - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, copy, - modify, merge, publish, distribute, sublicense, and/or sell copies - of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF - CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - - Except as contained in this notice, the names of the Centre for - Digital Music; Queen Mary, University of London; and Chris Cannam - shall not be used in advertising or otherwise to promote the sale, - use or other dealings in this Software without prior written - authorization. -*/ - -module test_frequency; - -mat = load may.matrix; -vec = load may.vector; -win = load may.signal.window; -mm = load may.mathmisc; -cm = load may.matrix.complex; -syn = load may.stream.syntheticstream; -plot = load may.plot; - -{ compare } = load may.test; - -{ cqt } = load cqt; - -// Test with a single windowed sinusoid, repeating at various frequencies - -sinTestStream sampleRate duration signalFreq = // duration is in samples - (sin = syn.sinusoid sampleRate signalFreq; - chunk = mat.getRow 0 (sin.read duration); - syn.precalculatedMono sampleRate (win.windowed win.hann chunk)); - -// We want to make a CQ transform spanning more than one octave, but -// not going all the way to fs/2 so we can test it also with -// frequencies above and below its extents - -sampleRate = 100; - -// fs/2 = 50 so 10->40 gives us 2 octaves -cqmin = 10; -cqmax = 40; -bpo = 4; // fairly arbitrary - -testFreqs = map (* 5) [ 0..10 ]; -duration = sampleRate * 2; - -threshold = 0.08; - -streamBuilder = sinTestStream sampleRate duration; - -binForFreq f = - mm.round (bpo * mm.log2 (f / cqmin)) - 1; - -report message matrix = - (eprintln message; - eprintln "matrix is:"; - mat.eprint matrix); -// chart = plot.plot [Grid matrix]; -// sleep 100; -// chart#dispose()); - -tests = mapIntoHash - do f: "freq_\(f)" done - do f: \( - str = streamBuilder f; - cq = cqt { maxFreq = cqmax, minFreq = cqmin, binsPerOctave = bpo } str; - spec = cq.cqSpectrogram; - rightSize = all id - (map do s: - compare (mat.size s) { - rows = cq.kernel.binsPerOctave * cq.octaves, - columns = cq.kernel.atomsPerFrame * mm.pow 2 (cq.octaves - 1) - } - done spec); - m = mat.concatHorizontal spec; -// println "binFrequencies = \(cq.kernel.binFrequencies)"; -// println "binForFreq \(f) = \(binForFreq f)"; - var colno = 0; - success = all id - (rightSize :: map do c: - // The test passes for this column if: - // - // * the max bin is the expected one, or - // - // * the expected max is out of range entirely (but - // we need to test _something_ in this case -- - // what?), or - // - // * all bins are below a threshold, or - // - // * this is an odd column and the expected max is in - // the lower octave - // - // We should also check that all values in the lower - // octave are zero for odd columns. - // - expected = binForFreq f; - good = - (expected < 0 or expected >= vec.length c) or - ((colno % 2 == 1) and expected < (vec.length c / 2)) or - (vec.max c < threshold) or - (vec.maxindex c == binForFreq f); - if (not good) then - report " * bad! maxindex \(vec.maxindex c) != expected \(binForFreq f) for freq \(f) in column \(colno) of \(mat.width m): column is \(vec.list c)" m; - fi; - colno := colno + 1; - good; - done (mat.asColumns m)); - success; - ) done - testFreqs; - -tests is hash<string, () -> boolean>