src/lib/audio_filter.cc

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