[Python] MFCC 구현

ㅇ Process

1. 파일 불러오기

2. 전처리 (증폭)

3. Framing (프레임화)

4. Windowing (프레임 적용)

5. Fourier Transform (FT)

6. Mel-Filter 적용

7. Discrete Cosine Transform (DST)

8. MFCC

9. 출력

 

---

1. 파일 불러오기

 

ㅇ 음성파일을 byte 코드로 읽어서 사용

 

- 아래 예는 8000 sample_rate, 16bit 음성 데이터 사용

- 즉 signal의 index는 시간(1/진동수), value는 magnitude를 의미

f=open(target, 'rb')
buf = bytearray(f.read())
pcm_data = numpy.frombuffer(buf, dtype = 'int16')
signal = librosa.util.buf_to_float(x=pcm_data, n_bytes=2)   # pcm to wav
x_axis = [i/sample_rate if i>0 else i for i in range(0, len(signal))]   # convert to second(time) level
sample_len = len(signal)/sample_rate        # total time

 

---

 

2. 전처리 (증폭)

 

ㅇ 음성 인식 증대 및 SNR (Signal to Noise) 감소 목적으로 고주파 영역 증폭

y(t)=x(t)−αx(t−1)

pre_emphasis		= 0.97
emphasized_signal	= numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])
signal_length		= len(emphasized_signal)

 

---

3. Framing

 

ㅇ 프레임 크기

- 프레임 크기(size)는 일반적으로 20~40ms (0.02~0.04초) 로 정의, 사람이 말하는 음소가 바뀔 수 있는 최소 시간이기 때문

 

- frame_stride = 프레임을 얼마의 시간 간격으로 겹쳐나갈 것인지 결정

frame_size      = 0.025
frame_stride    = 0.01

 

ㅇ 프레임 단위 변환

- 앞서 프레임 크기를 시간으로 정했으나 실제 음성 데이터는 진동수를 기반으로 정의되어 있기에, 정의된 시간에 맞는 진동수로 변환

- 비례식 사용

1 sec : 8000 Hz = frame_size : frame_length

frame_length    = int(round(frame_size * sample_rate))		# 0.025초 동안의 진동수
frame_step      = int(round(frame_stride * sample_rate))	# 0.01초 동안의 진동수
num_frames	= int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))

 

ㅇ 최적 프레임 수

- 최소 1개의 프레임 (frame_length)을 보장하고, 나머지 구역을 frame_step (Hz 단위)로 나누었을 때 몇 개의 프레임으로 나눌 수 있는지 계산

num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))

 

ㅇ 프레임 동기화

- 최적 프레임 수(num_frames)로 나누는 과정에서 정확히 나누어 떨어지지 않는 경우가 발생하기에, 해당 구역 또한 cover하기 위한 목적으로 패딩 구역 반영

pad_signal_length = num_frames * frame_step + frame_length
diff_len = numpy.zeros((pad_signal_length - signal_length))
pad_signal = numpy.append(emphasized_signal, diff_len)

 

---

4. Windowing

 

ㅇ 프레임에 맞는 윈도우 생성 및 Hamming 적용

- indices_frame = 프레임을 1 간격으로 분할 (즉, [0, ... , frame_length] 인 배열을 'num_frames"만큼 행으로 가지는 배열

- indices_step = 전체 프레임 영역을 'frame_step'으로 분할한 배열을 'frame_length'만큼 행으로 가지는 배열

- indices = 전체 프레임 구역에서 'frame_step' 간격으로 'frame_length' 크기의 프레임들을 겹친다고 할 때, 순차로 나열되는 프레임의 순번이 index이고, 내부 값은 해당 프레임 시작점부터 'frame_length'의 크기를 가지는 배열

예) frame_step=80, frame_length=200, signal

 

 

을 나눌 때, 각 프레임 번호가 row index에 해당하고, 해당 프레임

indices_frame = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1))
indices_step = numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1))
indices = indices_frame + indices_step.T
frames = pad_signal[indices.astype(numpy.int32, copy=False)]
frames *= numpy.hamming(frame_length)

 

 

 

ㅇ 임시 코드 (정리중...)

import os
import numpy
import librosa
import librosa.display
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from matplotlib.figure import Figure
from scipy.fftpack import dct

target = os.getcwd()+"/sample/ffd0af0e-607e-4c5a-b_20220119100421824_1684161253_55_011_0_combined.pcm"
# target = os.getcwd()+"/sample/2f24974a-2c3a-4b9c-b_20220119094833820_0930447499_46_020_0_combined.pcm"
sample_rate = 8000

