语言识别的matlab仿真程序

 

matlab代写程序   网址：http://item.taobao.com/item.htm?id=17898804306  加QQ：1046233911

matlab代写程序   网址：http://item.taobao.com/item.htm?id=17898804306  加QQ：1046233911

一、加窗和分帧函数

% enframe.m     加窗和分帧分帧函数，用在语音处理中的。相当于加了一个移动的窗

%    一般调用格式为

%     y=enframe(x,framelength,step);

nx=length(x);

nwin=length(win);

if (nwin == 1)

   len = win;

else

   len = nwin;

end

if (nargin < 3)

   inc = len;

end

nf = fix((nx-len+inc)/inc);

f=zeros(nf,len);

indf= inc*(0:(nf-1)).';

inds = (1:len);

f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));

if (nwin > 1)

    w = win(:)';

    f = f .* w(ones(nf,1),:);

end

二、语音信号端点检测程序

function [x1,x2]=vad(k,fs,i)

%% [X1,X2]=vad(k,fs)

%%     语音信号端点检测程序，k为语音信号，

%%     fs为其采样频率，程序绘制出语音信号

%%     相关波形图并返回起止端点所对帧号

%close all;

% 幅度归一化到[-1,1]

k=double(k);

k=k./max(abs(k));     %   k>0,k=1;k<=-1.

%------------------------------

%

此处录取部分代码

(3,1,3);

plot(t1,k1);

axis([0.4,0.5,min(k),max(k)]);

z=sprintf('(III)    “�.wav”语音结束处放大波形图',i);

title(z);

xlabel('Time:s');

ylabel('Amplitude(normalized)');

%------------------------------

%        计算短时过零率

%------------------------------

FrameLen=240;

FrameInc=80;

FrameTemp1=enframe(k(1:end-1),FrameLen,FrameInc);

FrameTemp2=enframe(k(2:end),FrameLen,FrameInc);

signs=(FrameTemp1.*FrameTemp2)<0;

diffs=abs(FrameTemp1-FrameTemp2)>0.01;

zcr=sum(signs.*diffs,2);

disp('显示原始波形图……');

figure,subplot(3,1,1);

plot(t,k);

axis([0,(length(k)-1)/fs,min(k),max(k)]);

z=sprintf('(I)  “�.wav”语音信号波形',i);

title(z);

xlabel('Time:s');

ylabel('Amplitude(normalized)');

disp('显示过短时零率……')

zcrInd=1:length(zcr);

subplot(3,1,2);

plot(zcrInd,zcr);

axis([0,length(zcr),0,max(zcr)]);

title('(II)  短时过零率');

xlabel('Frame');

ylabel('Zcr');

%------------------------------

%         计算短时能量

%------------------------------

amp=sum(abs(enframe(filter([1 -0.9375], 1, k), FrameLen, FrameInc)), 2);

disp('显示短时能量……')

ampInd=1:length(amp);

subplot(3,1,3);

plot(ampInd,amp);

axis([0,length(amp),0,max(amp)]);

title('(III)  短时能量');

xlabel('Frame');

ylabel('Energy');

%------------------------------

%         设置门限

%------------------------------

disp('设置门限……');

ZcrLow=max([round(mean(zcr)*0.1),3]);                   %过零率低门限

ZcrHigh=max([round(max(zcr)*0.1),5]);                   %过零率高门限

AmpLow=min([min(amp)*10,mean(amp)*0.2,max(amp)*0.1]);   %能量低门限

AmpHigh=max([min(amp)*10,mean(amp)*0.2,max(amp)*0.1]);  %能量高门限

%------------------------------

%         端点检测

%------------------------------

MaxSilence=8;   %最长语音间隙时间

MinAudio=16;    %最短语音时间

Status=0;       %状态：0静音段,1过渡段,2语音段,3结束段

HoldTime=0;     %语音持续时间

SilenceTime=0;  %语音间隙时间

disp('开始端点检测……');

for n=1:length(zcr)

    switch Status

        case{0,1}

            if amp(n)>AmpHigh || zcr(n)>ZcrHigh

               x1=n-HoldTime;

                Status=2;

                HoldTime=HoldTime+1;

                SilenceTime=0;

            elseif amp(n)>AmpLow || zcr(n)>ZcrLow

                Status=1;

                HoldTime=HoldTime+1;

 

此处略去部分代码     

   HoldTime=HoldTime+1;

            else

               SilenceTime=SilenceTime+1;

               if SilenceTime<MaxSilence

                   HoldTime=HoldTime+1;

               elseif (HoldTime-SilenceTime)<MinAudio

                   Status=0;

                   HoldTime=0;

                   SilenceTime=0;

               else

                   Status=3;

               end

            end

        case 3,

            break;

    end

    if Status==3

        break;

    end

end

HoldTime=HoldTime-SilenceTime;

x2=x1+HoldTime;

disp('显示端点……');

figure,subplot(3,1,1);

plot(k);

axis([1,length(k),min(k),max(k)]);

z=sprintf('(I)  “�.wav”语音信号',i);

title(z);

xlabel('Sample');

ylabel('Speech');

line([x1*FrameInc,x1*FrameInc],[min(k),max(k)],'color','red');

line([x2*FrameInc,x2*FrameInc],[min(k),max(k)],'color','red');

subplot(3,1,2);

plot(zcr);

axis([1,length(zcr),0,max(zcr)]);

title('(II)  短时过零率');

xlabel('Frame');

ylabel('ZCR');

line([x1,x1],[0,max(zcr)],'Color','red');

line([x2,x2],[0,max(zcr)],'Color','red');

subplot(3,1,3);

plot(amp);

axis([1,length(amp),0,max(amp)]);

title('(III)  短时能量');

xlabel('Frame');

ylabel('Energy');

line([x1,x1],[0,max(amp)],'Color','red');

line([x2,x2],[0,max(amp)],'Color','red');

三、计算MFCC系数

function cc=mfcc(k)

%------------------------------

% cc=mfcc(k)计算语音k的MFCC系数

%------------------------------

% M为滤波器个数，N为一帧语音采样点数

M=24; N=256;

% 归一化mel滤波器组系数

bank=melbankm(M,N,8000,0,0.5,'m');

figure;

plot(linspace(0,N/2,129),bank);

title('Mel-Spaced Filterbank');

xlabel('Frequency [Hz]');

bank=full(bank);

bank=bank/max(bank(:));

% DCT系数,12*24

for i=1:12

  j=0:23;

  dctcoef(i,:)=cos((2*j+1)*i*pi/(2*24));

end

% 归一化倒谱提升窗口

w=1+6*sin(pi*[1:12]./12);

w=w/max(w);

% 预加重

AggrK=double(k);

AggrK=filter([1,-0.9375],1,AggrK);

% 分帧

FrameK=enframe(AggrK,N,80);

% 加窗

for i=1:size(FrameK,1)

    FrameK(i,:)=(FrameK(i,:))'.*hamming(N);

end

FrameK=FrameK';

% 计算功率谱

S=(abs(fft(FrameK))).^2;

disp('显示功率谱……')

figure;

plot(S);

axis([1,size(S,1),0,2]);

title('Power Spectrum (M=24, N=256)');

xlabel('Frame');

ylabel('Frequency [Hz]');

colorbar;

% 将功率谱通过滤波器组

P=bank*S(1:129,:);

% 取对数后作离散余弦变换

D=dctcoef*log(P);

% 倒谱提升窗

for i=1:size(D,2)

    m(i,:)=(D(:,i).*w')';

end

% 差分系数

dtm=zeros(size(m));

for i=3:size(m,1)-2

  dtm(i,:)=-2*m(i-2,:)-m(i-1,:)+m(i+1,:)+2*m(i+2,:);

end

dtm=dtm/3;

%合并mfcc参数和一阶差分mfcc参数

cc=[m,dtm];

%去除首尾两帧，因为这两帧的一阶差分参数为0

cc=cc(3:size(m,1)-2,:);

四、计算MFCC时所需Mel尺度滤波器组

function [x,mn,mx]=melbankm(p,n,fs,fl,fh,w)

%MELBANKM determine matrix for a mel-spaced filterbank [X,MN,MX]=(P,N,FS,FL,FH,W)

if nargin < 6

  w='tz';

  if nargin < 5

    fh=0.5;

    if nargin < 4

      fl=0;

    end

  end

end

f0=700/fs;

fn2=floor(n/2);

lr=log((f0+fh)/(f0+fl))/(p+1);

% convert to fft bin numbers with 0 for DC term

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

b2=ceil(bl(2));

 

此处略去部分代码

 

  v=2*[0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5];

  mn=1;

  mx=fn2+1;

else

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

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

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

  fp=floor(pf);

  pm=pf-fp;

  k2=b2-b1+1;

  k3=b3-b1+1;

  k4=b4-b1+1;

  r=[fp(k2:k4) 1+fp(1:k3)];

  c=[k2:k4 1:k3];

  v=2*[1-pm(k2:k4) pm(1:k3)];

  mn=b1+1;

  mx=b4+1;

end

if any(w=='n')

  v=1-cos(v*pi/2);

elseif any(w=='m')

  v=1-0.92/1.08*cos(v*pi/2);

end

if nargout > 1

  x=sparse(r,c,v);

else

  x=sparse(r,c+mn-1,v,p,1+fn2);

end

五、DTW算法

function dist = dtw(test, ref)

global x y_min y_max

global t r

global D d

global m n

t = test;

r = ref;

n = size(t,1);

m = size(r,1);

 

d = zeros(m,1);

D =  ones(m,1) * realmax;

D(1) = 0;

% 如果两个模板长度相差过多，匹配失败

if (2*m-n<3) | (2*n-m<2)

         dist = realmax;

         return

end

 

% 计算匹配区域

xa = round((2*m-n)/3);

xb = round((2*n-m)*2/3);

 

if xb>xa

         %xb>xa, 按下面三个区域匹配

         %        1   :xa

         %        xa+1:xb

         %        xb+1:N

         for x = 1:xa

                   y_max = 2*x;

                   y_min = round(0.5*x);

                   warp

         end

         for x = (xa+1):xb

                   y_max = round(0.5*(x-n)+m);

                   y_min = round(0.5*x);

                   warp

         end

         for x = (xb+1):n

                   y_max = round(0.5*(x-n)+m);

                   y_min = round(2*(x-n)+m);

                   warp

         end

elseif xa>xb

         %xa>xb, 按下面三个区域匹配

         %        0   :xb

         %        xb+1:xa

         %        xa+1:N

         for x = 1:xb

                   y_max = 2*x;

                   y_min = round(0.5*x);

                   warp

         end

         for x = (xb+1):xa

                   y_max = 2*x;

                   y_min = round(2*(x-n)+m);

                   warp

         end

         for x = (xa+1):n

                   y_max = round(0.5*(x-n)+m);

                   y_min = round(2*(x-n)+m);

                   warp

         end

elseif xa==xb

         %xa=xb, 按下面两个区域匹配

         %        0   :xa

         %        xa+1:N

         for x = 1:xa

                   y_max = 2*x;

                   y_min = round(0.5*x);

                   warp

         end

         for x = (xa+1):n

                   y_max = round(0.5*(x-n)+m);

                   y_min = round(2*(x-n)+m);

                   warp

         end

end

 

%返回匹配分数

dist = D(m);

 

function warp

global x y_min y_max

global t r

global D d

global m n

 

d = D;

for y = y_min:y_max

         D1 = D(y);

         if y>1

                   D2 = D(y-1);

         else

        D2 = realmax;

         end

         if y>2

                   D3 = D(y-2);

         else

        D3 = realmax;

         end

    d(y) = sum((t(x,:)-r(y,:)).^2) + min([D1,D2,D3]);

end

D = d;

六、语音训练阶段程序

%% train.m

disp('正在生成训练参数……');

for i=0:3

    fname=sprintf('train\\�.wav',i );

    [k,fs]=wavread(fname);

    [x1,x2]=vad(k,fs,i);

    cc=mfcc(k);

    cc=cc(x1-2:x2-2,:);

    ref(i+1).mfcc=cc;

end

disp('正在存储模板库……');

save 'mfcc.mat' ref;

%   close all;

clc;

七、基于DTW的语音识别程序

对test中所有语音进行识别

% test1.m

clear;

close all;

clc;

disp('正在导入参考模板参数...');

load mfcc.mat;

disp('正在计算测试模板的参数...')

  for i=0:3

         fname = sprintf('test\\�.wav',i);

    [k,fs]=wavread(fname);

    [x1,x2]=vad(k,fs,i);

    cc=mfcc(k);

    cc=cc(x1-2:x2-2,:);

    test(i+1).mfcc=cc;

 

此处略去部分代码

        dist(i,j) = dtw(test(i).mfcc, ref(j).mfcc);

    end

 end

disp('正在计算匹配结果...')

for i=1:4

         [d,j] = min(dist(i,:));  %  d为返回值，j为列数

   if j==1

       fprintf('测试模板 �.wav 的识别结果为:前进\n',i-1);end

 

       if j==2

       fprintf('测试模板 �.wav 的识别结果为:后退\n',i-1);   end

 

       if j==3

       fprintf('测试模板 �.wav 的识别结果为:向左\n',i-1);  end

 

        if j==4

       fprintf('测试模板 �.wav 的识别结果为:向右\n',i-1);  end 

end

  close all;

 

对test文件夹中任一个语音进行识别

% test2.m

clear;

close all;

clc;

disp('正在导入参考模板参数...');

load mfcc.mat;

disp('正在计算测试模板的参数...')

    i=1;

%     i=2;

%       i=3;

%    i=4;

         fname = sprintf('test\\�.wav',i-1);

    [k,fs]=wavread(fname);

    [x1,x2]=vad(k,fs,i-1);

    cc=mfcc(k);

    cc=cc(x1-2:x2-2,:);

    test(i).mfcc=cc;

disp('正在进行模板匹配...');

dist = zeros(1,4);

for j=1:4;

    dist(j)=dtw(test(i).mfcc,ref(j).mfcc);

end

disp('正在计算匹配结果...');    

m=min(dist);

[row,column]=find(dist==m);

j=column;

if j==1;        fprintf('测试模板 �.wav 的识别结果为:前进\n',i-1);  end

if j==2;       fprintf('测试模板 �.wav 的识别结果为:后退\n',i-1);   end

if j==3;       fprintf('测试模板 �.wav 的识别结果为:向左\n',i-1);   end

if j==4;      fprintf('测试模板 �.wav 的识别结果为:向右\n',i-1);   end

%  close all;