c@116
|
1 /*
|
c@116
|
2 Constant-Q library
|
c@116
|
3 Copyright (c) 2013-2014 Queen Mary, University of London
|
c@116
|
4
|
c@116
|
5 Permission is hereby granted, free of charge, to any person
|
c@116
|
6 obtaining a copy of this software and associated documentation
|
c@116
|
7 files (the "Software"), to deal in the Software without
|
c@116
|
8 restriction, including without limitation the rights to use, copy,
|
c@116
|
9 modify, merge, publish, distribute, sublicense, and/or sell copies
|
c@116
|
10 of the Software, and to permit persons to whom the Software is
|
c@116
|
11 furnished to do so, subject to the following conditions:
|
c@116
|
12
|
c@116
|
13 The above copyright notice and this permission notice shall be
|
c@116
|
14 included in all copies or substantial portions of the Software.
|
c@116
|
15
|
c@116
|
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
c@116
|
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
c@116
|
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
c@116
|
19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
|
c@116
|
20 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
|
c@116
|
21 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
|
c@116
|
22 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
c@116
|
23
|
c@116
|
24 Except as contained in this notice, the names of the Centre for
|
c@116
|
25 Digital Music; Queen Mary, University of London; and Chris Cannam
|
c@116
|
26 shall not be used in advertising or otherwise to promote the sale,
|
c@116
|
27 use or other dealings in this Software without prior written
|
c@116
|
28 authorization.
|
c@116
|
29 */
|
c@116
|
30
|
c@116
|
31 module icqt;
|
c@116
|
32
|
c@116
|
33 cqt = load cqt;
|
c@116
|
34 cm = load may.matrix.complex;
|
c@116
|
35 framer = load may.stream.framer;
|
c@116
|
36 win = load may.signal.window;
|
c@116
|
37 mm = load may.mathmisc;
|
c@116
|
38 vec = load may.vector;
|
c@116
|
39 resamp = load may.stream.resample;
|
c@116
|
40 manip = load may.stream.manipulate;
|
c@116
|
41 mat = load may.matrix;
|
c@116
|
42
|
c@116
|
43 icqt cq =
|
c@116
|
44 (kdata = cq.kernel;
|
c@116
|
45
|
c@116
|
46 // kdata.kernel is the kernel matrix for a single octave. It has
|
c@116
|
47 // width kdata.fftSize and height kdata.binsPerOctave *
|
c@116
|
48 // kdata.atomsPerFrame.
|
c@116
|
49 //
|
c@116
|
50 // kdata.fftHop is the overlap between kernel matrices in a single
|
c@116
|
51 // octave.
|
c@116
|
52 //
|
c@116
|
53 // cq.sampleRate is the output stream sample rate and cq.octaves
|
c@116
|
54 // the number of octaves.
|
c@116
|
55 //
|
c@116
|
56 // cq.cqComplex is the list of complex matrices containing the CQ
|
c@116
|
57 // output. Each has width kdata.atomsPerFrame * 2^(cq.octaves-1)
|
c@116
|
58 // and height kdata.binsPerOctave * cq.octaves.
|
c@116
|
59
|
c@116
|
60 bpo = kdata.binsPerOctave;
|
c@116
|
61 atomsPerFrame = kdata.atomsPerFrame;
|
c@116
|
62 octaves = cq.octaves;
|
c@116
|
63
|
c@116
|
64 // transform a single block, all octaves tall, into an array
|
c@116
|
65 // (indexed by octave) of lists of individual columns (valid
|
c@116
|
66 // values for that octave only)
|
c@116
|
67 decomposeOctaves mat = array
|
c@116
|
68 (map do oct:
|
c@116
|
69 inverted = cm.fromRows (reverse (cm.asRows mat));
|
c@116
|
70 octMat = cm.rowSlice inverted (bpo * oct) (bpo * (oct + 1));
|
c@116
|
71 gap = (mm.pow 2 oct) - 1;
|
c@116
|
72 pickFrom cols =
|
c@116
|
73 case cols of
|
c@116
|
74 c::cs: c::(pickFrom (drop gap cs));
|
c@116
|
75 _: [];
|
c@116
|
76 esac;
|
c@116
|
77 pickFrom (cm.asColumns octMat);
|
c@116
|
78 done [0..octaves-1]);
|
c@116
|
79
|
c@116
|
80 // transform a list of the arrays produced by decomposeOctaves
|
c@116
|
81 // into a single array (indexed by octave) of lists of the
|
c@116
|
82 // individual columns
|
c@116
|
83 flattenOctaves decomposed =
|
c@116
|
84 (flattenAux acc decomposed =
|
c@116
|
85 case decomposed of
|
c@116
|
86 chunk::rest:
|
c@116
|
87 flattenAux
|
c@116
|
88 (array
|
c@116
|
89 (map do oct:
|
c@116
|
90 acc[oct] ++ chunk[oct]
|
c@116
|
91 done [0..octaves-1]))
|
c@116
|
92 rest;
|
c@116
|
93 _: acc;
|
c@116
|
94 esac;
|
c@116
|
95 flattenAux (array (map \[] [0..octaves-1])) decomposed);
|
c@116
|
96
|
c@116
|
97 // given a list of columns, deinterleave and pile up each sequence
|
c@116
|
98 // of atomsPerFrame columns into a single tall column and return
|
c@116
|
99 // the resulting list of tall columns
|
c@116
|
100 pile cols =
|
c@116
|
101 (pileAux acc cols =
|
c@116
|
102 if (length cols) >= atomsPerFrame then
|
c@116
|
103 atoms = take atomsPerFrame cols;
|
c@116
|
104 juggled = concat (reverse (cm.asRows (cm.fromColumns atoms)));
|
c@116
|
105 pileAux (juggled :: acc) (drop atomsPerFrame cols);
|
c@116
|
106 else
|
c@116
|
107 reverse acc
|
c@116
|
108 fi;
|
c@116
|
109 pileAux [] cols);
|
c@116
|
110
|
c@116
|
111 octaveColumnLists =
|
c@116
|
112 map pile (list (flattenOctaves (map decomposeOctaves cq.cqComplex)));
|
c@116
|
113
|
c@116
|
114 for octaveColumnLists do l: println "octave column list length: \(length l)" done;
|
c@116
|
115
|
c@116
|
116 kernel = cm.transposed kdata.kernel; // right way around for the multiply
|
c@116
|
117
|
c@116
|
118 spectra =
|
c@116
|
119 map do l:
|
c@116
|
120 map do col:
|
c@116
|
121 res = cm.transposed (cm.product kernel (cm.newComplexColumnVector col));
|
c@116
|
122 cm.columnSlice res 0 (kdata.fftSize / 2 + 1)
|
c@116
|
123 done l;
|
c@116
|
124 done octaveColumnLists;
|
c@116
|
125
|
c@116
|
126 eprintln "calculated spectra, now to ifft, overlap-add...";
|
c@116
|
127
|
c@116
|
128 rates = map do oct: cq.sampleRate / (mm.pow 2 oct) done [0..cq.octaves-1];
|
c@116
|
129
|
c@116
|
130 resynthesised =
|
c@116
|
131 map2 do frames rate:
|
c@116
|
132 framer.complexStreamed rate kdata.fftSize
|
c@116
|
133 [ FrequencyDomain true, Window win.boxcar, Hop kdata.fftHop ]
|
c@116
|
134 frames;
|
c@116
|
135 done spectra rates;
|
c@116
|
136
|
c@116
|
137 eprintln "have streams:";
|
c@116
|
138 for resynthesised eprintln;
|
c@116
|
139
|
c@116
|
140 resampled = map (resamp.resampledTo cq.sampleRate) resynthesised;
|
c@116
|
141
|
c@116
|
142 manip.sum resampled;
|
c@116
|
143 );
|
c@116
|
144
|
c@116
|
145 {icqt}
|