src/lib/audio_filter.cc

   1 /*
   2     Copyright (C) 2014-2021 Carl Hetherington <cth@carlh.net>
   3
   4     This file is part of DCP-o-matic.
   5
   6     DCP-o-matic is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     DCP-o-matic is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with DCP-o-matic.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 */
  20
  21
  22 #include "audio_filter.h"
  23 #include "audio_buffers.h"
  24 #include "util.h"
  25 #include <cmath>
  26
  27
  28 using std::make_shared;
  29 using std::min;
  30 using std::shared_ptr;
  31
  32
  33 /** @return array of floats which the caller must destroy with delete[] */
  34 float *
  35 AudioFilter::sinc_blackman (float cutoff, bool invert) const
  36 {
  37         float* ir = new float[_M + 1];
  38
  39         /* Impulse response */
  40
  41         for (int i = 0; i <= _M; ++i) {
  42                 if (i == (_M / 2)) {
  43                         ir[i] = 2 * M_PI * cutoff;
  44                 } else {
  45                         /* sinc */
  46                         ir[i] = sin (2 * M_PI * cutoff * (i - _M / 2)) / (i - _M / 2);
  47                         /* Blackman window */
  48                         ir[i] *= (0.42 - 0.5 * cos (2 * M_PI * i / _M) + 0.08 * cos (4 * M_PI * i / _M));
  49                 }
  50         }
  51
  52         /* Normalise */
  53
  54         float sum = 0;
  55         for (int i = 0; i <= _M; ++i) {
  56                 sum += ir[i];
  57         }
  58
  59         for (int i = 0; i <= _M; ++i) {
  60                 ir[i] /= sum;
  61         }
  62
  63         /* Frequency inversion (swapping low-pass for high-pass, or whatever) */
  64
  65         if (invert) {
  66                 for (int i = 0; i <= _M; ++i) {
  67                         ir[i] = -ir[i];
  68                 }
  69                 ir[_M / 2] += 1;
  70         }
  71
  72         return ir;
  73 }
  74
  75
  76 AudioFilter::~AudioFilter ()
  77 {
  78         delete[] _ir;
  79 }
  80
  81
  82 shared_ptr<AudioBuffers>
  83 AudioFilter::run (shared_ptr<const AudioBuffers> in)
  84 {
  85         auto out = make_shared<AudioBuffers>(in->channels(), in->frames());
  86
  87         if (!_tail) {
  88                 _tail.reset (new AudioBuffers (in->channels(), _M + 1));
  89                 _tail->make_silent ();
  90         }
  91
  92         int const channels = in->channels ();
  93         int const frames = in->frames ();
  94
  95         for (int i = 0; i < channels; ++i) {
  96                 auto tail_p = _tail->data (i);
  97                 auto in_p = in->data (i);
  98                 auto out_p = out->data (i);
  99                 for (int j = 0; j < frames; ++j) {
 100                         float s = 0;
 101                         for (int k = 0; k <= _M; ++k) {
 102                                 if ((j - k) < 0) {
 103                                         s += tail_p[j - k + _M + 1] * _ir[k];
 104                                 } else {
 105                                         s += in_p[j - k] * _ir[k];
 106                                 }
 107                         }
 108
 109                         out_p[j] = s;
 110                 }
 111         }
 112
 113         int const amount = min (in->frames(), _tail->frames());
 114         if (amount < _tail->frames ()) {
 115                 _tail->move (_tail->frames() - amount, amount, 0);
 116         }
 117         _tail->copy_from (in.get(), amount, in->frames() - amount, _tail->frames () - amount);
 118
 119         return out;
 120 }
 121
 122
 123 void
 124 AudioFilter::flush ()
 125 {
 126         _tail.reset ();
 127 }
 128
 129
 130 LowPassAudioFilter::LowPassAudioFilter (float transition_bandwidth, float cutoff)
 131         : AudioFilter (transition_bandwidth)
 132 {
 133         _ir = sinc_blackman (cutoff, false);
 134 }
 135
 136
 137 HighPassAudioFilter::HighPassAudioFilter (float transition_bandwidth, float cutoff)
 138         : AudioFilter (transition_bandwidth)
 139 {
 140         _ir = sinc_blackman (cutoff, true);
 141 }
 142
 143
 144 BandPassAudioFilter::BandPassAudioFilter (float transition_bandwidth, float lower, float higher)
 145         : AudioFilter (transition_bandwidth)
 146 {
 147         auto lpf = sinc_blackman (lower, false);
 148         auto hpf = sinc_blackman (higher, true);
 149
 150         delete[] _ir;
 151         _ir = new float[_M + 1];
 152         for (int i = 0; i <= _M; ++i) {
 153                 _ir[i] = lpf[i] + hpf[i];
 154         }
 155
 156         delete[] lpf;
 157         delete[] hpf;
 158
 159         /* We now have a band-stop, so invert for band-pass */
 160         for (int i = 0; i <= _M; ++i) {
 161                 _ir[i] = -_ir[i];
 162         }
 163
 164         _ir[_M / 2] += 1;
 165 }