DXR is a code search and navigation tool aimed at making sense of large projects. It supports full-text and regex searches as well as structural queries.

Header

Mercurial (c68fe15a81fc)

VCS Links

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
/*
 *  Copyright (c) 2015 The WebRTC project authors. All Rights Reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS.  All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */

// An implementation of a 3-band FIR filter-bank with DCT modulation, similar to
// the proposed in "Multirate Signal Processing for Communication Systems" by
// Fredric J Harris.
//
// The idea is to take a heterodyne system and change the order of the
// components to get something which is efficient to implement digitally.
//
// It is possible to separate the filter using the noble identity as follows:
//
// H(z) = H0(z^3) + z^-1 * H1(z^3) + z^-2 * H2(z^3)
//
// This is used in the analysis stage to first downsample serial to parallel
// and then filter each branch with one of these polyphase decompositions of the
// lowpass prototype. Because each filter is only a modulation of the prototype,
// it is enough to multiply each coefficient by the respective cosine value to
// shift it to the desired band. But because the cosine period is 12 samples,
// it requires separating the prototype even further using the noble identity.
// After filtering and modulating for each band, the output of all filters is
// accumulated to get the downsampled bands.
//
// A similar logic can be applied to the synthesis stage.

// MSVC++ requires this to be set before any other includes to get M_PI.
#define _USE_MATH_DEFINES

#include "modules/audio_processing/three_band_filter_bank.h"

#include <cmath>

#include "rtc_base/checks.h"

