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