MFCC를 하고있는데 왜 LPC를 하냐?

Matlab에서 포만트 관련 예제를 lpc로 하고있어서 입니다...

포만트에 대해서 간단히 언급한 후에 LPC를 통한 포만트 분석 후 이를 MFCC로 하겠습니다.

우선 소스

load mtlb;

segmentlen = 100;

noverlap = 90;

NFFT = 128;

[y,f,t,p] = spectrogram(mtlb,segmentlen,noverlap,NFFT,Fs);

surf(t,f,10*log10(abs(p)),'EdgeColor','none');

axis xy; axis tight; colormap(jet); view(0,90);

xlabel('Time');

ylabel('Frequency (Hz)');

dt = 1/Fs;

I0 = round(0.1/dt);

Iend = round(0.25/dt);

x = mtlb(I0:Iend);

x1 = x.*hamming(length(x));

preemph = [1 0.63];

x1 = filter(1,preemph,x1);

A = lpc(x1,8);

rts = roots(A);

rts = rts(imag(rts)>=0);

angz = atan2(imag(rts),real(rts));

[frqs,indices] = sort(angz.*(Fs/(2*pi)));

bw = -1/2*(Fs/(2*pi))*log(abs(rts(indices)));

nn = 1;

for kk = 1:length(frqs)

    if (frqs(kk) > 90 && bw(kk) <400)

        formants(nn) = frqs(kk);

        nn = nn+1;

    end

end

formants

------------------------------------------------------------

짧습니다. 실제 함수 하나하나 구현하려면 어렵겠지만요.

mtlb : 음원

segmentation : 분할... 이라는데. 해밍윈도우와 관련. 

noverlap : 앞서 MFCC에서의 overlap 부분.

NFFT : FFT resolution. 주파수 분해능, 

dt : 주기

IO와 len : 0.1~0.25초 때의 번호를 뜻합니다

x : mtlb의 0.1~0.25초 사이의 신호정보 입니다.

round는 가장 가까운 정수쪽으로 반올림을 합니다.

x1 : 해밍윈도우를 통과한 신호정보 입니다.

preemph : 프리엠퍼시스 필터의 적용입니다. 프리엠퍼시스 필터는 highpass all-pole(AR(1))필터라는데 all-pole이 무슨 뜻인지 모르겠네요. 이를 x1에 적용하였습니다.

필터관련은 나중에 더 공부하기로 하죠

A=lpc(x1,8) 을 통해서 신호 x1에 대하여 8개의  lpc ㄱ

이것저것 설치할 때, 관리자 권한을 일일히 획득하기 귀찮다.

제어판에서 사용자계정에 들어간 후 계정유형을 관리자로 변경하자.

cmd에서의 net user 계정명 /active: on off 사용가능

만약 관리자계정을 잘못만져서 아무 계정에도 접근하지 못하게 되었을 때의 해결방법.

windows시디로 부팅 후 shift+f10으로 cmd 접근

c드라이버로 접근 cd:\ c

windows\system32

copy utilman.exe utilman.exe.1

del utilman.exe

copy cmd.exe utilman.exe

재부팅후 계정 선택창에서 왼쪽 아래의 시계모양 클릭

net user administrator 비번 입력

DISKPART

컴퓨터 포맷 후 shift+F10 으로 명령프롬프트에 접근 가능하다.

디스크 포맷 및 파티션 관리는 다음 사이트에 자세히 기재되어 있다.

http://cappleblog.co.kr/

한시간이면 될것을 몇일 걸린 이유...

1. 선이 끊어져 있는데 이걸 발견하지 못했다..    :단선 확인은 필수다

2. GND끼리 연결하지 않았다.

MCU랑 연계해서 값을 읽는다고 한다면 하나의 회로를 이루어야 한다. 전원은 따로더라도 엔코더쪽 GND따로 MCU쪽 GND따로가 아니라는 말씀. 두 GND를 연결해주자. 그러지 않으면 값이 불규칙적이다.

엔코더가 제대로 돌아가는지 확인하는법.

엔코더에 전원을 연결해 주고, 엔코더의 A B상 출력과 GND사이에 각각 LED를 연결해준다.

엔코더를 돌릴때 두 LED가 깜빡여야 한다. 계속 켜져있거나 계속 꺼져있으면 불량이다.

//******************************************************************

헤더

#include <stm32f10x_it.h>