namespace webrtc {
namespace {

const size_t kNumBands = 3;
const size_t kSparsity = 4;

// Factors to take into account when choosing |kNumCoeffs|:
//   1. Higher |kNumCoeffs|, means faster transition, which ensures less
//      aliasing. This is especially important when there is non-linear
//      processing between the splitting and merging.
//   2. The delay that this filter bank introduces is
//      |kNumBands| * |kSparsity| * |kNumCoeffs| / 2, so it increases linearly
//      with |kNumCoeffs|.
//   3. The computation complexity also increases linearly with |kNumCoeffs|.
const size_t kNumCoeffs = 4;

// The Matlab code to generate these |kLowpassCoeffs| is:
//
// N = kNumBands * kSparsity * kNumCoeffs - 1;
// h = fir1(N, 1 / (2 * kNumBands), kaiser(N + 1, 3.5));
// reshape(h, kNumBands * kSparsity, kNumCoeffs);
//
// Because the total bandwidth of the lower and higher band is double the middle
// one (because of the spectrum parity), the low-pass prototype is half the
// bandwidth of 1 / (2 * |kNumBands|) and is then shifted with cosine modulation
// to the right places.
// A Kaiser window is used because of its flexibility and the alpha is set to
// 3.5, since that sets a stop band attenuation of 40dB ensuring a fast
// transition.
const float kLowpassCoeffs[kNumBands * kSparsity][kNumCoeffs] =
    {{-0.00047749f, -0.00496888f, +0.16547118f, +0.00425496f},
     {-0.00173287f, -0.01585778f, +0.14989004f, +0.00994113f},
     {-0.00304815f, -0.02536082f, +0.12154542f, +0.01157993f},
     {-0.00383509f, -0.02982767f, +0.08543175f, +0.00983212f},
     {-0.00346946f, -0.02587886f, +0.04760441f, +0.00607594f},
     {-0.00154717f, -0.01136076f, +0.01387458f, +0.00186353f},
     {+0.00186353f, +0.01387458f, -0.01136076f, -0.00154717f},
     {+0.00607594f, +0.04760441f, -0.02587886f, -0.00346946f},
     {+0.00983212f, +0.08543175f, -0.02982767f, -0.00383509f},
     {+0.01157993f, +0.12154542f, -0.02536082f, -0.00304815f},
     {+0.00994113f, +0.14989004f, -0.01585778f, -0.00173287f},
     {+0.00425496f, +0.16547118f, -0.00496888f, -0.00047749f}};

// Downsamples |in| into |out|, taking one every |kNumbands| starting from
// |offset|. |split_length| is the |out| length. |in| has to be at least
// |kNumBands| * |split_length| long.
void Downsample(const float* in,
                size_t split_length,
                size_t offset,
                float* out) {
  for (size_t i = 0; i < split_length; ++i) {
    out[i] = in[kNumBands * i + offset];
  }
}

// Upsamples |in| into |out|, scaling by |kNumBands| and accumulating it every
// |kNumBands| starting from |offset|. |split_length| is the |in| length. |out|
// has to be at least |kNumBands| * |split_length| long.
void Upsample(const float* in, size_t split_length, size_t offset, float* out) {
  for (size_t i = 0; i < split_length; ++i) {
    out[kNumBands * i + offset] += kNumBands * in[i];
  }
}

}  // namespace

// Because the low-pass filter prototype has half bandwidth it is possible to
// use a DCT to shift it in both directions at the same time, to the center
// frequencies [1 / 12, 3 / 12, 5 / 12].
ThreeBandFilterBank::ThreeBandFilterBank(size_t length)
    : in_buffer_(rtc::CheckedDivExact(length, kNumBands)),
      out_buffer_(in_buffer_.size()) {
  for (size_t i = 0; i < kSparsity; ++i) {
    for (size_t j = 0; j < kNumBands; ++j) {
      analysis_filters_.push_back(
          std::unique_ptr<SparseFIRFilter>(new SparseFIRFilter(
              kLowpassCoeffs[i * kNumBands + j], kNumCoeffs, kSparsity, i)));
      synthesis_filters_.push_back(
          std::unique_ptr<SparseFIRFilter>(new SparseFIRFilter(
              kLowpassCoeffs[i * kNumBands + j], kNumCoeffs, kSparsity, i)));
    }
  }
  dct_modulation_.resize(kNumBands * kSparsity);
  for (size_t i = 0; i < dct_modulation_.size(); ++i) {
    dct_modulation_[i].resize(kNumBands);
    for (size_t j = 0; j < kNumBands; ++j) {
      dct_modulation_[i][j] =
          2.f * cos(2.f * M_PI * i * (2.f * j + 1.f) / dct_modulation_.size());
    }
  }
}

ThreeBandFilterBank::~ThreeBandFilterBank() = default;

// The analysis can be separated in these steps:
//   1. Serial to parallel downsampling by a factor of |kNumBands|.
//   2. Filtering of |kSparsity| different delayed signals with polyphase
//      decomposition of the low-pass prototype filter and upsampled by a factor
//      of |kSparsity|.
//   3. Modulating with cosines and accumulating to get the desired band.
void ThreeBandFilterBank::Analysis(const float* in,
                                   size_t length,
                                   float* const* out) {
  RTC_CHECK_EQ(in_buffer_.size(), rtc::CheckedDivExact(length, kNumBands));
  for (size_t i = 0; i < kNumBands; ++i) {
    memset(out[i], 0, in_buffer_.size() * sizeof(*out[i]));
  }
  for (size_t i = 0; i < kNumBands; ++i) {
    Downsample(in, in_buffer_.size(), kNumBands - i - 1, &in_buffer_[0]);
    for (size_t j = 0; j < kSparsity; ++j) {
      const size_t offset = i + j * kNumBands;
      analysis_filters_[offset]->Filter(&in_buffer_[0],
                                        in_buffer_.size(),
                                        &out_buffer_[0]);
      DownModulate(&out_buffer_[0], out_buffer_.size(), offset, out);
    }
  }
}

// The synthesis can be separated in these steps:
//   1. Modulating with cosines.
//   2. Filtering each one with a polyphase decomposition of the low-pass
//      prototype filter upsampled by a factor of |kSparsity| and accumulating
//      |kSparsity| signals with different delays.
//   3. Parallel to serial upsampling by a factor of |kNumBands|.
void ThreeBandFilterBank::Synthesis(const float* const* in,
                                    size_t split_length,
                                    float* out) {
  RTC_CHECK_EQ(in_buffer_.size(), split_length);
  memset(out, 0, kNumBands * in_buffer_.size() * sizeof(*out));
  for (size_t i = 0; i < kNumBands; ++i) {
    for (size_t j = 0; j < kSparsity; ++j) {
      const size_t offset = i + j * kNumBands;
      UpModulate(in, in_buffer_.size(), offset, &in_buffer_[0]);
      synthesis_filters_[offset]->Filter(&in_buffer_[0],
                                         in_buffer_.size(),
                                         &out_buffer_[0]);
      Upsample(&out_buffer_[0], out_buffer_.size(), i, out);
    }
  }
}


// Modulates |in| by |dct_modulation_| and accumulates it in each of the
// |kNumBands| bands of |out|. |offset| is the index in the period of the
// cosines used for modulation. |split_length| is the length of |in| and each
// band of |out|.
void ThreeBandFilterBank::DownModulate(const float* in,
                                       size_t split_length,
                                       size_t offset,
                                       float* const* out) {
  for (size_t i = 0; i < kNumBands; ++i) {
    for (size_t j = 0; j < split_length; ++j) {
      out[i][j] += dct_modulation_[offset][i] * in[j];
    }
  }
}

// Modulates each of the |kNumBands| bands of |in| by |dct_modulation_| and
// accumulates them in |out|. |out| is cleared before starting to accumulate.
// |offset| is the index in the period of the cosines used for modulation.
// |split_length| is the length of each band of |in| and |out|.
void ThreeBandFilterBank::UpModulate(const float* const* in,
                                     size_t split_length,
                                     size_t offset,
                                     float* out) {
  memset(out, 0, split_length * sizeof(*out));
  for (size_t i = 0; i < kNumBands; ++i) {
    for (size_t j = 0; j < split_length; ++j) {
      out[j] += dct_modulation_[offset][i] * in[i][j];
    }
  }
}

}  // namespace webrtc