Line data Source code
1 : /*
2 : * Copyright (c) 2015 The WebRTC project authors. All Rights Reserved.
3 : *
4 : * Use of this source code is governed by a BSD-style license
5 : * that can be found in the LICENSE file in the root of the source
6 : * tree. An additional intellectual property rights grant can be found
7 : * in the file PATENTS. All contributing project authors may
8 : * be found in the AUTHORS file in the root of the source tree.
9 : */
10 :
11 : // An implementation of a 3-band FIR filter-bank with DCT modulation, similar to
12 : // the proposed in "Multirate Signal Processing for Communication Systems" by
13 : // Fredric J Harris.
14 : //
15 : // The idea is to take a heterodyne system and change the order of the
16 : // components to get something which is efficient to implement digitally.
17 : //
18 : // It is possible to separate the filter using the noble identity as follows:
19 : //
20 : // H(z) = H0(z^3) + z^-1 * H1(z^3) + z^-2 * H2(z^3)
21 : //
22 : // This is used in the analysis stage to first downsample serial to parallel
23 : // and then filter each branch with one of these polyphase decompositions of the
24 : // lowpass prototype. Because each filter is only a modulation of the prototype,
25 : // it is enough to multiply each coefficient by the respective cosine value to
26 : // shift it to the desired band. But because the cosine period is 12 samples,
27 : // it requires separating the prototype even further using the noble identity.
28 : // After filtering and modulating for each band, the output of all filters is
29 : // accumulated to get the downsampled bands.
30 : //
31 : // A similar logic can be applied to the synthesis stage.
32 :
33 : // MSVC++ requires this to be set before any other includes to get M_PI.
34 : #ifndef _USE_MATH_DEFINES
35 : #define _USE_MATH_DEFINES
36 : #endif
37 :
38 : #include "webrtc/modules/audio_processing/three_band_filter_bank.h"
39 :
40 : #include <cmath>
41 :
42 : #include "webrtc/base/checks.h"
43 :
44 : namespace webrtc {
45 : namespace {
46 :
47 : const size_t kNumBands = 3;
48 : const size_t kSparsity = 4;
49 :
50 : // Factors to take into account when choosing |kNumCoeffs|:
51 : // 1. Higher |kNumCoeffs|, means faster transition, which ensures less
52 : // aliasing. This is especially important when there is non-linear
53 : // processing between the splitting and merging.
54 : // 2. The delay that this filter bank introduces is
55 : // |kNumBands| * |kSparsity| * |kNumCoeffs| / 2, so it increases linearly
56 : // with |kNumCoeffs|.
57 : // 3. The computation complexity also increases linearly with |kNumCoeffs|.
58 : const size_t kNumCoeffs = 4;
59 :
60 : // The Matlab code to generate these |kLowpassCoeffs| is:
61 : //
62 : // N = kNumBands * kSparsity * kNumCoeffs - 1;
63 : // h = fir1(N, 1 / (2 * kNumBands), kaiser(N + 1, 3.5));
64 : // reshape(h, kNumBands * kSparsity, kNumCoeffs);
65 : //
66 : // Because the total bandwidth of the lower and higher band is double the middle
67 : // one (because of the spectrum parity), the low-pass prototype is half the
68 : // bandwidth of 1 / (2 * |kNumBands|) and is then shifted with cosine modulation
69 : // to the right places.
70 : // A Kaiser window is used because of its flexibility and the alpha is set to
71 : // 3.5, since that sets a stop band attenuation of 40dB ensuring a fast
72 : // transition.
73 : const float kLowpassCoeffs[kNumBands * kSparsity][kNumCoeffs] =
74 : {{-0.00047749f, -0.00496888f, +0.16547118f, +0.00425496f},
75 : {-0.00173287f, -0.01585778f, +0.14989004f, +0.00994113f},
76 : {-0.00304815f, -0.02536082f, +0.12154542f, +0.01157993f},
77 : {-0.00383509f, -0.02982767f, +0.08543175f, +0.00983212f},
78 : {-0.00346946f, -0.02587886f, +0.04760441f, +0.00607594f},
79 : {-0.00154717f, -0.01136076f, +0.01387458f, +0.00186353f},
80 : {+0.00186353f, +0.01387458f, -0.01136076f, -0.00154717f},
81 : {+0.00607594f, +0.04760441f, -0.02587886f, -0.00346946f},
82 : {+0.00983212f, +0.08543175f, -0.02982767f, -0.00383509f},
83 : {+0.01157993f, +0.12154542f, -0.02536082f, -0.00304815f},
84 : {+0.00994113f, +0.14989004f, -0.01585778f, -0.00173287f},
85 : {+0.00425496f, +0.16547118f, -0.00496888f, -0.00047749f}};
86 :
87 : // Downsamples |in| into |out|, taking one every |kNumbands| starting from
88 : // |offset|. |split_length| is the |out| length. |in| has to be at least
89 : // |kNumBands| * |split_length| long.
90 0 : void Downsample(const float* in,
91 : size_t split_length,
92 : size_t offset,
93 : float* out) {
94 0 : for (size_t i = 0; i < split_length; ++i) {
95 0 : out[i] = in[kNumBands * i + offset];
96 : }
97 0 : }
98 :
99 : // Upsamples |in| into |out|, scaling by |kNumBands| and accumulating it every
100 : // |kNumBands| starting from |offset|. |split_length| is the |in| length. |out|
101 : // has to be at least |kNumBands| * |split_length| long.
102 0 : void Upsample(const float* in, size_t split_length, size_t offset, float* out) {
103 0 : for (size_t i = 0; i < split_length; ++i) {
104 0 : out[kNumBands * i + offset] += kNumBands * in[i];
105 : }
106 0 : }
107 :
108 : } // namespace
109 :
110 : // Because the low-pass filter prototype has half bandwidth it is possible to
111 : // use a DCT to shift it in both directions at the same time, to the center
112 : // frequencies [1 / 12, 3 / 12, 5 / 12].
113 0 : ThreeBandFilterBank::ThreeBandFilterBank(size_t length)
114 : : in_buffer_(rtc::CheckedDivExact(length, kNumBands)),
115 0 : out_buffer_(in_buffer_.size()) {
116 0 : for (size_t i = 0; i < kSparsity; ++i) {
117 0 : for (size_t j = 0; j < kNumBands; ++j) {
118 0 : analysis_filters_.push_back(
119 0 : std::unique_ptr<SparseFIRFilter>(new SparseFIRFilter(
120 0 : kLowpassCoeffs[i * kNumBands + j], kNumCoeffs, kSparsity, i)));
121 0 : synthesis_filters_.push_back(
122 0 : std::unique_ptr<SparseFIRFilter>(new SparseFIRFilter(
123 0 : kLowpassCoeffs[i * kNumBands + j], kNumCoeffs, kSparsity, i)));
124 : }
125 : }
126 0 : dct_modulation_.resize(kNumBands * kSparsity);
127 0 : for (size_t i = 0; i < dct_modulation_.size(); ++i) {
128 0 : dct_modulation_[i].resize(kNumBands);
129 0 : for (size_t j = 0; j < kNumBands; ++j) {
130 0 : dct_modulation_[i][j] =
131 0 : 2.f * cos(2.f * M_PI * i * (2.f * j + 1.f) / dct_modulation_.size());
132 : }
133 : }
134 0 : }
135 :
136 : ThreeBandFilterBank::~ThreeBandFilterBank() = default;
137 :
138 : // The analysis can be separated in these steps:
139 : // 1. Serial to parallel downsampling by a factor of |kNumBands|.
140 : // 2. Filtering of |kSparsity| different delayed signals with polyphase
141 : // decomposition of the low-pass prototype filter and upsampled by a factor
142 : // of |kSparsity|.
143 : // 3. Modulating with cosines and accumulating to get the desired band.
144 0 : void ThreeBandFilterBank::Analysis(const float* in,
145 : size_t length,
146 : float* const* out) {
147 0 : RTC_CHECK_EQ(in_buffer_.size(), rtc::CheckedDivExact(length, kNumBands));
148 0 : for (size_t i = 0; i < kNumBands; ++i) {
149 0 : memset(out[i], 0, in_buffer_.size() * sizeof(*out[i]));
150 : }
151 0 : for (size_t i = 0; i < kNumBands; ++i) {
152 0 : Downsample(in, in_buffer_.size(), kNumBands - i - 1, &in_buffer_[0]);
153 0 : for (size_t j = 0; j < kSparsity; ++j) {
154 0 : const size_t offset = i + j * kNumBands;
155 0 : analysis_filters_[offset]->Filter(&in_buffer_[0],
156 : in_buffer_.size(),
157 0 : &out_buffer_[0]);
158 0 : DownModulate(&out_buffer_[0], out_buffer_.size(), offset, out);
159 : }
160 : }
161 0 : }
162 :
163 : // The synthesis can be separated in these steps:
164 : // 1. Modulating with cosines.
165 : // 2. Filtering each one with a polyphase decomposition of the low-pass
166 : // prototype filter upsampled by a factor of |kSparsity| and accumulating
167 : // |kSparsity| signals with different delays.
168 : // 3. Parallel to serial upsampling by a factor of |kNumBands|.
169 0 : void ThreeBandFilterBank::Synthesis(const float* const* in,
170 : size_t split_length,
171 : float* out) {
172 0 : RTC_CHECK_EQ(in_buffer_.size(), split_length);
173 0 : memset(out, 0, kNumBands * in_buffer_.size() * sizeof(*out));
174 0 : for (size_t i = 0; i < kNumBands; ++i) {
175 0 : for (size_t j = 0; j < kSparsity; ++j) {
176 0 : const size_t offset = i + j * kNumBands;
177 0 : UpModulate(in, in_buffer_.size(), offset, &in_buffer_[0]);
178 0 : synthesis_filters_[offset]->Filter(&in_buffer_[0],
179 : in_buffer_.size(),
180 0 : &out_buffer_[0]);
181 0 : Upsample(&out_buffer_[0], out_buffer_.size(), i, out);
182 : }
183 : }
184 0 : }
185 :
186 :
187 : // Modulates |in| by |dct_modulation_| and accumulates it in each of the
188 : // |kNumBands| bands of |out|. |offset| is the index in the period of the
189 : // cosines used for modulation. |split_length| is the length of |in| and each
190 : // band of |out|.
191 0 : void ThreeBandFilterBank::DownModulate(const float* in,
192 : size_t split_length,
193 : size_t offset,
194 : float* const* out) {
195 0 : for (size_t i = 0; i < kNumBands; ++i) {
196 0 : for (size_t j = 0; j < split_length; ++j) {
197 0 : out[i][j] += dct_modulation_[offset][i] * in[j];
198 : }
199 : }
200 0 : }
201 :
202 : // Modulates each of the |kNumBands| bands of |in| by |dct_modulation_| and
203 : // accumulates them in |out|. |out| is cleared before starting to accumulate.
204 : // |offset| is the index in the period of the cosines used for modulation.
205 : // |split_length| is the length of each band of |in| and |out|.
206 0 : void ThreeBandFilterBank::UpModulate(const float* const* in,
207 : size_t split_length,
208 : size_t offset,
209 : float* out) {
210 0 : memset(out, 0, split_length * sizeof(*out));
211 0 : for (size_t i = 0; i < kNumBands; ++i) {
212 0 : for (size_t j = 0; j < split_length; ++j) {
213 0 : out[j] += dct_modulation_[offset][i] * in[i][j];
214 : }
215 : }
216 0 : }
217 :
218 : } // namespace webrtc
|
