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 cqtkernel; Chris@366: Chris@366: vec = load may.vector; Chris@366: complex = load may.complex; Chris@366: window = load may.signal.window; Chris@366: fft = load may.transform.fft; Chris@366: cm = load may.matrix.complex; Chris@366: Chris@366: { pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc; Chris@366: Chris@366: makeKernel { sampleRate, maxFreq, binsPerOctave } = Chris@366: (q = 1; Chris@366: atomHopFactor = 0.25; Chris@366: thresh = 0.0005; Chris@366: minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave)); Chris@366: bigQ = q / ((pow 2 (1/binsPerOctave)) - 1); Chris@366: Chris@366: maxNK = round(bigQ * sampleRate / minFreq); Chris@366: minNK = round(bigQ * sampleRate / Chris@366: (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave)))); Chris@366: Chris@366: atomHop = round(minNK * atomHopFactor); Chris@366: Chris@366: firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop)); Chris@366: Chris@366: fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2)); Chris@366: Chris@366: // eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)"; Chris@366: // eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)"; Chris@366: Chris@366: winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1; Chris@366: Chris@366: lastCentre = firstCentre + (winNr - 1) * atomHop; Chris@366: Chris@366: fftHop = (lastCentre + atomHop) - firstCentre; Chris@366: Chris@366: // eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)"; Chris@366: Chris@366: fftFunc = fft.forward fftSize; Chris@366: Chris@366: // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating Chris@366: // function. We can't do that here; we need to generate real and imag Chris@366: // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x). Chris@366: Chris@366: binFrequencies = array []; Chris@366: Chris@366: kernels = map do k: Chris@366: Chris@366: nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave)))); Chris@366: Chris@366: // the cq MATLAB toolbox uses a symmetric window for Chris@366: // blackmanharris -- which is odd because it uses a periodic one Chris@366: // for other types. Oh well Chris@366: win = vec.divideBy nk Chris@366: (vec.sqrt Chris@366: (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk)); Chris@366: Chris@366: fk = minFreq * (pow 2 ((k-1)/binsPerOctave)); Chris@366: Chris@366: push binFrequencies fk; Chris@366: Chris@366: genKernel f = vec.multiply Chris@366: [win, Chris@366: vec.fromList Chris@366: (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])]; Chris@366: Chris@366: reals = genKernel cos; Chris@366: imags = genKernel sin; Chris@366: Chris@366: atomOffset = firstCentre - ceil(nk/2); Chris@366: Chris@366: map do i: Chris@366: Chris@366: shift = vec.zeros (atomOffset + ((i-1) * atomHop)); Chris@366: Chris@366: specKernel = fftFunc Chris@366: (complex.complexArray Chris@366: (vec.concat [shift, reals]) Chris@366: (vec.concat [shift, imags])); Chris@366: Chris@366: map do c: Chris@366: if complex.magnitude c <= thresh then complex.zero else c fi Chris@366: done specKernel; Chris@366: Chris@366: done [1..winNr]; Chris@366: Chris@366: done [1..binsPerOctave]; Chris@366: Chris@366: kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels))); Chris@366: Chris@366: // eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))"; Chris@366: Chris@366: // Normalisation Chris@366: Chris@366: wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat)); Chris@366: wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat)); Chris@366: Chris@366: subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1)); Chris@366: square = cm.product (cm.conjugateTransposed subset) subset; Chris@366: Chris@366: diag = complex.magnitudes (cm.getDiagonal 0 square); Chris@366: wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2); Chris@366: Chris@366: weight = (fftHop / fftSize) / (vec.mean (vec.abs wK)); Chris@366: weight = sqrt(weight); Chris@366: Chris@366: kernel = cm.scaled weight kmat; Chris@366: Chris@366: // eprintln "weight = \(weight)"; Chris@366: Chris@366: { Chris@366: kernel, Chris@366: fftSize, Chris@366: fftHop, Chris@366: binsPerOctave, Chris@366: atomsPerFrame = winNr, Chris@366: atomSpacing = atomHop, Chris@366: firstCentre, Chris@366: maxFrequency = maxFreq, Chris@366: minFrequency = minFreq, Chris@366: binFrequencies, Chris@366: bigQ Chris@366: }); Chris@366: Chris@366: { Chris@366: makeKernel Chris@366: } Chris@366: