mel-computations.cc (15064B)
1 /** 2 * Copyright (c) 2022 Xiaomi Corporation (authors: Fangjun Kuang) 3 * 4 * See LICENSE for clarification regarding multiple authors 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 // This file is copied/modified from kaldi/src/feat/mel-computations.cc 20 21 #include "kaldi-native-fbank/csrc/mel-computations.h" 22 23 #include <stdio.h> 24 25 #include <algorithm> 26 #include <sstream> 27 #include <vector> 28 29 #include "kaldi-native-fbank/csrc/feature-window.h" 30 #include "kaldi-native-fbank/csrc/kaldi-math.h" 31 #include "kaldi-native-fbank/csrc/log.h" 32 33 namespace knf { 34 35 std::ostream &operator<<(std::ostream &os, const MelBanksOptions &opts) { 36 os << opts.ToString(); 37 return os; 38 } 39 40 float MelBanks::VtlnWarpFreq( 41 float vtln_low_cutoff, // upper+lower frequency cutoffs for VTLN. 42 float vtln_high_cutoff, 43 float low_freq, // upper+lower frequency cutoffs in mel computation 44 float high_freq, float vtln_warp_factor, float freq) { 45 /// This computes a VTLN warping function that is not the same as HTK's one, 46 /// but has similar inputs (this function has the advantage of never producing 47 /// empty bins). 48 49 /// This function computes a warp function F(freq), defined between low_freq 50 /// and high_freq inclusive, with the following properties: 51 /// F(low_freq) == low_freq 52 /// F(high_freq) == high_freq 53 /// The function is continuous and piecewise linear with two inflection 54 /// points. 55 /// The lower inflection point (measured in terms of the unwarped 56 /// frequency) is at frequency l, determined as described below. 57 /// The higher inflection point is at a frequency h, determined as 58 /// described below. 59 /// If l <= f <= h, then F(f) = f/vtln_warp_factor. 60 /// If the higher inflection point (measured in terms of the unwarped 61 /// frequency) is at h, then max(h, F(h)) == vtln_high_cutoff. 62 /// Since (by the last point) F(h) == h/vtln_warp_factor, then 63 /// max(h, h/vtln_warp_factor) == vtln_high_cutoff, so 64 /// h = vtln_high_cutoff / max(1, 1/vtln_warp_factor). 65 /// = vtln_high_cutoff * min(1, vtln_warp_factor). 66 /// If the lower inflection point (measured in terms of the unwarped 67 /// frequency) is at l, then min(l, F(l)) == vtln_low_cutoff 68 /// This implies that l = vtln_low_cutoff / min(1, 1/vtln_warp_factor) 69 /// = vtln_low_cutoff * max(1, vtln_warp_factor) 70 71 if (freq < low_freq || freq > high_freq) 72 return freq; // in case this gets called 73 // for out-of-range frequencies, just return the freq. 74 75 KNF_CHECK_GT(vtln_low_cutoff, low_freq); 76 KNF_CHECK_LT(vtln_high_cutoff, high_freq); 77 78 float one = 1.0f; 79 float l = vtln_low_cutoff * std::max(one, vtln_warp_factor); 80 float h = vtln_high_cutoff * std::min(one, vtln_warp_factor); 81 float scale = 1.0f / vtln_warp_factor; 82 float Fl = scale * l; // F(l); 83 float Fh = scale * h; // F(h); 84 KNF_CHECK(l > low_freq && h < high_freq); 85 // slope of left part of the 3-piece linear function 86 float scale_left = (Fl - low_freq) / (l - low_freq); 87 // [slope of center part is just "scale"] 88 89 // slope of right part of the 3-piece linear function 90 float scale_right = (high_freq - Fh) / (high_freq - h); 91 92 if (freq < l) { 93 return low_freq + scale_left * (freq - low_freq); 94 } else if (freq < h) { 95 return scale * freq; 96 } else { // freq >= h 97 return high_freq + scale_right * (freq - high_freq); 98 } 99 } 100 101 float MelBanks::VtlnWarpMelFreq( 102 float vtln_low_cutoff, // upper+lower frequency cutoffs for VTLN. 103 float vtln_high_cutoff, 104 float low_freq, // upper+lower frequency cutoffs in mel computation 105 float high_freq, float vtln_warp_factor, float mel_freq) { 106 return MelScale(VtlnWarpFreq(vtln_low_cutoff, vtln_high_cutoff, low_freq, 107 high_freq, vtln_warp_factor, 108 InverseMelScale(mel_freq))); 109 } 110 111 MelBanks::MelBanks(const MelBanksOptions &opts, 112 const FrameExtractionOptions &frame_opts, 113 float vtln_warp_factor) { 114 if (opts.is_librosa) { 115 InitLibrosaMelBanks(opts, frame_opts, vtln_warp_factor); 116 } else { 117 InitKaldiMelBanks(opts, frame_opts, vtln_warp_factor); 118 } 119 } 120 121 void MelBanks::InitKaldiMelBanks(const MelBanksOptions &opts, 122 const FrameExtractionOptions &frame_opts, 123 float vtln_warp_factor) { 124 htk_mode_ = opts.htk_mode; 125 int32_t num_bins = opts.num_bins; 126 if (num_bins < 3) { 127 KNF_LOG(FATAL) << "Must have at least 3 mel bins"; 128 } 129 130 float sample_freq = frame_opts.samp_freq; 131 int32_t window_length_padded = frame_opts.PaddedWindowSize(); 132 KNF_CHECK_EQ(window_length_padded % 2, 0); 133 134 int32_t num_fft_bins = window_length_padded / 2; 135 float nyquist = 0.5f * sample_freq; 136 137 float low_freq = opts.low_freq, high_freq; 138 if (opts.high_freq > 0.0f) { 139 high_freq = opts.high_freq; 140 } else { 141 high_freq = nyquist + opts.high_freq; 142 } 143 144 if (low_freq < 0.0f || low_freq >= nyquist || high_freq <= 0.0f || 145 high_freq > nyquist || high_freq <= low_freq) { 146 KNF_LOG(FATAL) << "Bad values in options: low-freq " << low_freq 147 << " and high-freq " << high_freq << " vs. nyquist " 148 << nyquist; 149 } 150 151 float fft_bin_width = sample_freq / window_length_padded; 152 // fft-bin width [think of it as Nyquist-freq / half-window-length] 153 154 float mel_low_freq = MelScale(low_freq); 155 float mel_high_freq = MelScale(high_freq); 156 157 debug_ = opts.debug_mel; 158 159 // divide by num_bins+1 in next line because of end-effects where the bins 160 // spread out to the sides. 161 float mel_freq_delta = (mel_high_freq - mel_low_freq) / (num_bins + 1); 162 163 float vtln_low = opts.vtln_low, vtln_high = opts.vtln_high; 164 if (vtln_high < 0.0f) { 165 vtln_high += nyquist; 166 } 167 168 if (vtln_warp_factor != 1.0f && 169 (vtln_low < 0.0f || vtln_low <= low_freq || vtln_low >= high_freq || 170 vtln_high <= 0.0f || vtln_high >= high_freq || vtln_high <= vtln_low)) { 171 KNF_LOG(FATAL) << "Bad values in options: vtln-low " << vtln_low 172 << " and vtln-high " << vtln_high << ", versus " 173 << "low-freq " << low_freq << " and high-freq " << high_freq; 174 } 175 176 bins_.resize(num_bins); 177 178 for (int32_t bin = 0; bin < num_bins; ++bin) { 179 float left_mel = mel_low_freq + bin * mel_freq_delta, 180 center_mel = mel_low_freq + (bin + 1) * mel_freq_delta, 181 right_mel = mel_low_freq + (bin + 2) * mel_freq_delta; 182 183 if (vtln_warp_factor != 1.0f) { 184 left_mel = VtlnWarpMelFreq(vtln_low, vtln_high, low_freq, high_freq, 185 vtln_warp_factor, left_mel); 186 center_mel = VtlnWarpMelFreq(vtln_low, vtln_high, low_freq, high_freq, 187 vtln_warp_factor, center_mel); 188 right_mel = VtlnWarpMelFreq(vtln_low, vtln_high, low_freq, high_freq, 189 vtln_warp_factor, right_mel); 190 } 191 192 // this_bin will be a vector of coefficients that is only 193 // nonzero where this mel bin is active. 194 std::vector<float> this_bin(num_fft_bins); 195 196 int32_t first_index = -1, last_index = -1; 197 for (int32_t i = 0; i < num_fft_bins; ++i) { 198 float freq = (fft_bin_width * i); // Center frequency of this fft 199 // bin. 200 float mel = MelScale(freq); 201 if (mel > left_mel && mel < right_mel) { 202 float weight; 203 if (mel <= center_mel) { 204 weight = (mel - left_mel) / (center_mel - left_mel); 205 } else { 206 weight = (right_mel - mel) / (right_mel - center_mel); 207 } 208 this_bin[i] = weight; 209 if (first_index == -1) { 210 first_index = i; 211 } 212 last_index = i; 213 } 214 } 215 KNF_CHECK(first_index != -1 && last_index >= first_index && 216 "You may have set num_mel_bins too large."); 217 218 bins_[bin].first = first_index; 219 int32_t size = last_index + 1 - first_index; 220 bins_[bin].second.insert(bins_[bin].second.end(), 221 this_bin.begin() + first_index, 222 this_bin.begin() + first_index + size); 223 224 // Replicate a bug in HTK, for testing purposes. 225 if (opts.htk_mode && bin == 0 && mel_low_freq != 0.0f) { 226 bins_[bin].second[0] = 0.0; 227 } 228 } // for (int32_t bin = 0; bin < num_bins; ++bin) { 229 230 if (debug_) { 231 std::ostringstream os; 232 for (size_t i = 0; i < bins_.size(); i++) { 233 os << "bin " << i << ", offset = " << bins_[i].first << ", vec = "; 234 for (auto k : bins_[i].second) os << k << ", "; 235 os << "\n"; 236 } 237 KNF_LOG(INFO) << os.str(); 238 } 239 } 240 241 void MelBanks::InitLibrosaMelBanks(const MelBanksOptions &opts, 242 const FrameExtractionOptions &frame_opts, 243 float vtln_warp_factor) { 244 htk_mode_ = opts.htk_mode; 245 int32_t num_bins = opts.num_bins; 246 if (num_bins < 3) { 247 KNF_LOG(FATAL) << "Must have at least 3 mel bins"; 248 } 249 250 float sample_freq = frame_opts.samp_freq; 251 int32_t window_length_padded = frame_opts.PaddedWindowSize(); 252 KNF_CHECK_EQ(window_length_padded % 2, 0); 253 254 int32_t num_fft_bins = window_length_padded / 2; 255 float nyquist = 0.5f * sample_freq; 256 257 float low_freq = opts.low_freq, high_freq; 258 if (opts.high_freq > 0.0f) { 259 high_freq = opts.high_freq; 260 } else { 261 high_freq = nyquist + opts.high_freq; 262 } 263 264 if (low_freq < 0.0f || low_freq >= nyquist || high_freq <= 0.0f || 265 high_freq > nyquist || high_freq <= low_freq) { 266 KNF_LOG(FATAL) << "Bad values in options: low-freq " << low_freq 267 << " and high-freq " << high_freq << " vs. nyquist " 268 << nyquist; 269 } 270 271 float fft_bin_width = sample_freq / window_length_padded; 272 273 float mel_low_freq = MelScaleSlaney(low_freq); 274 float mel_high_freq = MelScaleSlaney(high_freq); 275 276 debug_ = opts.debug_mel; 277 278 // divide by num_bins+1 in next line because of end-effects where the bins 279 // spread out to the sides. 280 float mel_freq_delta = (mel_high_freq - mel_low_freq) / (num_bins + 1); 281 282 bool slaney_norm = false; 283 if (!opts.norm.empty()) { 284 if (opts.norm != "slaney") { 285 KNF_LOG(FATAL) << "Unsupported norm: " << opts.norm; 286 } 287 slaney_norm = true; 288 } 289 290 bins_.resize(num_bins); 291 for (int32_t bin = 0; bin < num_bins; ++bin) { 292 float left_mel = mel_low_freq + bin * mel_freq_delta; 293 float center_mel = mel_low_freq + (bin + 1) * mel_freq_delta; 294 float right_mel = mel_low_freq + (bin + 2) * mel_freq_delta; 295 296 float left_hz = InverseMelScaleSlaney(left_mel); 297 float center_hz = InverseMelScaleSlaney(center_mel); 298 float right_hz = InverseMelScaleSlaney(right_mel); 299 300 // this_bin will be a vector of coefficients that is only 301 // nonzero where this mel bin is active. 302 // 303 // It is not an error to use num_fft_bins + 1 here. It is different 304 // from Kaldi. 305 std::vector<float> this_bin(num_fft_bins + 1); 306 307 int32_t first_index = -1, last_index = -1; 308 for (int32_t i = 0; i < num_fft_bins + 1; ++i) { 309 float hz = (fft_bin_width * i); // Center frequency of this fft bin. 310 if (hz > left_hz && hz < right_hz) { 311 float weight; 312 if (hz <= center_hz) { 313 weight = (hz - left_hz) / (center_hz - left_hz); 314 } else { 315 weight = (right_hz - hz) / (right_hz - center_hz); 316 } 317 318 if (slaney_norm) { 319 weight *= 2 / (right_hz - left_hz); 320 } 321 322 this_bin[i] = weight; 323 if (first_index == -1) { 324 first_index = i; 325 } 326 327 last_index = i; 328 } 329 } // for (int32_t i = 0; i < num_fft_bins + 1; ++i) 330 331 KNF_CHECK(first_index != -1 && last_index >= first_index && 332 "You may have set num_mel_bins too large."); 333 334 bins_[bin].first = first_index; 335 int32_t size = last_index + 1 - first_index; 336 bins_[bin].second.insert(bins_[bin].second.end(), 337 this_bin.begin() + first_index, 338 this_bin.begin() + first_index + size); 339 } // for (int32_t bin = 0; bin < num_bins; ++bin) 340 341 if (debug_) { 342 std::ostringstream os; 343 for (size_t i = 0; i < bins_.size(); i++) { 344 os << "bin " << i << ", offset = " << bins_[i].first << ", vec = "; 345 for (auto k : bins_[i].second) os << k << ", "; 346 os << "\n"; 347 } 348 fprintf(stderr, "%s\n", os.str().c_str()); 349 } 350 } 351 352 MelBanks::MelBanks(const float *weights, int32_t num_rows, int32_t num_cols) 353 : debug_(false), htk_mode_(false) { 354 bins_.resize(num_rows); 355 for (int32_t bin = 0; bin < num_rows; ++bin) { 356 const float *this_bin = weights + bin * num_cols; 357 358 int32_t first_index = -1, last_index = -1; 359 360 for (int32_t i = 0; i < num_cols; ++i) { 361 if (this_bin[i] == 0) { 362 continue; 363 } 364 if (first_index == -1) first_index = i; 365 last_index = i; 366 } 367 368 KNF_CHECK(first_index != -1 && last_index >= first_index && 369 "You have an incorrect weight matrix."); 370 371 bins_[bin].first = first_index; 372 int32_t size = last_index + 1 - first_index; 373 374 bins_[bin].second.insert(bins_[bin].second.end(), this_bin + first_index, 375 this_bin + first_index + size); 376 } 377 } 378 379 // "power_spectrum" contains fft energies. 380 void MelBanks::Compute(const float *power_spectrum, 381 float *mel_energies_out) const { 382 int32_t num_bins = bins_.size(); 383 384 for (int32_t i = 0; i < num_bins; i++) { 385 int32_t offset = bins_[i].first; 386 const auto &v = bins_[i].second; 387 float energy = 0; 388 for (int32_t k = 0; k != v.size(); ++k) { 389 energy += v[k] * power_spectrum[k + offset]; 390 } 391 392 // HTK-like flooring- for testing purposes (we prefer dither) 393 if (htk_mode_ && energy < 1.0) { 394 energy = 1.0; 395 } 396 397 mel_energies_out[i] = energy; 398 399 // The following assert was added due to a problem with OpenBlas that 400 // we had at one point (it was a bug in that library). Just to detect 401 // it early. 402 KNF_CHECK_EQ(energy, energy); // check that energy is not nan 403 } 404 405 if (debug_) { 406 fprintf(stderr, "MEL BANKS:\n"); 407 for (int32_t i = 0; i < num_bins; i++) 408 fprintf(stderr, " %f", mel_energies_out[i]); 409 fprintf(stderr, "\n"); 410 } 411 } 412 413 void ComputeLifterCoeffs(float Q, std::vector<float> *coeffs) { 414 // Compute liftering coefficients (scaling on cepstral coeffs) 415 // coeffs are numbered slightly differently from HTK: the zeroth 416 // index is C0, which is not affected. 417 for (int32_t i = 0; i != static_cast<int32_t>(coeffs->size()); ++i) { 418 (*coeffs)[i] = 1.0 + 0.5 * Q * sin(M_PI * i / Q); 419 } 420 } 421 422 } // namespace knf