Chris@366: /* Chris@366: Constant-Q library Chris@366: Copyright (c) 2013-2014 Queen Mary, University of London Chris@366: Chris@366: Permission is hereby granted, free of charge, to any person Chris@366: obtaining a copy of this software and associated documentation Chris@366: files (the "Software"), to deal in the Software without Chris@366: restriction, including without limitation the rights to use, copy, Chris@366: modify, merge, publish, distribute, sublicense, and/or sell copies Chris@366: of the Software, and to permit persons to whom the Software is Chris@366: furnished to do so, subject to the following conditions: Chris@366: Chris@366: The above copyright notice and this permission notice shall be Chris@366: included in all copies or substantial portions of the Software. Chris@366: Chris@366: THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, Chris@366: EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF Chris@366: MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND Chris@366: NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY Chris@366: CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF Chris@366: CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION Chris@366: WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. Chris@366: Chris@366: Except as contained in this notice, the names of the Centre for Chris@366: Digital Music; Queen Mary, University of London; and Chris Cannam Chris@366: shall not be used in advertising or otherwise to promote the sale, Chris@366: use or other dealings in this Software without prior written Chris@366: authorization. Chris@366: */ Chris@366: Chris@366: module icqt; Chris@366: Chris@366: cqt = load cqt; Chris@366: cm = load may.matrix.complex; Chris@366: framer = load may.stream.framer; Chris@366: win = load may.signal.window; Chris@366: mm = load may.mathmisc; Chris@366: vec = load may.vector; Chris@366: resamp = load may.stream.resample; Chris@366: manip = load may.stream.manipulate; Chris@366: mat = load may.matrix; Chris@366: Chris@366: icqt cq = Chris@366: (kdata = cq.kernel; Chris@366: Chris@366: // kdata.kernel is the kernel matrix for a single octave. It has Chris@366: // width kdata.fftSize and height kdata.binsPerOctave * Chris@366: // kdata.atomsPerFrame. Chris@366: // Chris@366: // kdata.fftHop is the overlap between kernel matrices in a single Chris@366: // octave. Chris@366: // Chris@366: // cq.sampleRate is the output stream sample rate and cq.octaves Chris@366: // the number of octaves. Chris@366: // Chris@366: // cq.cqComplex is the list of complex matrices containing the CQ Chris@366: // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1) Chris@366: // and height kdata.binsPerOctave * cq.octaves. Chris@366: Chris@366: bpo = kdata.binsPerOctave; Chris@366: atomsPerFrame = kdata.atomsPerFrame; Chris@366: octaves = cq.octaves; Chris@366: Chris@366: // transform a single block, all octaves tall, into an array Chris@366: // (indexed by octave) of lists of individual columns (valid Chris@366: // values for that octave only) Chris@366: decomposeOctaves mat = array Chris@366: (map do oct: Chris@366: inverted = cm.fromRows (reverse (cm.asRows mat)); Chris@366: octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1)); Chris@366: gap = (mm.pow 2 oct) - 1; Chris@366: pickFrom cols = Chris@366: case cols of Chris@366: c::cs: c::(pickFrom (drop gap cs)); Chris@366: _: []; Chris@366: esac; Chris@366: pickFrom (cm.asColumns octMat); Chris@366: done [0..octaves-1]); Chris@366: Chris@366: // transform a list of the arrays produced by decomposeOctaves Chris@366: // into a single array (indexed by octave) of lists of the Chris@366: // individual columns Chris@366: flattenOctaves decomposed = Chris@366: (flattenAux acc decomposed = Chris@366: case decomposed of Chris@366: chunk::rest: Chris@366: flattenAux Chris@366: (array Chris@366: (map do oct: Chris@366: acc[oct] ++ chunk[oct] Chris@366: done [0..octaves-1])) Chris@366: rest; Chris@366: _: acc; Chris@366: esac; Chris@366: flattenAux (array (map \[] [0..octaves-1])) decomposed); Chris@366: Chris@366: // given a list of columns, deinterleave and pile up each sequence Chris@366: // of atomsPerFrame columns into a single tall column and return Chris@366: // the resulting list of tall columns Chris@366: pile cols = Chris@366: (pileAux acc cols = Chris@366: if (length cols) >= atomsPerFrame then Chris@366: atoms = take atomsPerFrame cols; Chris@366: juggled = concat (reverse (cm.asRows (cm.fromColumns atoms))); Chris@366: pileAux (juggled :: acc) (drop atomsPerFrame cols); Chris@366: else Chris@366: reverse acc Chris@366: fi; Chris@366: pileAux [] cols); Chris@366: Chris@366: octaveColumnLists = Chris@366: map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex))); Chris@366: Chris@366: for octaveColumnLists do l: println "octave column list length: \(length l)" done; Chris@366: Chris@366: kernel = cm.transposed kdata.kernel; // right way around for the multiply Chris@366: Chris@366: spectra = Chris@366: map do l: Chris@366: map do col: Chris@366: res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col)); Chris@366: cm.columnSlice res 0 (kdata.fftSize / 2 + 1) Chris@366: done l; Chris@366: done octaveColumnLists; Chris@366: Chris@366: eprintln "calculated spectra, now to ifft, overlap-add..."; Chris@366: Chris@366: rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1]; Chris@366: Chris@366: resynthesised = Chris@366: map2 do frames rate: Chris@366: framer.complexStreamed rate kdata.fftSize Chris@366: [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ] Chris@366: frames; Chris@366: done spectra rates; Chris@366: Chris@366: eprintln "have streams:"; Chris@366: for resynthesised eprintln; Chris@366: Chris@366: resampled = map (resamp.resampledTo cq.sampleRate) resynthesised; Chris@366: Chris@366: manip.sum resampled; Chris@366: ); Chris@366: Chris@366: {icqt}