void Encoder_RCC(void);

void Encoder_GPIO(void);

void Encoder_TIM(void);

void Encoder_Setting(void);

void Encoder_NVIC(void);

void Encoder_Disable(void);

저번과 같은 맥락의 함수들이다

//*********************************************************************

본문

1. TIM3 사용, PortC로 Remapping할것임. 본래 타이머 외의 용도를 사용하기 위해서 AFIO를 Enable

void Encoder_RCC(void){

  //PortC TIM3

  RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOC|RCC_APB2Periph_AFIO,ENABLE);

  RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM3,ENABLE);

}

2. TIM3를 Full-Remap하여 PortC의 핀6~9를 TIM3로 사용

void Encoder_GPIO(void){

  GPIO_InitTypeDef GPIO_InitStructure;

  GPIO_PinRemapConfig ( GPIO_FullRemap_TIM3, ENABLE );

  GPIO_InitStructure.GPIO_Pin = GPIO_Pin_6|GPIO_Pin_7|GPIO_Pin_8|GPIO_Pin_9;

  GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;

  GPIO_InitStructure.GPIO_Mode = GPIO_Mode_IPU;//GPIO_Mode_IN_FLOATING;//

  GPIO_Init(GPIOC,&GPIO_InitStructure); 

}

3.

TIM_EncoderInterfaceConfig : 엔코더 설정 함수. 타이머, 엔코더 모드, A B 상의 펄스읽는 방식 설정

TIM_ARRPreloadConfig : Enables or disables the TIM peripheral Preload register on ARR.

//타이머카운터의 기본기능이 아닌 다른 기능을 사용하기 위한 함수라는데 모르겠다.

TIM_ClockDivision : 분주비. 본인은 4분주 사용

TIM_ITConfig : 타이머 인터럽트 활성화

void Encoder_TIM(void){

  TIM_EncoderInterfaceConfig(TIM3, TIM_EncoderMode_TI12, TIM_ICPolarity_Rising, TIM_ICPolarity_Rising);

  TIM_ARRPreloadConfig(TIM3, ENABLE);

  TIM_Cmd(TIM3, ENABLE);

  //***************************************************************

  TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;

  TIM_TimeBaseStructure.TIM_Period = 500;//260;//624;

  TIM_TimeBaseStructure.TIM_Prescaler = 1;//576;

  //

  TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV4;//TIM_CKD_DIV1;//TIM_CKD_DIV4;//

  TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up; //upcounting mode

  //TIM_TimeBaseStructure.TIM_RepetitionCounter = 0; 

  TIM_TimeBaseInit(TIM3, &TIM_TimeBaseStructure);

  TIM_ITConfig(TIM3, TIM_IT_Update, ENABLE); 

}

void Encoder_NVIC(void){

  NVIC_InitTypeDef NVIC_InitStructure;

    NVIC_InitStructure.NVIC_IRQChannel = TIM3_IRQn;            

    //NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 6;

    //NVIC_InitStructure.NVIC_IRQChannelSubPriority = 0;

    NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;

    NVIC_Init(&NVIC_InitStructure);

}

void Encoder_Setting(void){

  Encoder_RCC();

  Encoder_GPIO();

  //Encoder_NVIC(); //인터럽트 없었음.. 그 자체로 모터의 펄스값을 읽어들였음. 전에 시험테스트 해볼때도

  Encoder_TIM();

}

void Encoder_Disable(void){

  RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOC|RCC_APB2Periph_AFIO,DISABLE);

  RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM3,DISABLE);

  TIM_ARRPreloadConfig(TIM3, DISABLE);

  TIM_Cmd(TIM3, DISABLE);

  TIM_ITConfig(TIM3, TIM_IT_Update, DISABLE); 

}

지금 하고있는 졸작 프로젝트에 맞춘것임. 용도불명인 내용이 많을 수 있음.

Bluetooth.zip

헤더부터 살펴보면

1.

#define GPIO_Bluetooth  GPIOA

#define GPIO_Bluetooth_Rx_Pin GPIO_Pin_3

#define GPIO_Bluetooth_Tx_Pin GPIO_Pin_2

PortA의 2번과 3번핀을 사용하였다. Cortex회로도를 보면 알겠지만 각각 USART2의 TX와 RX이다.

//-----------------Buffer Size-------------------------

#define U2_BUFFER_SIZE 100

//----------------------------------------------------

