语音信号处理matlab实现

短时能量分析matlab源程序：

x=wavread('4.wav');

%计算N=50，帧移=50时的语音能量

s=fra(50,50,x);%对输入的语音信号进行分帧，其中帧长50，帧移50

s2=s.^2;%一帧内各种点的能量

energy=sum(s2,2);%求一帧能量

subplot(2,2,1);

plot(energy)

xlabel('帧数');

ylabel('短时能量E');

legend('N=50');

axis([0,500,0,30])

%计算N=100，帧移=100时的语音能量

s=fra(100,100,x);

s2=s.^2;

energy=sum(s2,2);

subplot(2,2,2);

plot(energy)

xlabel('帧数');

ylabel('短时能量E');

legend('N=100');

axis([0,300,0,30])

%计算N=400，帧移=400时的语音能量

s=fra(400,400,x);

s2=s.^2;

energy=sum(s2,2);

subplot(2,2,3);

plot(energy)

xlabel('帧数');

ylabel('短时能量E');

legend('N=400');

axis([0,60,0,100])

%计算N=800，帧移=800时的语音能量

s=fra(800,800,x);

s2=s.^2;

energy=sum(s2,2);

subplot(2,2,4);

plot(energy)

xlabel('帧数');

ylabel('短时能量E');

legend('N=800');

axis([0,30,0,200])

分帧子函数：

function f=fra(len,inc,x)

%对读入语音分帧，len为帧长，inc为帧重叠样点数，x为输入语音数据

fh=fix(((size(x,1)-len)/inc)+1);%计算帧数

f=zeros(fh,len);%设一个零矩阵，行为帧数，列为帧长

i=1;n=1;

while i<=fh %帧间循环

j=1;

while j<=len %帧内循环

f(i,j)=x(n);

j=j+1;n=n+1;

end

n=n-len+inc;%下一帧开始位置

i=i+1;

end

运行结果：

短时自相关分析matlab源程序：

x=wavread('4.wav');

s1=x(1:5000);

N=320;%选择的窗长，加N=320的矩形窗

A=[];

for k=1:320

sum=0;

for m=1:N-(k-1)

sum=sum+s1(m)*s1(m+k-1);%计算自相关 end

A(k)=sum;

end

for k=1:320

A1(k)=A(k)/A(1);%归一化A（k）

end

N=160;

B=[];

for k=1:160

sum=0;

for m=1:N-(k-1)

sum=sum+s1(m)*s1(m+k-1);%计算自相关 end

B(k)=sum;

end

for k=1:160

B1(k)=B(k)/B(1);

end

N=70;

C=[];

for k=1:70

sum=0;

for m=1:N-(k-1)

sum=sum+s1(m)*s1(m+k-1);

end

C(k)=sum;

end

for k=1:70

C1(k)=C(k)/C(1);

end

figure(1)

subplot(3,1,1)

plot(A1)

xlabel('延时k')

ylabel('R(k)')

axis([0,350,-1,1]);

legend('N=320') subplot(3,1,2)

plot(B1)

xlabel('延时k')

ylabel('R(k)')

axis([0,350,-1,1]);

legend('N=160')

subplot(3,1,3)

plot(C1)

xlabel('延时k')

ylabel('R(k)')

axis([0,350,-1,1]);

legend('N=70')

修正的自相关matlab源程序：

b=wavread('4.wav');

b1=b(1:5000);

N=320;%选择的窗长，加N=320的矩形窗

A=[];

for k=1:320

sum=0;

for m=1:N

sum=sum+b1(m)*b1(m+k-1);%计算自相关 end

A(k)=sum;

end

for k=1:320

A1(k)=A(k)/A(1);%归一化A（k）

end

figure(1)

subplot(3,1,1)

plot(A1)

xlabel('延时k')

ylabel('R(k)')

axis([0,400,-0.5,1]);

legend('N=320')

b2=b(1:5000);

N=160;%选择的窗长，加N=160的矩形窗

B=[];

for k=1:160

sum=0;

for m=1:N

sum=sum+b2(m)*b2(m+k-1);%计算自相关 end

B(k)=sum;

end

for k=1:160

B1(k)=B(k)/B(1);%归一化A（k）end

figure(1)

subplot(3,1,2)

plot(B1)

xlabel('延时k')

ylabel('R(k)')

axis([0,400,-0.5,1]);

legend('N=160')

b3=b(1:2500);

N=70;%选择的窗长，加N=70的矩形窗C=[];

for k=1:70

sum=0; for m=1:N

sum=sum+b3(m)*b3(m+k-1);%计算自相关 end

C(k)=sum;

end

for k=1:70

C1(k)=C(k)/C(1);%归一化A（k）

end

figure(1)

subplot(3,1,3)

plot(C1)

xlabel('延时k')

ylabel('R(k)')

axis([0,400,-0.5,1]);

legend('N=70')

短时自相关运行结果：

修正自相关运行结果：

自相关法matlab源程序：

%本程序运行结果为中心削波前后的语音波形，以及削波前后的自相关波形

%读入数据采样fs=8khz 采样位数16bit

a=wavread('4.wav'); %打开语音文件

%读语音文件

L=length(a)%测定语音的长度

m=max(a)

for i=1:L

a(i)=a(i)/m;%数据归一化

end

m=max(a)%找到最大正值

n=min(a)%找到最小负值

ht=(m+n)/2;%为保证幅度值与横坐标对称

for i=1:L%数据中心下移，保持和横坐标轴对称

a(i)=a(i)-ht;

end

figure(1)

subplot(2,1,1)

plot(a)

axis([0,10000,-1,1]);

title('中心削波前语音波形')

xlabel('样点数')

ylabel('幅度')

coeff=0.7;%中心削波函数系数取0.7

th0=max(a)*coeff;%求中心削波函数门限

for k=1:L%中心削波

if a(k)>=th0

a(k)=a(k)-th0;

elseif a(k)<=(-th0)

a(k)=a(k)+th0;

else a(k)=0;

end

end

m=max(a);

for i=1:L%中心削波函数幅度的归一化

a(i)=a(i)/m;

end

subplot(212)

plot(a)

axis([0,10000,-2,2]);

title('中心削波后语音波形')

xlabel('样点数')

ylabel('幅度')

%没有经过中心削波的修正自相关计算b=wavread('4.wav');

N=320;%选择的窗长，加N=320的矩形窗

A=[];

for k=1:320%选择延迟长度

sum=0;

for m=1:N

sum=sum+b(m+7500)*b(m+7500+k-1);%计算自相关

end

A(k)=sum;

end

for k=1:320

B(k)=A(k)/A(1);%归一化A（k）

end

figure(2)

subplot(2,1,1)

plot(B)

axis([0,400,-1,1]);

title('中心削波前修正自相关')

xlabel('延时k')

ylabel('幅度')

%中心削波函数和修正的自相关方法结合

N=320;%选择的窗长，加N=320的矩形窗

A=[];

for k=1:320

sum=0;

for m=1:N

sum=sum+a(m+7500)*a(m+7500+k-1);%对削波后的函数计算自相关

end

A(k)=sum;

end

for k=1:320

C(k)=A(k)/A(1)%归一化A（k）

end

figure(2)

subplot(2,1,2)

plot(C)

axis([0,400,-1,1]);

title('中心削波后修正自相关')

xlabel('延时k')

ylabel('幅度')

运行结果：

AMDF法matlab源程序：

b=wavread('4.wav');

b1=b(3500:6000);

N=320;%选择的窗长，加N=320的矩形窗

A=[];

for k=1:320

sum=0;

for m=1:N

sum=sum+abs(b1(m)-b1(m+k-1));%计算自相关

end

A(k)=sum;

end s=b(3500:10000); figure(1)

subplot(211)

plot(s)

xlabel('样点') ylabel('幅度')

axis([0,2500,-1,1]); subplot(212)

plot(A)

xlabel('延时k') ylabel('AMDF')

axis([0,400,0,200]);

运行结果：

DFT谱估计matlab源程序：x=wavread('4.wav');

f=fra(256,50,x);

for i=1:size(f,1)

y=f(i,:);

end

s=y'.*hamming(256);

k=0:255;

kk=k/256;

sk=fft(s,256);

ssk=log(abs(sk));

subplot(211)

plot(s);

title('语音信号');

subplot(212)

plot(kk,ssk);

title('DFT谱');

运行结果：

LPC谱估计matlab源程序：

x=wavread('4.wav');

f=fra(256,50,x);

for i=1:size(f,1)

y=f(i,:);

end

s=y'.*hamming(256);

[a,g]=lpc(s,12);

[h,f]=freqz(g,a,100,8000); plot(f,20*log10(abs(h))); 运行结果：倒谱法估计matlab源程序：

y=wavread('4.wav');

f=fra(256,50,y);

for i=1:size(f,1)

x=f(i,:);

end

s=x'.*hamming(256);

sk=fft(s,256);

ssk=log(abs(sk));

cn=ifft(ssk,256); %求出语音信号倒谱c(n)

l(1)=1;

for i=2:63

l(i)=2;

end % 产生低时窗l( )

for i=64:256

l(i)=0;

end

y=cn'.*l %从f(n)中截取出h(n)

yk=fft(y,256);%求h(n)的DFT变换

k=0:255;

kk=k/256;

plot(kk,yk)

运行结果：

LPC 复倒谱估计matlab 源程序: y=wavread('4.wav'); f=fra(256,50,y); for i=1:size(f,1) x=f(i,:); end

s=x'.*hamming(256);

lpcc=lpc2lpcc(s);%调用子函数lpc21pcc 计算LPCC 系数

lpcc2=fliplr(lpcc);

cn=(lpcc+lpcc2)/2;%求出倒谱系数 lk=fft(cn,256);

llk=log(abs(lk));%求其频谱 k=0:255; kk=k/256; plot(kk,llk)

lpc2lpcc 子函数源程序：

function lpcc=lpc2lpcc(lpc) n_lpc=30; n_lpcc=40;

lpcc=zeros(n_lpcc,1); lpcc(1)=lpc(1); for n= 2:n_lpc

lpcc(n)=lpc(n);

for l= 1:n-1

lpcc(n)=lpcc(n)+lpc(1)*lpcc(n-1)*(n-1)/n;

end end

for n= n_lpc+ 1:n_lpcc lpcc(n)=0; for l=1:n_lpc

lpcc(n)=lpcc(n)+lpc(l)*lpcc(n-1)*(n-1)/n;

end end

lpcc= -lpcc;

运行结果：