Merge remote-tracking branch 'origin/master' into 2.0
[dcpomatic.git] / 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 <cmath>
21 #include "audio_filter.h"
22 #include "audio_buffers.h"
23
24 using std::vector;
25 using std::min;
26 using boost::shared_ptr;
27
28 vector<float>
29 AudioFilter::sinc_blackman (float cutoff, bool invert) const
30 {
31         vector<float> ir (_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 shared_ptr<AudioBuffers>
70 AudioFilter::run (shared_ptr<AudioBuffers> in)
71 {
72         shared_ptr<AudioBuffers> out (new AudioBuffers (in->channels(), in->frames()));
73         
74         if (!_tail) {
75                 _tail.reset (new AudioBuffers (in->channels(), _M + 1));
76                 _tail->make_silent ();
77         }
78         
79         for (int i = 0; i < in->channels(); ++i) {
80                 for (int j = 0; j < in->frames(); ++j) {
81                         float s = 0;
82                         for (int k = 0; k <= _M; ++k) {
83                                 if ((j - k) < 0) {
84                                         s += _tail->data(i)[j - k + _M + 1] * _ir[k];
85                                 } else {
86                                         s += in->data(i)[j - k] * _ir[k];
87                                 }
88                         }
89                         
90                         out->data(i)[j] = s;
91                 }
92         }
93         
94         int const amount = min (in->frames(), _tail->frames());
95         if (amount < _tail->frames ()) {
96                 _tail->move (amount, 0, _tail->frames() - amount);
97         }
98         _tail->copy_from (in.get(), amount, in->frames() - amount, _tail->frames () - amount);
99         
100         return out;
101 }
102
103 void
104 AudioFilter::flush ()
105 {
106         _tail.reset ();
107 }
108
109 LowPassAudioFilter::LowPassAudioFilter (float transition_bandwidth, float cutoff)
110         : AudioFilter (transition_bandwidth)
111 {
112         _ir = sinc_blackman (cutoff, false);
113 }
114
115
116 HighPassAudioFilter::HighPassAudioFilter (float transition_bandwidth, float cutoff)
117         : AudioFilter (transition_bandwidth)
118 {
119         _ir = sinc_blackman (cutoff, true);
120 }
121
122 BandPassAudioFilter::BandPassAudioFilter (float transition_bandwidth, float lower, float higher)
123         : AudioFilter (transition_bandwidth)
124 {
125         vector<float> lpf = sinc_blackman (lower, false);
126         vector<float> hpf = sinc_blackman (higher, true);
127         
128         _ir.resize (_M + 1);
129         for (int i = 0; i <= _M; ++i) {
130                 _ir[i] = lpf[i] + hpf[i];
131         }
132         
133         /* We now have a band-stop, so invert for band-pass */
134         for (int i = 0; i <= _M; ++i) {
135                 _ir[i] = -_ir[i];
136         }
137         
138         _ir[_M / 2] += 1;
139 }