void Blutooth_RCC_Configuration(void);

void Blutooth_NVIC_Configuration(void);

void Blutooth_GPIO_Configuration(void);

void Blutooth_UART2_Configuration(void);

void Blutooth_Setting(void);    //위의 RCC~UART2까지를 함수하나로 부르려고 한것이다.

//void USART2_IRQHandler(void);

void PutChar_blue(uint8_t c);        //블루투스에 한문자를 전송

void PutString_blue(uint8_t *s);    //문자단위로 문자열을 전송

void Send_Command_test(uint16_t Command);    //전송테스트용. 패킷을 위루어 전송

void Transmit_data(uint8_t *s);    //실질적으로 문자단위로 전송하는 부분

//********************Status Test**************************

bool input_check(void);

bool Connection_Check(void);

//********************Send Speed***************************

void Send_EncoderState(uint16_t Speed);

그냥보면은 용도를 알 수 없는 함수가 꽤 있다. 실제 내가 복습하려고 만들어 놓은 내용이니..

몇곳만 짚어 적어둔다.

1.NVIC : 실제로는 부르지 않게 해놓았다. 이걸 부르면 어쩐지 시리얼 통신부분에서 진행하질 못한다. 딱히 할 필요성은 못느낌

2.USART_ITConfig     수신에 대한 인터럽트를 가능하게 해준단다.

3.USART_SendData(USART2,c);    데이터 보낼 때

4.

intData=0;

  Command_Data[4]=(0x00ff)&intData;

  Command_Data[3]=((0xff00)&intData)>>8;

intData에 저장한 값을 PC가 받게 하려고 할 때, 다음과 같이 넣어서 보내면 의도한 값을 컴퓨터에서 받을 수 있다.

나같이 아무것도 모르는 채로 맨땅에 헤딩하는 사람들을 위한 블루투스 사용기...

대상은 펌테크의 FB155BC_SMD

메뉴얼을 보면 어떠한 명령이 있다던가 다 나와있지만, '정말 아무것도 모르는 사람'들을 위한 개략적인 설명.

인터페이스 보드가 없다. 하지만 설정은 하고싶다. 이때 AT명령어를 이용한 설정을 하면 된다고 한다.

AT명령어 관련은 관련 메뉴얼을 보면 되지만 'AT명령어를 어떻게 사용하지?'라는 것에 대한 언급은 없다...

디바이스마트에서 판매하는 단가 6000원? 정도의 RS232 통신모듈을 사용하면 된다. 또는 소자 몇개 가지고 스스로 만들면 된다.(본인의 경우 원치않게 만들어 사용했음) 모듈의 TX, RX를 블루투스의 RX, TX에 각각 연결해준다. 물론 블루투스와 RS232모듈 둘다 전원이 들어가 있어야 겠지... 테라텀으로는 자신이 전송하는 내용을 보는 방법을 모르겠으니,  ComportMaster 또는 하이퍼터미널로 신호를 보낸다. 하이퍼 터미널에서 자신이 보내는 내용을 보는 방법은 네이버에서 쉽게 찾을 수 있다(Property로 들어가면 된다) ComportMaster사용을 추천.

나머지는 메뉴얼대로 신호를 주고 받으면 된다.

팁1.

후배에게서 받은 팁.

공장에서 막 생산된 상태의 공장초기화를 하려면?

3번 Status에서 핀을 뽑아서, GND에 잠시 연결한 후 뗀다. 그러면 모든 설정이 초기화 된다.

팁2.

이 역시 후배에게서 받은 팁.

인터페이스 보드에서 나오는 메뉴를 사용하는 설정을 하고 싶은 경우.

3번 Status와 GND를 잠시 연결한 후 뗀다.

위에꺼랑 뭐가 다르냐면...

팁1도 팁2도 물론 전원이 들어간 상태에서 하는건데. 팁2의경우, Serial통신 창을 열어놓은 상태라면 알겠지만 GND에 데었다 떼면 인터페이스보드를 사용할때와 같은 메뉴가 나온다(이미 공장초기화 된 상태). 그상태에서 전원의 재인가 없이 설정을 진행한 후 종료했다 켜면 설정이 적용된 상태로 사용 가능한것. 팁1의 경우 전원을 재인가 해서 설정없이 말그대로 공장초기화 된것

스텝모터 돌리기 전에... 회로에 대한 기초 이론이 부족했기에.. 조금 조사해둔 내용

2상 여자라던가 이런부분은 알고있으니, 다른 내용 정리

D-CUT SHAFT : 모터 단독 제품. 그렇지 않을 경우 감속기 부착용

모터드라이버를 이용하여 모터에 전압과 전류를 공급한다.

그런데 공급하는 전압과 전류가 규정치를 넘어서 걱정이 된다. 이에 대한 내용

모터 스펙에 나와있는 전압은 정격전류를 흘리기 위한 최소전압이라 한다.

보통 표시된 전압의 4배정도의 전압을 공급한다.

표시된 전압만큼만 공급해도 규정 토크는 나오지만 최고 회전속도를 약간 더 올려주기 위해서,

그리고 높은 회전속도에서 토크저하를 줄이기 위해서 표시된 전압보다 4~8배 정도 높은 전압을 주는것이 유리하다.

모터드라이버 회로에서 모터의 규정전압을 넘는 전압을 인가하여도 토크는 동일하다.

예를들어 2V 1A모터에 16V를 인가한다 하더라도 전류가 1A를 넘는 순간 전압을 끊고, 조금 기다렸다가 다시 전원을 인가하는 동작(Chopping에 의한 전류제한)을 하기 때문에 항상 모터에는 1A의 평균/실효 전류가 흐른다.

모터의 자속은 전류X권선 회수에 직접 비례한다. 전압의 요소가 없다는 점을 상기하자

다만 고전압을 인가햇을 때 1A까지 전류가 상승하는 시간이 저전압보다 짧기 때문에 상대적으로 약간 더 토크가 크게 느껴지고 고속 회전에 유리해진다.

스텝모터는 정전류 구동형이다. 코일전압은 큰 의미가 없다.

전압이 높으면 홀딩토크는 커질것이나 일반적으로 크기 필요하지는 않다.

전압이 낮은쪽이 역기전력이 작다. 따라서 낮은 전압의 고 전류급 모터를 추천한다고 한다.(3V이상 되지 않도록)

파워서플라이 등으로 스텝모터를 구동할때.. 이런 질문이 있었다고 한다.

5V 20A를 회로에 공급하면 과전류로 인해서 회로에 열이 과다하게 발생하거나 혹은 손상이 가지 않는지요?

그리고 스텝모터 구동드라이버에도 12V, 17A면 전류 과다로 모터나 회로에 이상이 발생하지 않는가요?

모터 구동하는데 전류가 중요한데 보통 3A면 되는데 17A까지 공급되어도 되나요?

:전류는 필요한 만큼만 사용됩니다. 공급하는 측에서 남아돈다고 해서 무조건 다 흘러 들어가는 것이 아닙니다.

전류가 모자라면 문제가 됩니다만 남는 것은 사용되지 않습니다.

----------------------

스테핑 모터 스펙이 5V라도 정격으로 주어지는 전압을 사용하는게 좋다. 모터 컨트롤러의 정격을 알아보고 거기에 맞는 전압을 넣어주자.  적거나 무부하시 10 또는 12V라도 잘 돌아가는 이유는 스텝모터 구동부 내부에서 24V를 생성해 주기 때문이지만, 이런 전압변환에서는 전압이 낮아지는 대신 내부에 흐르는 전류가 높아지기 때문에 과열이 날 수 있다. 따라서 24V를 공급하는 것이 좋다. W=IV인데 I를 맞추고 V가 낮으면 돌아갈수는 있어도 내부전류가 높아져서 과열발생, 전류가 충분히 공급되지 않으면 전압이 불안정해지면서 구동부에 열 발생. 방열판 추가 장착이 필요할 수 있음.

모터 구동부로 가는 선이 너무 얇으면 충분한 전류가 전달될 수 없어서 모터 구동부에서 열이 발생

결론:

주어진 전류값은 지켜준다.

전압값은 4~8배가 될 수 있도록.

정격전압 24V정도로 맞춰주자.

모터 구동부로 가기 위한 선은 굵게 해주자

앞서했던 STFT에 Mel_Frequency Filter Bank를 적용시키는 부분입니다.

Log Operation 까지 적용했습니다.

%% Adapt mel_filterbank

filter_demension=30;

MFB=melfilterbank(filter_demension, N, FS);

%Frequency to Mel-frequency

%주파수->멜주파수

Aspec=MFB*PSPEC;%Aspec:Audio Spectrogram

