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