matlab freq2mel,[Kaldi] 特征提取--MFCC(二)

20180701qzd

本章讲解mfcc理论知识

一 基本含义

MFCC是Mel-Frequency Cepstral Coefficients的缩写，顾名思义MFCC特征提取包含两个关键步骤：转化到梅尔频率，然后进行倒谱分析。

1. 梅尔频率

梅尔刻度是一种基于人耳对等距的音高(pitch)变化的感官判断而定的非线性频率刻度。和频率的赫兹的关系如下:

所以当在梅尔刻度上面上是均匀分度的话，对于的赫兹之间的距离将会越来越大，所以梅尔刻度的滤波器组的尺度变化如下：

梅尔刻度的滤波器组在低频部分的分辨率高，跟人耳的听觉特性是相符的，这也是梅尔刻度的物理意义所在。

这一步的含义是：首先对时域信号进行傅里叶变换转换到频域，然后再利用梅尔频率刻度的滤波器组对应频域信号进行切分，最后每个频率段对应一个数值。

2. 倒谱分析

倒谱的含义是：对时域信号做傅里叶变换，然后取log，然后再进行反傅里叶变换。可以分为复倒谱、实倒谱和功率倒谱，我们用的是功率倒谱。

倒谱分析可用于将信号分解，两个信号的卷积转化为两个信号的相加。举例如下：

假设上面的频率谱X(k)，时域信号为x(n)那么满足

考虑将频域X(k)拆分为两部分的乘积：

假设两部分对应的时域信号分别是h(n)和e(n)，那么满足：

此时我们是无法区分开h(n)和e(n)。

对频域两边取log：

然后进行反傅里叶变换：

假设此时得到的时域信号如下：

虽然此时获得时域信号x’(n)即为倒谱，已经和原始的时域信号x(n)不一样，但是可以把时域信号的卷积关系转化为了线性加关系。

对应上图的频域信号，可以拆分成两部分的乘积：频谱的包络和频谱的细节。频谱的峰值即为共振峰，它决定了信号频域的包络，是辨别声音的重要信息，所以进行倒谱分析目的就是获得频谱的包络信息。包络部分对应的是频谱的低频信息，而细节部分对应的是频谱的高频信息。倒谱分析已经将两部分对应的时域信号的卷积关系转化为了线性加关系，所以只需要将倒谱通过一个低通滤波器即可获得包络部分对应的时域信号h’(t)。

二 基本流程

1. 预加重

目的

为了消除发声过程中，声带和嘴唇造成的效应，来补偿语音信号受到发音系统所压抑的高频部分。并且能突显高频的共振峰。

简单理解就是在频域上面都乘以一个系数，这个系数跟频率成正相关，所以高频的幅值会有所提升。实际上就是通过了一个H(z)=1−kz−1H(z)=1−kz−1高通滤波器。

实现

2. 加窗

目的

用于平滑信号，使用汉明窗加以平滑的话，相比于矩形窗函数，会减弱FFT以后旁瓣大小以及频谱泄露。

实现

使用汉明窗对信号进行加窗处理

3. 频域转换

目的

将时域信号转化到频域进行后续的频率分析

实现

幅度谱：

功率谱：

4. 使用梅尔刻度滤波器组过滤

目的

因为频域信号有很多冗余，滤波器组可以对频域的幅值进行精简，每一个频段用一个值来表示。

实现

对于FFT得到的幅度谱，分别跟每一个滤波器进行频率相乘累加，得到的值即为该帧数据在在该滤波器对应频段的能量值。如果滤波器的个数为22，那么此时应该得到22个能量值

5. 能量值取log

由于人耳对声音的感知并不是线性的，用log这种非线性关系更好描述。取完log以后才可以进行倒谱分析。

6. 离散余弦变换

DCT和DFT类似，但是只使用实数，不涉及复数运算。信号经过DCT变换以后，能量会集中到低频部分，可以用于图像压缩。

目的

按照倒谱的定义，该步需要进行反傅里叶变换然后通过低通滤波器获得最后的低频信号。这里使用DCT直接就可以获取频率谱的低频信息。

由于滤波器之间是有重叠的，所以前面的获得的能量值之间是具有相关性的，DCT还可以对数据进行降维压缩和抽象，获得最后的特征参数。相比于傅里叶变换，离散余弦变换的结果没有虚部，更好计算。

实现

7. 差分

目的

由于语音信号是时域连续的，分帧提取的特征信息只反应了本帧语音的特性，为了使特征更能体现时域连续性，可以在特征维度增加前后帧信息的维度。常用的是一阶差分和二阶差分。

实现

附： matlab代码

%% clear

clear all; close all; clc;

%% Pre-process

% parameter definition

filename = 'sx438.wav';

frame_length_ms = 25;

frame_shift_ms = 10;

dc_offset = 300;

pre_emphasis_coef = 0.97;

rand_noise_factor = 0.6;

add_rand_noise = true;

% step 1 : read wav

wav = audioread(filename, 'native');

figure;

subplot(2,1,1); plot(wav); title('waveform'); hold on;

% step 2 : add dc offset

wav = wav + dc_offset;

subplot(2,1,2); plot(wav); title('add dc offset'); hold on;

% step 3 : extract frame

info = audioinfo(filename); info

num_of_frames = floor((info.Duration * 1000 - frame_length_ms) / frame_shift_ms); % avoid out of range