[mel_x,mel_y]=meshgrid(time,[1:filter_demension]);

figure;

mesh(mel_x,mel_y,Aspec);

%% Get Log Operation

LAspec=log(Aspec);

figure;

mesh(mel_x,mel_y,LAspec);

필터의 갯수를 30개로 하였습니다.

MFB는 멜 필터뱅크로 인수는 (필터갯수, NFFT(주파수 레졸루션, 2^9), FS(주파수 16000Hz))

Aspec는 오디오스펙트로그램으로 앞서 구한 MFB에 파워스펙트로그램인 PSPEC를 곱해줍니다.

AudioSpectrogram

LAspec는 로그를 적용한 오디오스펙트로그램입니다.

Loged AudioSpectrogram

앞서 spectrogram함수로 구한 스펙트로그램과 비슷하다는걸 알 수 있습니다. 정보의 용량이 대폭 줄었죠.

spectrogram

오늘 다루려고 하는것은 위의 melfilterbank(filter_demension, N, FS) 함수입니다.

사실 이 부분을 제대로 이해하지 못해서, 다른 matlab코드를 분석하면서 공부했죠...

함수 내용입니다.

-------------------------------------------------------------------------------

function m = melfilterbank(p, n, fs)

% The function computed the mel filter banks for robust speaker recognition

% The filter spectrum is such that the passband area remains the same yet

% the pasband frequencies decrease and the power increases, to emphasis the

% higher frequency components

% p number of filters in filterbank

% n length of fft

% fs sample rate in Hz

f0 = 700/fs;

fn2 = floor(n/2);                                   %16000/2=8000

lr = log(1 + 0.5/f0) / (p+1);                       %p+1만큼 등간격 등분

% convert to fft bin numbers with 0 for DC term

bl = n * (f0 * (exp([0 1 p p+1] * lr) - 1));

b1 = floor(bl(1)) + 1;

b2 = ceil(bl(2));

b3 = floor(bl(3));

b4 = min(fn2, ceil(bl(4))) - 1;

pf = log(1 + (b1:b4)/n/f0) / lr;

fp = floor(pf);

pm = pf - fp;

r = [fp(b2:b4) 1+fp(1:b3)];

c = [b2:b4 1:b3] + 1;

v = [1-pm(b2:b4) pm(1:b3)];

m= sparse(r,c, v, p, 1+fn2);

figure;

hold on

for num=1:p

    plot(m(num,:));

end

-------------------------------------------------------------------------------

사실 위 함수 말고도 다른 함수가 여럿 있었습니다만  이 함수가 제일 메모리를 적게 먹더군요.

그 이유가 바로 sparse함수 때문이었습니다.

매트릭스를 보면 0이 많은 경우가 많습니다. 사실 이런 데이터는 굳이 표현해 줄 필요가 없습니다.

예를들어 290*380의 매트릭스에서 60퍼센트 정도의 데이터가 0이면 이걸 다 표현해주면 290*380*0.6이라는 엄청난 메모리 낭비가 있습니다. 여기서 sparse를 이용해서 유효한 데이터의 값과 좌표데이터만 표현해 주자는 것이죠.

안드로이드 앱등의 플랫폼에 옮긴다고 할 때 메모리 문제는 굉장히 예민하게 다루어집니다. 이런 부분에서 세심한 배려가 제품의 가치를 높이는 것이겠죠.

그럼 함수 설명을 시작하겠습니다. 하나하나 짚어가면서 진행합니다.

일단 멜 필터뱅크 구축에 사용되는 함수는 두가지 입니다.

M(f)는 주파수를 Mel주파수로 바꿉니다.

M^(-1)(m)은 Mel주파수를 주파수로 바꿔줍니다.

여기서 주의, 1125라는 값을 보죠. M^(-1)(m)의 m에 M(f)를 대입시키면 1125/1125가 되어 f가 됩니다.

위 함수에서 1125값이 보이지 않는것은 이것 때문입니다.

음성 분석을 할 때, 주어지는 FS의 절반에 해당하는 주파수의 정보까지 얻을 수 있습니다.

FS=16000, 따라서 8000Hz까지 분석하게 되죠. 앞서 STFT에서 8000Hz까지 분석하는것도 같은 이유입니다.

lr 에서 /(p+1)을 제외하고 봅시다.

lr = log(1 + 0.5/f0);

