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