src/lib/audio_filter.cc

   1 /*
   2     Copyright (C) 2014 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 #include "audio_filter.h"
  22 #include "audio_buffers.h"
  23 #include <cmath>
  24
  25 using std::min;
  26 using boost::shared_ptr;
  27
  28 /** @return array of floats which the caller must destroy with delete[] */
  29 float *
  30 AudioFilter::sinc_blackman (float cutoff, bool invert) const
  31 {
  32         float* ir = new float[_M + 1];
  33
  34         /* Impulse response */
  35
  36         for (int i = 0; i <= _M; ++i) {
  37                 if (i == (_M / 2)) {
  38                         ir[i] = 2 * M_PI * cutoff;
  39                 } else {
  40                         /* sinc */
  41                         ir[i] = sin (2 * M_PI * cutoff * (i - _M / 2)) / (i - _M / 2);
  42                         /* Blackman window */
  43                         ir[i] *= (0.42 - 0.5 * cos (2 * M_PI * i / _M) + 0.08 * cos (4 * M_PI * i / _M));
  44                 }
  45         }
  46
  47         /* Normalise */
  48
  49         float sum = 0;
  50         for (int i = 0; i <= _M; ++i) {
  51                 sum += ir[i];
  52         }
  53
  54         for (int i = 0; i <= _M; ++i) {
  55                 ir[i] /= sum;
  56         }
  57
  58         /* Frequency inversion (swapping low-pass for high-pass, or whatever) */
  59
  60         if (invert) {
  61                 for (int i = 0; i <= _M; ++i) {
  62                         ir[i] = -ir[i];
  63                 }
  64                 ir[_M / 2] += 1;
  65         }
  66
  67         return ir;
  68 }
  69
  70 AudioFilter::~AudioFilter ()
  71 {
  72         delete[] _ir;
  73 }
  74
  75 shared_ptr<AudioBuffers>
  76 AudioFilter::run (shared_ptr<const AudioBuffers> in)
  77 {
  78         shared_ptr<AudioBuffers> out (new AudioBuffers (in->channels(), in->frames()));
  79
  80         if (!_tail) {
  81                 _tail.reset (new AudioBuffers (in->channels(), _M + 1));
  82                 _tail->make_silent ();
  83         }
  84
  85         int const channels = in->channels ();
  86         int const frames = in->frames ();
  87
  88         for (int i = 0; i < channels; ++i) {
  89                 float* tail_p = _tail->data (i);
  90                 float* in_p = in->data (i);
  91                 float* out_p = out->data (i);
  92                 for (int j = 0; j < frames; ++j) {
  93                         float s = 0;
  94                         for (int k = 0; k <= _M; ++k) {
  95                                 if ((j - k) < 0) {
  96                                         s += tail_p[j - k + _M + 1] * _ir[k];
  97                                 } else {
  98                                         s += in_p[j - k] * _ir[k];
  99                                 }
 100                         }
 101
 102                         out_p[j] = s;
 103                 }
 104         }
 105
 106         int const amount = min (in->frames(), _tail->frames());
 107         if (amount < _tail->frames ()) {
 108                 _tail->move (amount, 0, _tail->frames() - amount);
 109         }
 110         _tail->copy_from (in.get(), amount, in->frames() - amount, _tail->frames () - amount);
 111
 112         return out;
 113 }
 114
 115 void
 116 AudioFilter::flush ()
 117 {
 118         _tail.reset ();
 119 }
 120
 121 LowPassAudioFilter::LowPassAudioFilter (float transition_bandwidth, float cutoff)
 122         : AudioFilter (transition_bandwidth)
 123 {
 124         _ir = sinc_blackman (cutoff, false);
 125 }
 126
 127
 128 HighPassAudioFilter::HighPassAudioFilter (float transition_bandwidth, float cutoff)
 129         : AudioFilter (transition_bandwidth)
 130 {
 131         _ir = sinc_blackman (cutoff, true);
 132 }
 133
 134 BandPassAudioFilter::BandPassAudioFilter (float transition_bandwidth, float lower, float higher)
 135         : AudioFilter (transition_bandwidth)
 136 {
 137         float* lpf = sinc_blackman (lower, false);
 138         float* hpf = sinc_blackman (higher, true);
 139
 140         delete[] _ir;
 141         _ir = new float[_M + 1];
 142         for (int i = 0; i <= _M; ++i) {
 143                 _ir[i] = lpf[i] + hpf[i];
 144         }
 145
 146         delete[] lpf;
 147         delete[] hpf;
 148
 149         /* We now have a band-stop, so invert for band-pass */
 150         for (int i = 0; i <= _M; ++i) {
 151                 _ir[i] = -_ir[i];
 152         }
 153
 154         _ir[_M / 2] += 1;
 155 }