num_of_per_frame = (info.SampleRate * frame_length_ms / 1000);

frames = zeros(num_of_frames, num_of_per_frame);

for f = 1:num_of_frames

idx = (f - 1) * (frame_shift_ms * info.SampleRate / 1000) + 1; % addrees start from 1

frames(f,:) = wav(idx:idx+num_of_per_frame-1)';

if add_rand_noise

for i = 1:num_of_per_frame

frames(f,i) = frames(f,i) - 1.0 + 2 * rand(); % rand [-1, 1)

end

end

end

figure;

subplot(3,1,1); plot(frames(51,:)); title('pre-sub-mean'); hold on;

% step 4 : sub mean

for f = 1:num_of_frames

me = mean(frames(f,:));

frames(f,:) = frames(f,:) - me;

end

subplot(3,1,2); plot(frames(51,:)); title('post-sub-mean'); hold on;

% step 5 : pre-emphasis

for f = 1:num_of_frames

for i = num_of_per_frame:-1:2

frames(f,i) = frames(f,i) - pre_emphasis_coef * frames(f,i-1);

end

end

subplot(3,1,3); plot(frames(51,:)); title('post-pre-emph'); hold on;

% step 6 : hamming window

hamming = zeros(1, num_of_per_frame);

for i = 1:num_of_per_frame

hamming(1,i) = 0.54 - 0.46*cos(2*pi*(i-1) / (num_of_per_frame-1));

end

figure;

plot(hamming); title('Hamming Window Function'); hold on;

% step 7 : add window

figure;

subplot(2,1,1); plot(1:400, frames(51,:),'r', 1:400, hamming, 'g'); title('pre-add-win'); hold on;

for f = 1:num_of_frames

frames(f,:) = frames(f,:) .* hamming;

end

subplot(2,1,2); plot(1:400, frames(51,:)); title('post-add-win'); hold on;

%% FbankFilter Analysis

% fft

NFFT = 2^nextpow2(num_of_per_frame); % Next power of 2 from length of num_of_per_frame

frames_fft = complex(zeros(num_of_frames, NFFT));

for f = 1:num_of_frames

frames_fft(f,:) = fft(frames(f,:), NFFT); % pad 0

end

freq = info.SampleRate / 2 * linspace(0, 1, NFFT/2 + 1); % frequency of fft

figure; plot(freq, 2*abs(frames_fft(1,1:NFFT/2+1))); hold on;

xlabel('Frequency(Hz)'); ylabel('|frame\_fft|');

% compute energy

energy_frames = zeros(num_of_frames, NFFT/2);

for f = 1:num_of_frames

for i = 1:NFFT/2

energy_frames(f,i) = abs(frames_fft(f,i));

end

end

% mel filter coef

num_of_bins = 40;

low_freq = 20; high_freq = info.SampleRate / 2;

fft_bin_width = info.SampleRate / NFFT;

mel_low_freq = 2595*log10(1 + low_freq / 700);

mel_high_freq = 2595*log10(1 + high_freq / 700);

mel_freq_delta = (mel_high_freq - mel_low_freq) / (num_of_bins + 1);

mel_coef = zeros(num_of_bins, (NFFT/2)+3);

for b = 1:num_of_bins

left_mel = mel_low_freq + (b - 1) * mel_freq_delta;

center_mel = mel_low_freq + (b) * mel_freq_delta;

right_mel = mel_low_freq + (b + 1) * mel_freq_delta;

first_index = -1;

last_index = -1;

for f = 1:(NFFT/2)

tmp_freq = (fft_bin_width * (f - 1)) + 1;

tmp_mel_freq = 2595*log10(1 + tmp_freq / 700);

if tmp_mel_freq > left_mel && tmp_mel_freq < right_mel

weight = 0;

if tmp_mel_freq < center_mel

weight = (tmp_mel_freq - left_mel) / (center_mel - left_mel);

elseif tmp_mel_freq < right_mel

weight = (right_mel - tmp_mel_freq) / (right_mel - center_mel);

end

if first_index == -1

first_index = f;

end

last_index = f;

mel_coef(b, f) = weight;

end

end

mel_coef(b, (NFFT/2)+2) = first_index;

mel_coef(b, (NFFT/2)+3) = last_index;

end

figure;

for b = 1:num_of_bins

plot(mel_coef(b,1:(NFFT/2))); hold on;

end

xlabel('Frequency'); ylabel('mel_coef');

% take mel filter

mel_frames = zeros(num_of_frames, num_of_bins);

for f = 1:num_of_frames

for i = 1:num_of_bins

mel_frames(f,i) = sum(energy_frames(f,:) .* mel_coef(i,1:NFFT/2));

end

end

% take log

fbank_frames = log(mel_frames);

%% MFCC

N = num_of_bins; M = num_of_bins; % assume M == N == num_of_bins

mfcc_frames = zeros(num_of_frames, M);

dct_mat = discrete_cos_trans(N, M);

for f = 1:num_of_frames

mfcc_frames(f,:) = sqrt(2/N) * fbank_frames(f,:) * (dct_mat');

end

% mean

m = mean(mfcc_frames);

% variance

v = sqrt(mean(power(mfcc_frames, 2)) - power(m,2) + 0.0001); % avoid div 0

% cmvn

for f = 1:num_of_frames

mfcc_frames(f,:) = (mfcc_frames(f,:) - m) ./ v;

end