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 cqt;
c@116: 
c@116: cqtkernel = load cqtkernel;
c@116: resample = load may.stream.resample;
c@116: manipulate = load may.stream.manipulate;
c@116: mat = load may.matrix;
c@116: cm = load may.matrix.complex;
c@116: framer = load may.stream.framer;
c@116: cplx = load may.complex;
c@116: fft = load may.transform.fft;
c@116: vec = load may.vector;
c@116: ch = load may.stream.channels;
c@116: 
c@116: { pow, round, floor, ceil, log2, nextPowerOfTwo } = load may.mathmisc;
c@116: 
c@116: cqt { maxFreq, minFreq, binsPerOctave } str =
c@116:    (sampleRate = str.sampleRate;
c@116:     octaves = ceil (log2 (maxFreq / minFreq));
c@116: //    actualMinFreq = (maxFreq / (pow 2 octaves)) * (pow 2 (1/binsPerOctave));
c@116: 
c@116:     kdata = cqtkernel.makeKernel { sampleRate, maxFreq, binsPerOctave };
c@116: 
c@116: //    eprintln "sampleRate = \(sampleRate), maxFreq = \(maxFreq), minFreq = \(minFreq), actualMinFreq = \(actualMinFreq), octaves = \(octaves), binsPerOctave = \(binsPerOctave), fftSize = \(kdata.fftSize), hop = \(kdata.fftHop)";
c@116: 
c@116: //    eprintln "atomsPerFrame = \(kdata.atomsPerFrame)";
c@116: 
c@116:     padding = (kdata.fftSize * (pow 2 (octaves-1)));
c@116: 
c@116: //    eprintln "padding = \(padding)";
c@116: 
c@116:     str = manipulate.paddedBy padding str;
c@116: 
c@116:     streams = manipulate.duplicated octaves str;
c@116: 
c@116:     // forward transform uses the conjugate-transposed kernel, inverse
c@116:     // uses the original
c@116:     kernel = cm.transposed (cm.conjugateTransposed kdata.kernel);
c@116: 
c@116: //    eprintln "have kernel";
c@116: 
c@116:     fftFunc = fft.forward kdata.fftSize;
c@116: 
c@116:     cqblocks =
c@116:         map do octave:
c@116:             frames = map ch.mixedDown //!!! mono for now
c@116:                (framer.frames kdata.fftSize [ Hop kdata.fftHop, Padded false ]
c@116:                    (resample.decimated (pow 2 octave) streams[octave]));
c@116:             map do frame:
c@116:                 freq = fftFunc (cplx.complexArray frame (vec.zeros kdata.fftSize));
c@116: // eprintln "octave = \(octave), frame = \(vec.list frame)";
c@116: // eprintln "octave = \(octave), freq = \(freq)";
c@116:                 cm.product kernel (cm.newComplexColumnVector freq);
c@116:             done frames;
c@116:         done [0..octaves-1];
c@116: 
c@116:     // The cqblocks list is a list<list<matrix>>. Each top-level list
c@116:     // corresponds to an octave, from highest to lowest, each having
c@116:     // twice as many elements in its list as the next octave. The
c@116:     // sub-lists are sampled in time with an effective spacing of
c@116:     // fftSize * 2^(octave-1) audio frames, and the matrices are row
c@116:     // vectors with atomsPerFrame * binsPerOctave complex elements.
c@116:     //
c@116:     // ***
c@116:     // 
c@116:     // In a typical constant-Q structure, each (2^(octaves-1) *
c@116:     // fftHop) input frames gives us an output structure conceptually
c@116:     // like this:
c@116:     //
c@116:     // [][][][][][][][]   <- fftHop frames per highest-octave output value
c@116:     // [][][][][][][][]      layered as many times as binsPerOctave (here 2)
c@116:     // [--][--][--][--]   <- fftHop*2 frames for the next lower octave
c@116:     // [--][--][--][--]      etc
c@116:     // [------][------]
c@116:     // [------][------]
c@116:     // [--------------]
c@116:     // [--------------]
c@116:     //
c@116:     // ***
c@116:     //
c@116:     // But the kernel we're using here has more than one temporally
c@116:     // spaced atom; each individual cell is a row vector with
c@116:     // atomsPerFrame * binsPerOctave elements, but that actually
c@116:     // represents a rectangular matrix of result cells with width
c@116:     // atomsPerFrame and height binsPerOctave. The columns of this
c@116:     // matrix (the atoms) then need to be spaced by 2^(octave-1)
c@116:     // relative to those from the highest octave.
c@116: 
c@116:     // Reshape each row vector into the appropriate rectangular matrix
c@116:     // and split into single-atom columns
c@116: 
c@116:     emptyHops = kdata.firstCentre / kdata.atomSpacing; //!!! int? round?
c@116: //    maxDrop = emptyHops * (pow 2 (octaves-1)) - emptyHops;
c@116: //    eprintln "maxDrop = \(maxDrop)";
c@116: 
c@116:     cqblocks =
c@116:         map do octlist:
c@116:             concatMap do rv:
c@116:                 cm.asColumns
c@116:                    (cm.generate do row col:
c@116:                         cm.at rv ((row * kdata.atomsPerFrame) + col) 0
c@116:                     done {
c@116:                         rows = kdata.binsPerOctave,
c@116:                         columns = kdata.atomsPerFrame
c@116:                     })
c@116:             done octlist
c@116:         done cqblocks;
c@116: 
c@116:     cqblocks = array (map2 do octlist octave:
c@116:         d = emptyHops * (pow 2 (octaves-octave)) - emptyHops;
c@116: //        eprintln "dropping \(d)";
c@116:         drop d octlist;
c@116:     done cqblocks [1..octaves]);
c@116: 
c@116:     assembleBlock bits =
c@116:        (//eprintln "assembleBlock: structure of bits is:";
c@116:         //eprintln (map length bits);
c@116: 
c@116:         rows = octaves * kdata.binsPerOctave;
c@116:         columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
c@116: 
c@116:         cm.generate do row col:
c@116: 
c@116:             // bits structure: [1,2,4,8,...]
c@116: 
c@116:             // each elt of bits is a list of the chunks that should
c@116:             // make up this block in that octave (lowest octave first)
c@116: 
c@116:             // each chunk has atomsPerFrame * binsPerOctave elts in it
c@116: 
c@116:             // row is disposed with 0 at the top, highest octave (in
c@116:             // both pitch and index into bits structure)
c@116: 
c@116:             oct = int (row / binsPerOctave);
c@116:             binNo = row % kdata.binsPerOctave;
c@116: 
c@116:             chunks = pow 2 oct;
c@116:             colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
c@116:             atomNo = int (col / colsPerAtom);
c@116:             atomOffset = col % colsPerAtom;
c@116: 
c@116:             if atomOffset == 0 and atomNo < length bits[oct] then
c@116:                 bits[oct][atomNo][binNo];
c@116:             else
c@116:                 cplx.zero
c@116:             fi;
c@116: 
c@116:         done { rows, columns };
c@116:         );
c@116: 
c@116:     assembleBlockSpectrogram bits =
c@116:        (// As assembleBlock, but producing a dense magnitude
c@116:         // spectrogram (rather than a complex output). (todo:
c@116:         // interpolation, smoothing)
c@116: 
c@116:         //eprintln "assembleBlockSpectrogram: structure of bits is:";
c@116:         //eprintln (map length bits);
c@116: 
c@116:         rows = octaves * kdata.binsPerOctave;
c@116:         columns = (pow 2 (octaves - 1)) * kdata.atomsPerFrame;
c@116: 
c@116:         mat.generate do row col:
c@116: 
c@116:             oct = int (row / binsPerOctave);
c@116:             binNo = row % kdata.binsPerOctave;
c@116: 
c@116:             chunks = pow 2 oct;
c@116:             colsPerAtom = int (columns / (chunks * kdata.atomsPerFrame));
c@116:             atomNo = int (col / colsPerAtom);
c@116: 
c@116:             if atomNo < length bits[oct] then
c@116:                 cplx.magnitude bits[oct][atomNo][binNo];
c@116:             else 
c@116:                 0
c@116:             fi;
c@116: 
c@116:         done { rows, columns };
c@116:         );
c@116: 
c@116:     processOctaveLists assembler octs =
c@116:         case octs[0] of
c@116:         block::rest:
c@116:            (toAssemble = array 
c@116:                (map do oct:
c@116:                     n = kdata.atomsPerFrame * pow 2 oct;
c@116:                     if not empty? octs[oct] then
c@116:                         forBlock = array (take n octs[oct]);
c@116:                         octs[oct] := drop n octs[oct];
c@116:                         forBlock
c@116:                     else
c@116:                         array []
c@116:                     fi
c@116:                 done (keys octs));
c@116:             assembler toAssemble :. \(processOctaveLists assembler octs));
c@116:          _: []
c@116:         esac;
c@116: 
c@116: //eprintln "cqblocks has \(length cqblocks) entries";
c@116: 
c@116:     octaveLists = [:];
c@116: 
c@116:     cqblocks = array cqblocks;
c@116:     for [1..octaves] do oct:
c@116:         octaveLists[octaves - oct] := cqblocks[oct-1];
c@116:     done;
c@116: /*
c@116:     \() (map2 do octlist octave:
c@116: println "oct \(octaves) - \(octave) = \(octaves - octave)";
c@116:              octaveLists[octaves - octave] := octlist 
c@116:          done cqblocks [1..octaves]);
c@116: */
c@116: //eprintln "octaveLists keys are: \(keys octaveLists)";
c@116: 
c@116:     {
c@116:         kernel = kdata with {
c@116:             binFrequencies = array
c@116:                (concatMap do octave:
c@116:                     map do freq:
c@116:                         freq / (pow 2 octave);
c@116:                     done (reverse (list kdata.binFrequencies))
c@116:                 done [0..octaves-1])
c@116:         },
c@116:         sampleRate,
c@116:         octaves,
c@116:         get cqComplex () = processOctaveLists assembleBlock octaveLists,
c@116:         get cqSpectrogram () = processOctaveLists assembleBlockSpectrogram octaveLists,
c@116:     }
c@116:     );
c@116: 
c@116: { cqt }
c@116: