|
c@1
|
1
|
|
c@1
|
2 module cqtkernel;
|
|
c@1
|
3
|
|
c@3
|
4 vec = load may.vector;
|
|
c@3
|
5 bf = load may.vector.blockfuncs;
|
|
c@3
|
6 complex = load may.complex;
|
|
c@3
|
7 window = load may.signal.window;
|
|
c@3
|
8 fft = load may.transform.fft;
|
|
c@4
|
9 pl = load may.plot;
|
|
c@3
|
10
|
|
c@2
|
11 { pow, round, floor, ceil, nextPowerOfTwo } = load may.mathmisc;
|
|
c@1
|
12
|
|
c@1
|
13 fs = 48000;
|
|
c@1
|
14
|
|
c@1
|
15 fmax = fs/2;
|
|
c@1
|
16
|
|
c@1
|
17 bins = 24;
|
|
c@1
|
18
|
|
c@1
|
19 q = 1;
|
|
c@1
|
20
|
|
c@1
|
21 atomHopFactor = 0.25;
|
|
c@1
|
22
|
|
c@1
|
23 thresh = 0.0005;
|
|
c@1
|
24
|
|
c@1
|
25 fmin = (fmax/2) * (pow 2 (1/bins));
|
|
c@1
|
26
|
|
c@1
|
27 bigQ = q / ((pow 2 (1/bins)) - 1);
|
|
c@1
|
28
|
|
c@1
|
29 nk_max = round(bigQ * fs / fmin);
|
|
c@1
|
30
|
|
c@1
|
31 nk_min = round(bigQ * fs / (fmin * (pow 2 ((bins-1)/bins))));
|
|
c@1
|
32
|
|
c@1
|
33 atomHop = round(nk_min * atomHopFactor);
|
|
c@1
|
34
|
|
c@1
|
35 first_center = atomHop * Math#ceil(Math#ceil(nk_max/2) / atomHop);
|
|
c@1
|
36
|
|
c@1
|
37 fftLen = nextPowerOfTwo (first_center + Math#ceil(nk_max/2));
|
|
c@1
|
38
|
|
c@1
|
39 println "fs = \(fs), fmax = \(fmax), bins = \(bins), q = \(q), atomHopFactor = \(atomHopFactor), thresh = \(thresh)";
|
|
c@1
|
40
|
|
c@1
|
41 println "fmin = \(fmin), bigQ = \(bigQ), nk_max = \(nk_max), nk_min = \(nk_min), atomHop = \(atomHop), first_center = \(first_center), fftLen = \(fftLen)";
|
|
c@1
|
42
|
|
c@2
|
43 winNr = floor((fftLen - ceil(nk_max/2) - first_center) / atomHop) + 1;
|
|
c@2
|
44
|
|
c@2
|
45 last_center = first_center + (winNr - 1) * atomHop;
|
|
c@2
|
46
|
|
c@2
|
47 fftHop = (last_center + atomHop) - first_center;
|
|
c@2
|
48
|
|
c@2
|
49 fftOverlap = ((fftLen - fftHop) / fftLen) * 100; // as % -- why? just for diagnostics?
|
|
c@2
|
50
|
|
c@2
|
51 println "winNr = \(winNr), last_center = \(last_center), fftHop = \(fftHop), fftOverlap = \(fftOverlap)%";
|
|
c@2
|
52
|
|
c@3
|
53 fftFunc = fft.forward fftLen;
|
|
c@2
|
54
|
|
c@3
|
55 // Note the MATLAB uses exp(2*pi*1i*x) for a complex generating
|
|
c@3
|
56 // function. We can't do that here; we need to generate real and imag
|
|
c@3
|
57 // parts separately as real = cos(2*pi*x), imag = sin(2*pi*x).
|
|
c@2
|
58
|
|
c@3
|
59 kernels = map do k:
|
|
c@3
|
60
|
|
c@3
|
61 nk = round(bigQ * fs / (fmin * (pow 2 ((k-1)/bins))));
|
|
c@3
|
62 println "k = \(k) -> nk = \(nk)";
|
|
c@3
|
63
|
|
c@3
|
64 win = bf.divideBy nk (bf.sqrt (window.blackmanHarris nk));
|
|
c@3
|
65
|
|
c@3
|
66 fk = fmin * (pow 2 ((k-1)/bins));
|
|
c@3
|
67
|
|
c@3
|
68 println "fk = \(fk)";
|
|
c@3
|
69
|
|
c@3
|
70 genKernel f = bf.multiply win
|
|
c@3
|
71 (vec.fromList (map do i: f (2 * pi * fk * i / fs) done [0..nk-1]));
|
|
c@3
|
72
|
|
c@3
|
73 reals = genKernel cos;
|
|
c@3
|
74 imags = genKernel sin;
|
|
c@3
|
75
|
|
c@3
|
76 atomOffset = first_center - ceil(nk/2);
|
|
c@3
|
77
|
|
c@3
|
78 map do i:
|
|
c@3
|
79
|
|
c@3
|
80 shift = atomOffset + ((i-1) * atomHop);
|
|
c@3
|
81
|
|
c@4
|
82 println "shift = \(shift)";
|
|
c@4
|
83
|
|
c@3
|
84 specKernel = fftFunc
|
|
c@3
|
85 (complex.complexArray
|
|
c@3
|
86 (vec.concat [vec.zeros shift, reals])
|
|
c@3
|
87 (vec.concat [vec.zeros shift, imags]));
|
|
c@3
|
88
|
|
c@3
|
89 map do c:
|
|
c@4
|
90 if complex.magnitude c < thresh then complex.zero
|
|
c@4
|
91 else complex.scale (1/fftLen) c fi
|
|
c@3
|
92 done specKernel;
|
|
c@4
|
93
|
|
c@3
|
94 done [1..winNr];
|
|
c@3
|
95
|
|
c@3
|
96 done [1..bins];
|
|
c@3
|
97
|
|
c@4
|
98 array (map array (concat kernels));
|
|
c@2
|
99
|
|
c@1
|
100
|
