annotate src/PitchFilterbank.cpp @ 60:1ea2aed23d4a tip

Fix version
author Chris Cannam
date Thu, 13 Feb 2020 13:37:36 +0000
parents 14823d51a573
children
rev   line source
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