DSP - Lab 3: MFCC: Mel 频率的倒谱系数

# Load audio signal from disk.
y, sr = load_audio(p)

# Pre-emphasize and normalize the signal.
y = np.append(y[0], y[1:] - pre_emphasis * y[:-1])
y = y / np.max(np.abs(y))

# Get the window length for FFT.
n_window = t_window * sr // 1000
n_window = utils.next_pow2(n_window)

# Detect the ranges of voice activity.
ranges, i_starts, zcr = detect_voice_activity(y, n_window)

n_samples = y.shape[0]
i_starts = np.arange(0, n_samples, n_window // 2, dtype=int)
i_starts: np.ndarray = i_starts[i_starts + n_window < n_samples]

avg_amps = [np.average(np.abs(y[i:i+n_window])) for i in i_starts]

avg_zcrs = [np.sum(np.abs(
    np.sign(y[i:i+n_window-1]) - np.sign(y[i+1:i+n_window])
)) / 2 / t_window for i in i_starts]

# Step 1: Find the ranges by judging whether the average amplitude is
# higher than threshold `amp_th[1]`.
ranges_1: List[List[int]] = []
for k, avg_amp in enumerate(avg_amps):
    if avg_amp > amp_th[1]:
        if len(ranges_1) > 0 and k <= ranges_1[-1][1] + 2:  # overlaps
            ranges_1[-1][1] = k
        else:
            ranges_1.append([k, k])

# Step 2: Expand the ranges by judging whether the average amplitude is
# higher than threshold `amp_th[0]`.
ranges_2: List[List[int]] = []
for r in ranges_1:
    i_start, i_stop = r
    i_stop_prev = ranges_2[-1][1] if len(ranges_2) > 0 else 0
    while i_start > i_stop_prev and avg_amps[i_start] > amp_th[0]:
        i_start -= 1
    while i_stop < len(i_starts) - 1 and avg_amps[i_stop] > amp_th[0]:
        i_stop += 1
    if i_start <= i_stop_prev + 2 and i_stop_prev != 0:  # overlaps
        ranges_2[-1][1] = i_stop
    else:
        ranges_2.append([i_start, i_stop])

# Step 3: Expand the ranges by judging whether the average zero-crossing
# rate (ZCR) is higher than threshold `zcr_th`.
ranges_3: List[List[int]] = []
for r in ranges_2:
    i_start, i_stop = r
    i_stop_prev = ranges_3[-1][1] if len(ranges_3) > 0 else 0
    i_start_min = max(i_stop_prev, r[0] - zcr_step_th)
    i_stop_max = min(len(i_starts) - 1, r[1] + zcr_step_th)
    while i_start > i_start_min and avg_zcrs[i_start] > zcr_th:
        i_start -= 1
    while i_stop < i_stop_max and avg_zcrs[i_stop] > zcr_th:
        i_stop += 1
    if i_start <= i_stop_prev + 2 and i_stop_prev != 0:  # overlaps
        ranges_3[-1][1] = i_stop
    else:
        ranges_3.append([i_start, i_stop])

ranges = [[i_starts[r[0]], i_starts[r[1]] + n_window] for r in ranges_3]

# Obtain the Mel filter banks.
f_min, f_max = 20, sr // 2
filters = get_mel_filters(
    n_mel_filters, sr, n_window, f_min, f_max,
)

def mel_freq(f: np.ndarray) -> np.ndarray:
    return 2595 * np.log10(1 + f / 700)

mel_f_min, mel_f_max = mel_freq(f_min), mel_freq(f_max)

mel_f = np.linspace(mel_f_min, mel_f_max, n_filters + 2)
f = np.floor(i_mel_freq(mel_f) * n_window / sr).astype(int)

filter_len = n_window // 2
filters = np.array([np.concatenate([
    np.zeros(f[i - 1]),
    np.linspace(0, 1, f[i] - f[i - 1], endpoint=False),
    np.linspace(1, 0, f[i + 1] - f[i], endpoint=False),
    np.zeros(filter_len - f[i + 1]),
]) for i in range(1, n_filters + 1)])

# Get the spectrogram using STFT.
r = ranges[0]
spec, i_starts = create_spectrogram(y[r[0]:r[1]], n_window)
energy_spec = np.square(spec)

# Filter the energy spectrum with the Mel filter banks.
filtered_spec = np.dot(filters, energy_spec)
log_filtered_spec = 10 * np.log10(filtered_spec)

# Generate the MFCC.
cc = dct(log_filtered_spec, dim_mfcc + 1)[1:]

def dct(x: np.ndarray, d: int) -> np.ndarray:
    '''
    Perform a Discrete Cosine Transform of a 1D / 2D array.

    Args:
        `x`: source array, shape(n, l)
        `d`: dimension of the DCT matrix

    Returns:
        The result of DCT, shape(d, l)
    '''

    n = x.shape[0]
    c = np.sqrt(2 / n)
    s = np.pi / (2 * n) * np.arange(1, 2 * n, 2)
    dct_m = np.concatenate([
        [np.ones(n) / np.sqrt(n)],
        [c * np.cos(i * s) for i in range(1, d)],
    ])
    return np.dot(dct_m, x)

conda env update --name dsp --file environment.yml
conda activate dsp

python3 main.py