f0 = 700/fs; 까지 대입시켜 보면

lr = log(1+0.5*fs/700) 으로 (1)식에서 1125를 제외한 식이 됩니다.

여기서 /(p+1)을 하는 이유는 (p+1)등분 하기 위해서이죠. 30개의 필터를 만든다고 하면 p+1=31로 나눕니다.

앞서 영문자료 보면 10개 필터를 만들경우 포인트를 +2한 12포인터를 잡습니다. 일단 이렇게 알아둡시다.

bl을 분석해 봅시다.

bl = n * (f0 * (exp([0 1 p p+1] * lr) - 1));

bl = n/FS * (700*(exp(m)-1))

exp([0 1 p p+1])의 0에서 31은 32의 길이를 갖습니다. 30+2이죠?

위 식은 n/FS를 제외하면 (2)식에 1125를 뺀것과 같습니다. 

그럼 어째서 n/FS를 하느냐. 영문사이트 및 위에서 보는 그래프에 따르면 aplitude가 0~1로 표현되어 있습니다. 이와 같은 스케일로 표현하기 위해서 취해주었습니다.

실제로 bl의 결과가 0~256까지의 값이 나옵니다. 자연수로 해서 1단위로 보면 257의 길이를 갖습니다.

앞서 STFT에서 주파수 영역(0~8000)이 257개로 나뉜것과 같습니다.(STFT의 1+R/2=1+(2^9)/2과 같습니다.)

실제로 계산해보면 n*f/FS이기 때문에 (2^9)*(1/2)=256이 됩니다.

0 1 p p+1은 각각 용도가 있어서 넣은 값입니다.

b1 = floor(bl(1)) + 1;

b2 = ceil(bl(2));

b3 = floor(bl(3));

b4 = min(fn2, ceil(bl(4))) - 1;

bl(1)~bl(4)은 0 1 p p+1의 결과값입니다.

b1=1, b2=2, b3=234, b4=255가 됩니다.

----------------------------------------------------------------------여기까지 1단계로 봅시다.

pf를 식으로 표현하면 ln(1+ (b1:b4)*FS/(2^9)*1/700)/lr

pf가 가지는 값은 0에서 31 미만입니다.

fp는 버림을 적용함으로 0에서 30까지의 값을 가집니다.

pm은 pf에서 정수부를 제거한 값입니다. 이것이 가지는 의의는 영문 자료의 아래의 식을 구현했다는 것입니다.

이는 0Hz에서 8000Hz까지를 257등분하여 0에서 256까지의 257단계의 자연수 스케일로 표시하고 각 필터(0에서 30)마다 0에서 257까지의 스케일에 대하여 갖는 결과값(Hm, 0에서 1사이)을 갖는다. 이는 3차원 공간에의 자료이다.

Matlab으로 영문자료와 같은 그래프를 출력하려고 할때 그냥 plot(데이터)를 해버리면

이런 결과가 나온다. 

실제 3차원으로 보는 필터는 다음과 같다.

X,Z평면으로 보면 다음과 같다.

2차원 plot으로 hold on과 for을 이용하여 출력하면 다음과 같다.

필터를 몇개 제거하여보자.

for i=0:3:30으로 필터를 10개정도 제거하여 보자

필터를 통과한 성분들만 남아있는 것을 알 수 있다.

Mel_FilterBank에 사용된 식이 어떻게 나왔는가, Matlab 코드에서 좌표와 소수점 이하 자료로 나누는 부분 등 다루고 싶은 부분이 많으나 아직 이해가 부족하고 시간이 많이 소모되므로 프로젝트가 끝난 후로 미루도록 하겠다.

Mel Frequency Cepstral Coefficient (MFCC) tutorial

The first step in any automatic speech recognition system is to extract features i.e. identify the components of the audio signal that are good for identifying the linguistic content and discarding all the other stuff which carries information like background noise, emotion etc.

The main point to understand about speech is that the sounds generated by a human are filtered by the shape of the vocal tract including tongue, teeth etc. This shape determines what sound comes out. If we can determine the shape accurately, this should give us an accurate representation of the phoneme being produced. The shape of the vocal tract manifests itself in the envelope of the short time power spectrum, and the job of MFCCs is to accurately represent this envelope. This page will provide a short tutorial on MFCCs.