### PCM
f=open(target, 'rb')
buf = bytearray(f.read())
pcm_data = numpy.frombuffer(buf, dtype = 'int16')
signal = librosa.util.buf_to_float(x=pcm_data, n_bytes=2)   # wav
x_axis = [i/sample_rate if i>0 else i for i in range(0, len(signal))]   # convert to second(time) level
sample_len = len(signal)/sample_rate        # total time



### pre-Emphasis
pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])



### Framing
frame_size      = 0.025     # 20ms ~ 40ms
frame_stride    = 0.01      # frame_size's 50% +-10%

signal_length   = len(emphasized_signal)
# proportionality. original : frame
frame_length    = int(round(frame_size * sample_rate))      # 1(sec):sample_rate(Hz) = frame_size(sec):frame_length(Hz)
frame_step      = int(round(frame_stride * sample_rate))    # 1(sec):sample_rate(Hz) = frame_stride(sec):frame_step(Hz)

num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))  # Make sure that we have at least 1 frame
# Pad Signal to make sure that all frames have equal number of samples without truncating any samples from the original signal
pad_signal_length = num_frames * frame_step + frame_length
z = numpy.zeros((pad_signal_length - signal_length))
pad_signal = numpy.append(emphasized_signal, z)



### Overlap
indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
frames = pad_signal[indices.astype(numpy.int32, copy=False)]
frames *= numpy.hamming(frame_length)   # hamming


print("{} {}".format(num_frames, pad_signal_length))
print("signal_len {}, frame_len {}, frame_step {}, frame_num {}".format(signal_length, frame_length, frame_step, num_frames))
print("pad_signal_len {}".format(pad_signal_length))



### Fourier Transform
NFFT = 268  # 256 or 512
dft_frames = numpy.fft.rfft(frames, NFFT)           # Discrete Fourier Transform     # Complex number
# dft_frames = numpy.fft.rfft(frames)           # Discrete Fourier Transform     # Complex number
# NFFT = dft_frames.shape[1]
mag_frames = numpy.absolute(dft_frames)             # Magnitude of the FFT          # remove phase, olny left magnitude info.
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))   # Power Spectrum
print("{} : frames".format(frames.shape))
print("{} : dft_frames".format(dft_frames.shape))
print("{} : mag_frames".format(mag_frames.shape))
print("{} : pow_frames".format(pow_frames.shape))
print(mag_frames)

### Filter Bank
nfilt = 40      # typically use 40 filters
low_freq_mel, high_freq_mel = 0, (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # Convert Hz to Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)
hz_points = (700 * (10**(mel_points / 2595) - 1))           # Convert Mel to Hz
bin = numpy.floor((NFFT + 1) * hz_points / sample_rate)
fbank = numpy.zeros((nfilt, int(numpy.floor(NFFT / 2 + 1))))
for m in range(1, nfilt + 1):
    f_m_minus = int(bin[m - 1])   # left
    f_m = int(bin[m])             # center
    f_m_plus = int(bin[m + 1])    # right

    for k in range(f_m_minus, f_m):
        fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
    for k in range(f_m, f_m_plus):
        fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])

filter_banks = numpy.dot(pow_frames, fbank.T)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # Numerical Stability
filter_banks = 20 * numpy.log10(filter_banks)  # dB # log-mel spectrum
print("{} : filter_banks".format(filter_banks.shape))

### MFCC
num_ceps = 12
# mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1 : (num_ceps + 1)]   # Discrete Cosine Transform (DCT)
mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')
print("{} : mfcc".format(mfcc.shape))

###
fig, (ax_sig, ax_fft, ax_spec) = plt.subplots(3)
fig:Figure
ax_sig:Axes
ax_fft:Axes
ax_spec:Axes


ax_sig.plot(x_axis, signal)
ax_fft.plot(numpy.linspace(0, sample_rate, len(mag_frames)), mag_frames)


# spc = librosa.display.specshow(filter_banks.T, x_axis='time', y_axis='mel', sr=sample_rate ,ax=ax_spec)
spc = librosa.display.specshow(filter_banks.T, x_axis='time', y_axis='mel', sr=sample_rate ,ax=ax_spec)
# plt.colorbar(format='%+2.0f dB')
fig.colorbar(spc, ax=ax_spec, format='%+2.0f dB')
ax_spec.set_title("Spectrogram")


plt.show()