mcelp 1.0.1

Mitsubishi CELP speech codec: a 3.6 kbit/s speech encoder and decoder
Documentation 
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

//! Critical-band analysis of the magnitude spectrum.
//!
//! The 129 magnitude bins are grouped into twenty bands of increasing width,
//! averaged, and converted to decibels.  The band levels are what the rest of
//! the encoder's spectral machinery works on.

use crate::analysis::BINS;
pub use crate::fixed::DB_PER_OCTAVE;
use crate::fixed::{acc, exp, hi, low, mul32x16, sat, shift, square32};
use crate::tables::{BAND_EDGES, BAND_INV_WIDTH_HALVED, LOG2_FRACTION};
use std::ops::RangeInclusive;

/// Number of critical bands.
pub const BANDS: usize = 20;

/// Inclusive spectrum-bin range belonging to one critical band.
pub(crate) fn band_bins(band: usize) -> RangeInclusive<usize> {
    BAND_EDGES[2 * band] as usize..=BAND_EDGES[2 * band + 1] as usize
}

/// Floor applied to every band level.
const LEVEL_FLOOR: i16 = 2560;
/// Offset subtracted from the log before scaling.
const LOG_OFFSET: i16 = 4;

/// Base-2 logarithm of a positive accumulator.
///
/// Returns `log2(x)` with the integer part in the high half and a
/// table-interpolated fraction below it; zero for a non-positive input.
pub fn log2(value: i64) -> i64 {
    if value <= 0 {
        return 0;
    }
    let e = exp(value);
    let integer = 30 - e;
    let normalised = shift(value, e);

    let shifted = shift(normalised, -9);
    let index = low(shift(shifted, -16)) as usize;
    let frac = (shift(shifted, -1) & 32767) as i16;

    let base = (LOG2_FRACTION[index] as i64) << 16;
    let slope = acc(base - ((LOG2_FRACTION[index + 1] as i64) << 16));
    let interpolated = shift(acc(base - (hi(slope) as i64) * (frac as i64) * 2), -16);
    acc(((integer as i64) << 16) + shift(interpolated, 1))
}

/// Convert a linear magnitude or energy to the codec's clipped dB scale.
fn decibel_level(value: i64, offset: i16) -> i64 {
    let log = log2(value);
    let level = acc(shift(log, -1) - ((offset as i64) << 16));
    let db = sat(shift(mul32x16(level, DB_PER_OCTAVE), 12));
    db.max((LEVEL_FLOOR as i64) << 16)
}

/// Mean magnitude of one critical band, converted to decibels.
fn spectrum_level(spectrum: &[i16; BINS], band: usize) -> i64 {
    let mut sum = 0i64;
    for bin in band_bins(band) {
        sum = acc(sum + ((spectrum[bin] as i64) << 16));
    }
    let mean = mul32x16(sat(sum), BAND_INV_WIDTH_HALVED[band]);
    decibel_level(shift(mean, 1), LOG_OFFSET)
}

/// Mean level of each critical band, in decibels.
///
/// `spectrum` is the magnitude spectrum after it has been scaled back by the
/// analysis headroom.
pub fn levels(spectrum: &[i16; BINS]) -> [i64; BANDS] {
    std::array::from_fn(|band| spectrum_level(spectrum, band))
}

/// Offset subtracted from the log of a tracked band energy.
const TRACKED_OFFSET: i16 = 6;
/// Weight each band contributes to the summary level, Q15 (halved because the
/// multiply doubles).
const SUMMARY_WEIGHT: i16 = 1638;
/// Upper clip on the level difference.
const DIFF_MAX: i16 = 20480;
/// Lower clip on the level difference.
const DIFF_MIN: i16 = -5120;

/// Convert one band energy to decibels.
///
/// The same shaping as [`levels`], but applied to an energy that is already
/// tracked at full 32-bit precision rather than summed from the spectrum.
pub fn tracked_level(energy: i64) -> i64 {
    decibel_level(energy, TRACKED_OFFSET)
}

/// Convert a whole set of tracked band energies.
pub fn tracked_levels(energies: &[i64; BANDS]) -> [i64; BANDS] {
    std::array::from_fn(|band| tracked_level(energies[band]))
}

/// Weighted sum of a set of band levels.
pub fn summary(levels: &[i64; BANDS]) -> i64 {
    let mut total = 0i64;
    for &l in levels.iter() {
        total = acc(total + mul32x16(l, SUMMARY_WEIGHT));
    }
    total
}

/// Difference between the current and the tracked summary level, clipped.
pub fn summary_difference(current: i64, tracked: i64) -> i64 {
    let diff = acc(current - tracked);
    diff.min((DIFF_MAX as i64) << 16)
        .max((DIFF_MIN as i64) << 16)
}

/// Upper clip on a per-band excess level.
const EXCESS_MAX: i16 = 20480;

/// Per-band excess of the current level over a tracked reference, clipped to
/// `0..=20480`.
///
/// The two tracked sets give two different references, so this is run twice.
pub fn excess(levels: &[i64; BANDS], tracked: &[i64; BANDS]) -> [i64; BANDS] {
    let mut out = [0i64; BANDS];
    for i in 0..BANDS {
        let d = acc(levels[i] - tracked[i]);
        out[i] = d.min((EXCESS_MAX as i64) << 16).max(0);
    }
    out
}

/// Bands taking part in the SNR statistics (band 0 is left out).
const STAT_BANDS: usize = 18;
/// Weight applied when summing the per-band SNR, Q15.
const STAT_MEAN_WEIGHT: i16 = 8192;
/// `1/18` in Q15, used to turn the sums into a variance.
const STAT_INV_BANDS: i16 = 1820;

/// Weighted sum and variance of the per-band SNR.
///
/// Band 0 is excluded.  These are the spectral-flatness style features the
/// encoder's voice-activity logic works from.
pub fn snr_statistics(snr: &[i64; BANDS]) -> (i64, i64) {
    let mut sum = 0i64;
    let mut sum_sq = 0i64;
    for &v in snr[1..1 + STAT_BANDS].iter() {
        sum = acc(sum + mul32x16(v, STAT_MEAN_WEIGHT));
        sum_sq = acc(sum_sq + square32(v));
    }
    let scaled = shift(sum_sq, -4);
    let mean_sq = mul32x16(square32(sum), STAT_INV_BANDS);
    let variance = shift(mul32x16(acc(scaled - mean_sq), STAT_INV_BANDS), 4);
    (sum, variance)
}