Mel Frequency Cepstral Coefficents (MFCCs) are a feature widely used in automatic speech and speaker recognition. They were introduced by Davis and Mermelstein in the 1980's, and have been state-of-the-art ever since. Prior to the introduction of MFCCs, Linear Prediction Coefficients (LPCs) and Linear Prediction Cepstral Coefficients (LPCCs) and were the main feature type for automatic speech recognition (ASR). This page will go over the main aspects of MFCCs, why they make a good feature for ASR, and how to implement them.

Steps at a Glance 

We will give a high level intro to the implementation steps, then go in depth why we do the things we do. Towards the end we will go into a more detailed description of how to calculate MFCCs.

1. Frame the signal into short frames.
2. For each frame calculate the periodogram estimate of the power spectrum.
3. Apply the mel filterbank to the power spectra, sum the energy in each filter.
4. Take the logarithm of all filterbank energies.
5. Take the DCT of the log filterbank energies.
6. Keep DCT coefficients 2-13, discard the rest.

There are a few more things commonly done, sometimes the frame energy is appended to each feature vector. Delta and Delta-Delta features are usually also appended. Liftering is also commonly applied to the final features.

Why do we do these things? 

We will now go a little more slowly through the steps and explain why each of the steps is necessary.

An audio signal is constantly changing, so to simplify things we assume that on short time scales the audio signal doesn't change much (when we say it doesn't change, we mean statistically i.e. statistically stationary, obviously the samples are constantly changing on even short time scales). This is why we frame the signal into 20-40ms frames. If the frame is much shorter we don't have enough samples to get a reliable spectral estimate, if it is longer the signal changes too much throughout the frame.

The next step is to calculate the power spectrum of each frame. This is motivated by the human cochlea (an organ in the ear) which vibrates at different spots depending on the frequency of the incoming sounds. Depending on the location in the cochlea that vibrates (which wobbles small hairs), different nerves fire informing the brain that certain frequencies are present. Our periodogram estimate performs a similar job for us, identifying which frequencies are present in the frame.

The periodogram spectral estimate still contains a lot of information not required for Automatic Speech Recognition (ASR). In particular the cochlea can not discern the difference between two closely spaced frequencies. This effect becomes more pronounced as the frequencies increase. For this reason we take clumps of periodogram bins and sum them up to get an idea of how much energy exists in various frequency regions. This is performed by our Mel filterbank: the first filter is very narrow and gives an indication of how much energy exists near 0 Hertz. As the frequencies get higher our filters get wider as we become less concerned about variations. We are only interested in roughly how much energy occurs at each spot. The Mel scale tells us exactly how to space our filterbanks and how wide to make them. Seebelow for how to calculate the spacing.

Once we have the filterbank energies, we take the logarithm of them. This is also motivated by human hearing: we don't hear loudness on a linear scale. Generally to double the percieved volume of a sound we need to put 8 times as much energy into it. This means that large variations in energy may not sound all that different if the sound is loud to begin with. This compression operation makes our features match more closely what humans actually hear. Why the logarithm and not a cube root? The logarithm allows us to use cepstral mean subtraction, which is a channel normalisation technique.

The final step is to compute the DCT of the log filterbank energies. There are 2 main reasons this is performed. Because our filterbanks are all overlapping, the filterbank energies are quite correlated with each other. The DCT decorrelates the energies which means diagonal covariance matrices can be used to model the features in e.g. a HMM classifier. But notice that only 12 of the 26 DCT coefficients are kept. This is because the higher DCT coefficients represent fast changes in the filterbank energies and it turns out that these fast changes actually degrade ASR performance, so we get a small improvement by dropping them.

What is the Mel scale? 

The Mel scale relates perceived frequency, or pitch, of a pure tone to its actual measured frequency. Humans are much better at discerning small changes in pitch at low frequencies than they are at high frequencies. Incorporating this scale makes our features match more closely what humans hear.

The formula for converting from frequency to Mel scale is:

To go from Mels back to frequency:

Implementation steps 

We start with a speech signal, we'll assume sampled at 16kHz.

1. Frame the signal into 20-40 ms frames. 25ms is standard. This means the frame length for a 16kHz signal is 0.025*16000 = 400 samples. Frame step is usually something like 10ms (160 samples), which allows some overlap to the frames. The first 400 sample frame starts at sample 0, the next 400 sample frame starts at sample 160 etc. until the end of the speech file is reached. If the speech file does not divide into an even number of frames, pad it with zeros so that it does.

The next steps are applied to every single frame, one set of 12 MFCC coefficients is extracted for each frame. A short aside on notation: we call our time domain signal . Once it is framed we have where n ranges over 1-400 (if our frames are 400 samples) and  ranges over the number of frames. When we calculate the complex DFT, we get  - where the  denotes the frame number corresponding to the time-domain frame.  is then the power spectrum of frame .

2. To take the Discrete Fourier Transform of the frame, perform the following:

where  is an  sample long analysis window (e.g. hamming window), and  is the length of the DFT. The periodogram-based power spectral estimate for the speech frame  is given by:

This is called the Periodogram estimate of the power spectrum. We take the absolute value of the complex fourier transform, and square the result. We would generally perform a 512 point FFT and keep only the first 257 coefficents.

3. Compute the Mel-spaced filterbank. This is a set of 20-40 (26 is standard) triangular filters that we apply to the periodogram power spectral estimate from step 2. Our filterbank comes in the form of 26 vectors of length 257 (assuming the FFT settings fom step 2). Each vector is mostly zeros, but is non-zero for a certain section of the spectrum. To calculate filterbank energies we multiply each filterbank with the power spectrum, then add up the coefficents. Once this is performed we are left with 26 numbers that give us an indication of how much energy was in each filterbank. For a detailed explanation of how to calculate the filterbanks see below. Here is a plot to hopefully clear things up:

4. Take the log of each of the 26 energies from step 3. This leaves us with 26 log filterbank energies.

5. Take the Discrete Cosine Transform (DCT) of the 26 log filterbank energies to give 26 cepstral coefficents. For ASR, only the lower 12-13 of the 26 coefficients are kept.

The resulting features (12 numbers for each frame) are called Mel Frequency Cepstral Coefficients.

Computing the Mel filterbank 

To get the filterbanks shown in figure 1(a) we first have to choose a lower and upper frequency. Good values are 300Hz for the lower and 8000Hz for the upper frequency. Of course if the speech is sampled at 8000Hz our upper frequency is limited to 4000Hz. Then follow these steps:

1. Using equation 1, convert the upper and lower frequencies to Mels. In our case 300Hz is 401.25 Mels and 8000Hz is 2834.99 Mels.
2. For this example we will do 10 filterbanks, for which we need 12 points. This means we need 10 additional points spaced linearly between 401.25 and 2834.99. This comes out to:

   m(i) = 401.25, 622.50, 843.75, 1065.00, 1286.25, 1507.50, 1728.74, 
          1949.99, 2171.24, 2392.49, 2613.74, 2834.99

3. Now use equation 2 to convert these back to Hertz:

   f(i) = 300, 517.33, 781.90, 1103.97, 1496.04, 1973.32, 2554.33, 
          3261.62, 4122.63, 5170.76, 6446.70, 8000

   Notice that our start- and end-points are at the frequencies we wanted.
4. Now we create our filterbanks. The first filterbank will start at the first point, reach its peak at the second point, then return to zero at the 3rd point. The second filterbank will start at the 2nd point, reach its max at the 3rd, then be zero at the 4th etc. A formula for calculating these is as follows:
   
   where  is the number of filters we want, and  is the list of M+2 Mel-spaced frequencies.

The final plot of all 10 filters overlayed on each other is:

Deltas and Delta-Deltas 

Also known as differential and acceleration coefficients. The MFCC feature vector describes only the power spectral envelope of a single frame, but it seems like speech would also have information in the dynamics i.e. what are the trajectories of the MFCC coefficients over time. It turns out that calculating the MFCC trajectories and appending them to the original feature vector increases ASR performance by quite a bit (if we have 12 MFCC coefficients, we would also get 12 delta coefficients, which would combine to give a feature vector of length 24).

To calculate the delta coefficients, the following formula is used:

where  is a delta coefficient, from frame  computed in terms of the static coefficients  to . A typical value for  is 2. Delta-Delta (Acceleration) coefficients are calculated in the same way, but they are calculated from the deltas, not the static coefficients.

References 

Davis, S. Mermelstein, P. (1980) Comparison of Parametric Representations for Monosyllabic Word Recognition in Continuously Spoken Sentences. In IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. 28 No. 4, pp. 357-366

X. Huang, A. Acero, and H. Hon. Spoken Language Processing: A guide to theory, algorithm, and system development. Prentice Hall, 2001.