c@69: /* c@69: Constant-Q library c@69: Copyright (c) 2013-2014 Queen Mary, University of London c@69: c@69: Permission is hereby granted, free of charge, to any person c@69: obtaining a copy of this software and associated documentation c@69: files (the "Software"), to deal in the Software without c@69: restriction, including without limitation the rights to use, copy, c@69: modify, merge, publish, distribute, sublicense, and/or sell copies c@69: of the Software, and to permit persons to whom the Software is c@69: furnished to do so, subject to the following conditions: c@69: c@69: The above copyright notice and this permission notice shall be c@69: included in all copies or substantial portions of the Software. c@69: c@69: THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, c@69: EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF c@69: MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND c@69: NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY c@69: CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF c@69: CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION c@69: WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. c@69: c@69: Except as contained in this notice, the names of the Centre for c@69: Digital Music; Queen Mary, University of London; and Chris Cannam c@69: shall not be used in advertising or otherwise to promote the sale, c@69: use or other dealings in this Software without prior written c@69: authorization. c@69: */ c@1: c@1: module cqtkernel; c@1: c@3: vec = load may.vector; c@3: complex = load may.complex; c@3: window = load may.signal.window; c@3: fft = load may.transform.fft; c@6: cm = load may.matrix.complex; c@3: c@2: { pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc; c@1: c@9: makeKernel { sampleRate, maxFreq, binsPerOctave } = c@9: (q = 1; c@9: atomHopFactor = 0.25; c@9: thresh = 0.0005; c@9: minFreq = (maxFreq/2) * (pow 2 (1/binsPerOctave)); c@9: bigQ = q / ((pow 2 (1/binsPerOctave)) - 1); c@1: c@9: maxNK = round(bigQ * sampleRate / minFreq); c@9: minNK = round(bigQ * sampleRate / c@9: (minFreq * (pow 2 ((binsPerOctave-1) / binsPerOctave)))); c@1: c@9: atomHop = round(minNK * atomHopFactor); c@9: c@9: firstCentre = atomHop * (ceil ((ceil (maxNK/2)) / atomHop)); c@9: c@9: fftSize = nextPowerOfTwo (firstCentre + ceil (maxNK/2)); c@9: c@64: // eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), binsPerOctave = \(binsPerOctave), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)"; c@64: // eprintln "minFreq = \(minFreq), bigQ = \(bigQ), maxNK = \(maxNK), minNK = \(minNK), atomHop = \(atomHop), firstCentre = \(firstCentre), fftSize = \(fftSize)"; c@9: c@9: winNr = floor((fftSize - ceil(maxNK/2) - firstCentre) / atomHop) + 1; c@9: c@9: lastCentre = firstCentre + (winNr - 1) * atomHop; c@9: c@9: fftHop = (lastCentre + atomHop) - firstCentre; c@9: c@64: // eprintln "winNr = \(winNr), lastCentre = \(lastCentre), fftHop = \(fftHop)"; c@9: c@9: fftFunc = fft.forward fftSize; c@9: c@9: // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating c@9: // function. We can't do that here; we need to generate real and imag c@9: // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x). c@9: c@40: binFrequencies = array []; c@40: c@9: kernels = map do k: c@9: c@9: nk = round(bigQ * sampleRate / (minFreq * (pow 2 ((k-1)/binsPerOctave)))); c@23: c@9: // the cq MATLAB toolbox uses a symmetric window for c@9: // blackmanharris -- which is odd because it uses a periodic one c@9: // for other types. Oh well c@25: win = vec.divideBy nk c@25: (vec.sqrt c@9: (window.windowFunction (BlackmanHarris ()) [Symmetric true] nk)); c@23: c@9: fk = minFreq * (pow 2 ((k-1)/binsPerOctave)); c@23: c@40: push binFrequencies fk; c@40: c@25: genKernel f = vec.multiply win c@9: (vec.fromList c@9: (map do i: f (2 * pi * fk * i / sampleRate) done [0..nk-1])); c@9: c@9: reals = genKernel cos; c@9: imags = genKernel sin; c@9: c@9: atomOffset = firstCentre - ceil(nk/2); c@9: c@9: map do i: c@9: c@9: shift = vec.zeros (atomOffset + ((i-1) * atomHop)); c@9: c@9: specKernel = fftFunc c@9: (complex.complexArray c@9: (vec.concat [shift, reals]) c@9: (vec.concat [shift, imags])); c@23: c@9: map do c: c@9: if complex.magnitude c <= thresh then complex.zero else c fi c@9: done specKernel; c@23: c@9: done [1..winNr]; c@9: c@9: done [1..binsPerOctave]; c@9: c@43: kmat = cm.toSparse (cm.scaled (1/fftSize) (cm.fromRows (concat kernels))); c@9: c@64: // eprintln "density = \(cm.density kmat) (\(cm.nonZeroValues kmat) of \(cm.width kmat * cm.height kmat))"; c@9: c@9: // Normalisation c@9: c@25: wx1 = vec.maxindex (complex.magnitudes (cm.getRow 0 kmat)); c@25: wx2 = vec.maxindex (complex.magnitudes (cm.getRow (cm.height kmat - 1) kmat)); c@27: c@72: subset = cm.flipped (cm.columnSlice kmat wx1 (wx2+1)); c@9: square = cm.product (cm.conjugateTransposed subset) subset; c@27: c@9: diag = complex.magnitudes (cm.getDiagonal 0 square); c@9: wK = vec.slice diag (round(1/q)) (vec.length diag - round(1/q) - 2); c@27: c@25: weight = (fftHop / fftSize) / (vec.mean (vec.abs wK)); c@9: weight = sqrt(weight); c@1: c@64: // eprintln "weight = \(weight)"; c@23: c@9: { c@9: kernel = cm.scaled weight kmat, c@9: fftSize, c@9: fftHop, c@9: binsPerOctave, c@12: atomsPerFrame = winNr, c@12: atomSpacing = atomHop, c@13: firstCentre, c@40: maxFrequency = maxFreq, c@40: minFrequency = minFreq, c@40: binFrequencies, c@9: bigQ c@9: }); c@1: c@9: { c@9: makeKernel c@9: } c@1: