Chris@42
|
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
|
Chris@42
|
2
|
Chris@42
|
3 /*
|
Chris@42
|
4 Tipic
|
Chris@42
|
5
|
Chris@42
|
6 Centre for Digital Music, Queen Mary, University of London.
|
Chris@42
|
7
|
Chris@42
|
8 This program is free software; you can redistribute it and/or
|
Chris@42
|
9 modify it under the terms of the GNU General Public License as
|
Chris@42
|
10 published by the Free Software Foundation; either version 2 of the
|
Chris@42
|
11 License, or (at your option) any later version. See the file
|
Chris@42
|
12 COPYING included with this distribution for more information.
|
Chris@42
|
13 */
|
Chris@5
|
14
|
Chris@5
|
15 #include "PitchFilterbank.h"
|
Chris@5
|
16
|
Chris@42
|
17 #include "dsp/rateconversion/Resampler.h"
|
Chris@42
|
18 #include "dsp/signalconditioning/Filter.h"
|
Chris@5
|
19
|
Chris@5
|
20 #include "delays.h"
|
Chris@5
|
21 #include "filter-a.h"
|
Chris@5
|
22 #include "filter-b.h"
|
Chris@5
|
23
|
Chris@11
|
24 #include <deque>
|
Chris@11
|
25 #include <map>
|
Chris@6
|
26 #include <stdexcept>
|
Chris@11
|
27 #include <iostream>
|
Chris@59
|
28 #include <algorithm>
|
Chris@6
|
29
|
Chris@27
|
30 #include <cmath>
|
Chris@27
|
31
|
Chris@5
|
32 using namespace std;
|
Chris@5
|
33
|
Chris@5
|
34 static const int HIGHEST_FILTER_INDEX_AT_882 = 38; // MIDI pitch 59
|
Chris@5
|
35 static const int HIGHEST_FILTER_INDEX_AT_4410 = 74; // MIDI pitch 95
|
Chris@5
|
36 static const int HIGHEST_FILTER_INDEX_AT_22050 = 87; // MIDI pitch 108
|
Chris@5
|
37 static const int HIGHEST_FILTER_INDEX = HIGHEST_FILTER_INDEX_AT_22050;
|
Chris@5
|
38
|
Chris@5
|
39 class PitchFilterbank::D
|
Chris@5
|
40 {
|
Chris@5
|
41 public:
|
Chris@27
|
42 D(int sampleRate, double tuningFrequency) :
|
Chris@11
|
43 m_nfilters(HIGHEST_FILTER_INDEX + 1),
|
Chris@15
|
44 m_sampleRate(sampleRate),
|
Chris@28
|
45 m_tuningFrequency(tuningFrequency),
|
Chris@28
|
46 m_blockNo(0)
|
Chris@5
|
47 {
|
Chris@28
|
48 // To handle a non-440Hz tuning frequency, we resample the
|
Chris@28
|
49 // input and then adjust the output block timings
|
Chris@28
|
50 // accordingly. For a tuning freq >440 we want to lower the
|
Chris@28
|
51 // pitch of the input audio by slowing it down, therefore we
|
Chris@28
|
52 // want to pretend that it came in at a lower sample rate than
|
Chris@28
|
53 // it really did; for >440 the opposite applies. The effective
|
Chris@28
|
54 // input sample rate is the rate at which we pretend the audio
|
Chris@28
|
55 // was supplied. Rounding to the nearest int (because our
|
Chris@28
|
56 // resampler only supports integer rates) gives around 0.1Hz
|
Chris@28
|
57 // quantization close to 440Hz in 44.1kHz audio -- we could do
|
Chris@28
|
58 // better by using multiples of our source and target sample
|
Chris@28
|
59 // rates, but I think it probably isn't necessary.
|
Chris@29
|
60 m_effectiveInputRate =
|
Chris@28
|
61 int(round(m_sampleRate * (440.0 / m_tuningFrequency)));
|
Chris@27
|
62
|
Chris@18
|
63 //!!! nb "we use forward-backward filtering such that the
|
Chris@18
|
64 // resulting output signal has precisely zero phase distortion
|
Chris@18
|
65 // and a magnitude modified by the square of the filter’s
|
Chris@18
|
66 // magnitude response" -- we are not doing forward-backward
|
Chris@18
|
67 // here & so need to adapt magnitudes in compensation to match
|
Chris@18
|
68 // original
|
Chris@29
|
69
|
Chris@29
|
70 double snr = 50.0;
|
Chris@29
|
71 double bw = 0.05;
|
Chris@29
|
72 m_resamplers[882] = new Resampler(m_effectiveInputRate, 882, snr, bw);
|
Chris@29
|
73 m_resamplers[4410] = new Resampler(m_effectiveInputRate, 4410, snr, bw);
|
Chris@29
|
74 m_resamplers[22050] = new Resampler(m_effectiveInputRate, 22050, snr, bw);
|
Chris@11
|
75
|
Chris@11
|
76 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@5
|
77 int ix = i + 20;
|
Chris@5
|
78 int coeffs = sizeof(filter_a[0]) / sizeof(filter_a[0][0]);
|
Chris@5
|
79 vector<double> a(filter_a[ix], filter_a[ix] + coeffs);
|
Chris@5
|
80 vector<double> b(filter_b[ix], filter_b[ix] + coeffs);
|
Chris@10
|
81 m_filters.push_back(new Filter({ a, b }));
|
Chris@12
|
82 m_toCompensate.push_back(totalDelay(i));
|
Chris@5
|
83 }
|
Chris@11
|
84
|
Chris@11
|
85 m_filtered.resize(m_nfilters);
|
Chris@11
|
86 m_energies.resize(m_nfilters);
|
Chris@5
|
87 }
|
Chris@5
|
88
|
Chris@9
|
89 ~D() {
|
Chris@11
|
90 for (auto f: m_filters) delete f;
|
Chris@11
|
91 for (auto r: m_resamplers) delete r.second;
|
Chris@9
|
92 }
|
Chris@5
|
93
|
Chris@9
|
94 int getSampleRate() const { return m_sampleRate; }
|
Chris@27
|
95 double getTuningFrequency() const { return m_tuningFrequency; }
|
Chris@5
|
96
|
Chris@10
|
97 RealBlock process(const RealSequence &in) {
|
Chris@5
|
98
|
Chris@11
|
99 for (auto r: m_resamplers) {
|
Chris@11
|
100 m_resampled[r.first] = r.second->process(in.data(), in.size());
|
Chris@11
|
101 }
|
Chris@5
|
102
|
Chris@11
|
103 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@11
|
104 int rate = filterRate(i);
|
Chris@11
|
105 if (m_resampled.find(rate) == m_resampled.end()) {
|
Chris@11
|
106 throw logic_error("No resampled output for rate of filter");
|
Chris@10
|
107 }
|
Chris@13
|
108 pushFiltered(i, m_resampled[rate], false);
|
Chris@10
|
109 }
|
Chris@10
|
110
|
Chris@12
|
111 return energiesFromFiltered(false);
|
Chris@12
|
112 }
|
Chris@12
|
113
|
Chris@12
|
114 RealBlock getRemainingOutput() {
|
Chris@12
|
115
|
Chris@12
|
116 for (auto r: m_resamplers) {
|
Chris@12
|
117 RealSequence in(r.second->getLatency(), 0.0);
|
Chris@12
|
118 m_resampled[r.first] = r.second->process(in.data(), in.size());
|
Chris@12
|
119 }
|
Chris@12
|
120
|
Chris@12
|
121 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@12
|
122 int rate = filterRate(i);
|
Chris@13
|
123 pushFiltered(i, m_resampled[rate], true);
|
Chris@12
|
124 }
|
Chris@12
|
125
|
Chris@12
|
126 return energiesFromFiltered(true);
|
Chris@12
|
127 }
|
Chris@28
|
128
|
Chris@28
|
129 struct WindowPosition {
|
Chris@28
|
130 uint64_t start;
|
Chris@28
|
131 int size;
|
Chris@28
|
132 double factor;
|
Chris@28
|
133 };
|
Chris@28
|
134
|
Chris@28
|
135 WindowPosition windowPosition(int block, int i) {
|
Chris@28
|
136
|
Chris@32
|
137 int hop = 2205;
|
Chris@28
|
138
|
Chris@28
|
139 double rate = filterRate(i);
|
Chris@28
|
140 double topRate = 22050.0;
|
Chris@28
|
141 double rateRatio = topRate / rate;
|
Chris@29
|
142 double tuningRatio = m_sampleRate / double(m_effectiveInputRate);
|
Chris@28
|
143 double sizeRatio = tuningRatio / rateRatio;
|
Chris@28
|
144
|
Chris@32
|
145 uint64_t start(round((hop * uint64_t(block)) * sizeRatio));
|
Chris@28
|
146 int size(round((hop * 2) * sizeRatio));
|
Chris@28
|
147
|
Chris@28
|
148 return { start, size, rateRatio };
|
Chris@28
|
149 }
|
Chris@12
|
150
|
Chris@12
|
151 RealBlock energiesFromFiltered(bool drain) {
|
Chris@10
|
152
|
Chris@11
|
153 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@11
|
154
|
Chris@28
|
155 WindowPosition here = windowPosition(m_blockNo, i);
|
Chris@28
|
156 WindowPosition next = windowPosition(m_blockNo + 1, i);
|
Chris@28
|
157
|
Chris@28
|
158 int n = here.size;
|
Chris@28
|
159 int hop = next.start - here.start;
|
Chris@12
|
160
|
Chris@12
|
161 unsigned int minReq = n;
|
Chris@12
|
162 if (drain) minReq = hop;
|
Chris@12
|
163
|
Chris@31
|
164 // we use a separate buffer for each filter (not just one
|
Chris@31
|
165 // per resampling ratio) because each filter has a
|
Chris@31
|
166 // different delay, which we are compensating for when
|
Chris@31
|
167 // first filling the buffer.
|
Chris@31
|
168
|
Chris@12
|
169 while (m_filtered[i].size() >= minReq) {
|
Chris@28
|
170 double energy = calculateEnergy(m_filtered[i], n, here.factor);
|
Chris@11
|
171 m_energies[i].push_back(energy);
|
Chris@11
|
172 m_filtered[i] = RealSequence(m_filtered[i].begin() + hop,
|
Chris@11
|
173 m_filtered[i].end());
|
Chris@10
|
174 }
|
Chris@10
|
175 }
|
Chris@10
|
176
|
Chris@28
|
177 ++m_blockNo;
|
Chris@28
|
178
|
Chris@13
|
179 int minCols = 0, maxCols = 0;
|
Chris@11
|
180 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@11
|
181 int n = m_energies[i].size();
|
Chris@13
|
182 if (i == 0 || n < minCols) minCols = n;
|
Chris@13
|
183 if (i == 0 || n > maxCols) maxCols = n;
|
Chris@13
|
184 }
|
Chris@13
|
185
|
Chris@13
|
186 RealBlock out;
|
Chris@13
|
187
|
Chris@13
|
188 if (drain) {
|
Chris@13
|
189 out.resize(maxCols);
|
Chris@13
|
190 for (int j = 0; j < maxCols; ++j) {
|
Chris@13
|
191 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@13
|
192 if (m_energies[i].size() == 0) {
|
Chris@13
|
193 out[j].push_back(0.0);
|
Chris@13
|
194 } else {
|
Chris@13
|
195 out[j].push_back(m_energies[i][0]);
|
Chris@13
|
196 m_energies[i].pop_front();
|
Chris@13
|
197 }
|
Chris@13
|
198 }
|
Chris@11
|
199 }
|
Chris@13
|
200 } else {
|
Chris@13
|
201 out.resize(minCols);
|
Chris@13
|
202 for (int j = 0; j < minCols; ++j) {
|
Chris@13
|
203 for (int i = 0; i < m_nfilters; ++i) {
|
Chris@13
|
204 out[j].push_back(m_energies[i][0]);
|
Chris@13
|
205 m_energies[i].pop_front();
|
Chris@13
|
206 }
|
Chris@10
|
207 }
|
Chris@10
|
208 }
|
Chris@10
|
209
|
Chris@10
|
210 return out;
|
Chris@10
|
211 }
|
Chris@10
|
212
|
Chris@13
|
213 void pushFiltered(int i, const RealSequence &resampled, bool drain) {
|
Chris@12
|
214
|
Chris@13
|
215 RealSequence in(resampled);
|
Chris@13
|
216 if (drain) {
|
Chris@13
|
217 RealSequence pad(filterDelay(i), 0.0);
|
Chris@13
|
218 in.insert(in.end(), pad.begin(), pad.end());
|
Chris@13
|
219 }
|
Chris@13
|
220
|
Chris@13
|
221 int n = in.size();
|
Chris@10
|
222 RealSequence filtered(n, 0.0);
|
Chris@13
|
223 m_filters[i]->process(in.data(), filtered.data(), n);
|
Chris@12
|
224
|
Chris@12
|
225 int pushStart = 0, pushCount = n;
|
Chris@12
|
226
|
Chris@12
|
227 if (m_toCompensate[i] > 0) {
|
Chris@12
|
228 pushCount = max(pushCount - m_toCompensate[i], 0);
|
Chris@12
|
229 int compensating = n - pushCount;
|
Chris@12
|
230 pushStart += compensating;
|
Chris@12
|
231 m_toCompensate[i] -= compensating;
|
Chris@12
|
232 if (m_toCompensate[i] < 0) {
|
Chris@12
|
233 throw logic_error("Compensated for more latency than exists");
|
Chris@12
|
234 }
|
Chris@12
|
235 }
|
Chris@12
|
236
|
Chris@12
|
237 m_filtered[i].insert(m_filtered[i].end(),
|
Chris@12
|
238 filtered.begin() + pushStart,
|
Chris@12
|
239 filtered.begin() + pushStart + pushCount);
|
Chris@10
|
240 }
|
Chris@10
|
241
|
Chris@10
|
242 double calculateEnergy(const RealSequence &seq, int n, double factor) {
|
Chris@10
|
243 double energy = 0.0;
|
Chris@12
|
244 int sz = seq.size();
|
Chris@12
|
245 if (n > sz) n = sz;
|
Chris@10
|
246 for (int i = 0; i < n; ++i) {
|
Chris@11
|
247 energy += seq[i] * seq[i];
|
Chris@10
|
248 }
|
Chris@11
|
249 return energy * factor;
|
Chris@10
|
250 }
|
Chris@10
|
251
|
Chris@5
|
252 private:
|
Chris@11
|
253 int m_nfilters;
|
Chris@5
|
254 int m_sampleRate;
|
Chris@29
|
255 int m_effectiveInputRate;
|
Chris@27
|
256 double m_tuningFrequency;
|
Chris@5
|
257
|
Chris@5
|
258 // This vector is initialised with 88 filter instances.
|
Chris@5
|
259 // m_filters[n] (for n from 0 to 87) is for MIDI pitch 21+n, so we
|
Chris@5
|
260 // span the MIDI range from pitch 21 to 108. Then m_filters[n] has
|
Chris@5
|
261 // coefficients drawn from filter_a[20+n] and filter_b[20+n], and
|
Chris@5
|
262 // has effective delay filter_delay[20+n].
|
Chris@10
|
263 vector<Filter *> m_filters;
|
Chris@5
|
264
|
Chris@11
|
265 map<int, Resampler *> m_resamplers; // rate -> resampler
|
Chris@11
|
266 map<int, RealSequence> m_resampled;
|
Chris@12
|
267 vector<int> m_toCompensate; // latency remaining at start, per filter
|
Chris@11
|
268 vector<RealSequence> m_filtered;
|
Chris@11
|
269 vector<deque<double>> m_energies;
|
Chris@28
|
270 int m_blockNo;
|
Chris@16
|
271
|
Chris@11
|
272 Resampler *resamplerFor(int filterIndex) {
|
Chris@11
|
273 int rate = filterRate(filterIndex);
|
Chris@11
|
274 if (m_resamplers.find(rate) == m_resamplers.end()) {
|
Chris@11
|
275 throw std::logic_error("Filter index has unknown or unexpected rate");
|
Chris@11
|
276 }
|
Chris@11
|
277 return m_resamplers[rate];
|
Chris@11
|
278 }
|
Chris@11
|
279
|
Chris@11
|
280 int filterRate(int filterIndex) const {
|
Chris@11
|
281 if (filterIndex <= HIGHEST_FILTER_INDEX_AT_882) {
|
Chris@11
|
282 return 882;
|
Chris@11
|
283 } else if (filterIndex <= HIGHEST_FILTER_INDEX_AT_4410) {
|
Chris@11
|
284 return 4410;
|
Chris@11
|
285 } else {
|
Chris@11
|
286 return 22050;
|
Chris@11
|
287 }
|
Chris@5
|
288 }
|
Chris@12
|
289 int resamplerDelay(int filterIndex) const {
|
Chris@12
|
290 return const_cast<D *>(this)->resamplerFor(filterIndex)->getLatency();
|
Chris@12
|
291 }
|
Chris@6
|
292 int filterDelay(int filterIndex) const {
|
Chris@5
|
293 return filter_delay[20 + filterIndex];
|
Chris@5
|
294 }
|
Chris@6
|
295 int totalDelay(int filterIndex) const {
|
Chris@12
|
296 return resamplerDelay(filterIndex) + filterDelay(filterIndex);
|
Chris@5
|
297 }
|
Chris@5
|
298 };
|
Chris@5
|
299
|
Chris@27
|
300 PitchFilterbank::PitchFilterbank(int sampleRate, double tuningFrequency) :
|
Chris@15
|
301 m_d(new D(sampleRate, tuningFrequency))
|
Chris@9
|
302 {
|
Chris@9
|
303 }
|
Chris@9
|
304
|
Chris@9
|
305 PitchFilterbank::~PitchFilterbank()
|
Chris@9
|
306 {
|
Chris@9
|
307 delete m_d;
|
Chris@9
|
308 }
|
Chris@9
|
309
|
Chris@9
|
310 void
|
Chris@9
|
311 PitchFilterbank::reset()
|
Chris@9
|
312 {
|
Chris@9
|
313 int rate = m_d->getSampleRate();
|
Chris@27
|
314 double freq = m_d->getTuningFrequency();
|
Chris@9
|
315 delete m_d;
|
Chris@15
|
316 m_d = new D(rate, freq);
|
Chris@9
|
317 }
|
Chris@9
|
318
|
Chris@19
|
319 RealBlock
|
Chris@9
|
320 PitchFilterbank::process(const RealSequence &in)
|
Chris@9
|
321 {
|
Chris@9
|
322 return m_d->process(in);
|
Chris@9
|
323 }
|
Chris@9
|
324
|
Chris@19
|
325 RealBlock
|
Chris@9
|
326 PitchFilterbank::getRemainingOutput()
|
Chris@9
|
327 {
|
Chris@9
|
328 return m_d->getRemainingOutput();
|
Chris@9
|
329 }
|
Chris@9
|
330
|
Chris@27
|
331 void
|
Chris@27
|
332 PitchFilterbank::getPitchRange(int &minMidiPitch, int &maxMidiPitch)
|
Chris@27
|
333 {
|
Chris@27
|
334 minMidiPitch = 21;
|
Chris@27
|
335 maxMidiPitch = 108;
|
Chris@27
|
336 }
|
Chris@27
|
337
|
Chris@32
|
338 double
|
Chris@32
|
339 PitchFilterbank::getOutputSampleRate()
|
Chris@32
|
340 {
|
Chris@32
|
341 return 22050.0 / 2205.0;
|
Chris@32
|
342 }
|
Chris@27
|
343
|