语音信号处理matlab实现

短时能量分析matlab源程序：x=wavread('4.wav');%计算N=50，帧移=50时的语音能量s=fra(50,50,x);%对输入的语音信号进行分帧，其中帧长50，帧移50s2=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])分帧子函数：functionf=fra(len,inc,x)%对读入语音分帧，len为帧长，inc为帧重叠样点数，x为输入语音数据fh=fix(((size(x,1)-len)/inc)+1);%计算帧数f=zeros(fh,len);%设一个零矩阵，行为帧数，列为帧长i=1;n=1;whilei=fh%帧间循环j=1;whilej=len%帧内循环f(i,j)=x(n);j=j+1;n=n+1;endn=n-len+inc;%下一帧开始位置i=i+1;end运行结果：短时自相关分析matlab源程序：x=wavread('4.wav');s1=x(1:5000);N=320;%选择的窗长，加N=320的矩形窗A=[];fork=1:320sum=0;form=1:N-(k-1)sum=sum+s1(m)*s1(m+k-1);%计算自相关endA(k)=sum;endfork=1:320A1(k)=A(k)/A(1);%归一化A（k）endN=160;B=[];fork=1:160sum=0;form=1:N-(k-1)sum=sum+s1(m)*s1(m+k-1);%计算自相关endB(k)=sum;endfork=1:160B1(k)=B(k)/B(1);endN=70;C=[];fork=1:70sum=0;form=1:N-(k-1)sum=sum+s1(m)*s1(m+k-1);endC(k)=sum;endfork=1:70C1(k)=C(k)/C(1);endfigure(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=[];fork=1:320sum=0;form=1:Nsum=sum+b1(m)*b1(m+k-1);%计算自相关endA(k)=sum;endfork=1:320A1(k)=A(k)/A(1);%归一化A（k）endfigure(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=[];fork=1:160sum=0;form=1:Nsum=sum+b2(m)*b2(m+k-1);%计算自相关endB(k)=sum;endfork=1:160B1(k)=B(k)/B(1);%归一化A（k）endfigure(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=[];fork=1:70sum=0;form=1:Nsum=sum+b3(m)*b3(m+k-1);%计算自相关endC(k)=sum;endfork=1:70C1(k)=C(k)/C(1);%归一化A（k）endfigure(1)subplot(3,1,3)plot(C1)xlabel('延时k')ylabel('R(k)')axis([0,400,-0.5,1]);legend('N=70')短时自相关运行结果：修正自相关运行结果：自相关法matlab源程序：%本程序运行结果为中心削波前后的语音波形，以及削波前后的自相关波形%读入数据采样fs=8khz采样位数16bita=wavread('4.wav');%打开语音文件%读语音文件L=length(a)%测定语音的长度m=max(a)fori=1:La(i)=a(i)/m;%数据归一化endm=max(a)%找到最大正值n=min(a)%找到最小负值ht=(m+n)/2;%为保证幅度值与横坐标对称fori=1:L%数据中心下移，保持和横坐标轴对称a(i)=a(i)-ht;endfigure(1)subplot(2,1,1)plot(a)axis([0,10000,-1,1]);title('中心削波前语音波形')xlabel('样点数')ylabel('幅度')coeff=0.7;%中心削波函数系数取0.7th0=max(a)*coeff;%求中心削波函数门限fork=1:L%中心削波ifa(k)=th0a(k)=a(k)-th0;elseifa(k)=(-th0)a(k)=a(k)+th0;elsea(k)=0;endendm=max(a);fori=1:L%中心削波函数幅度的归一化a(i)=a(i)/m;endsubplot(212)plot(a)axis([0,10000,-2,2]);title('中心削波后语音波形')xlabel('样点数')ylabel('幅度')%没有经过中心削波的修正自相关计算b=wavread('4.wav');N=320;%选择的窗长，加N=320的矩形窗A=[];fork=1:320%选择延迟长度sum=0;form=1:Nsum=sum+b(m+7500)*b(m+7500+k-1);%计算自相关endA(k)=sum;endfork=1:320B(k)=A(k)/A(1);%归一化A（k）endfigure(2)subplot(2,1,1)plot(B)axis([0,400,-1,1]);title('中心削波前修正自相关')xlabel('延时k')ylabel('幅度')%中心削波函数和修正的自相关方法结合N=320;%选择的窗长，加N=320的矩形窗A=[];fork=1:320sum=0;form=1:Nsum=sum+a(m+7500)*a(m+7500+k-1);%对削波后的函数计算自相关endA(k)=sum;endfork=1:320C(k)=A(k)/A(1)%归一化A（k）endfigure(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=[];fork=1:320sum=0;form=1:Nsum=sum+abs(b1(m)-b1(m+k-1));%计算自相关endA(k)=sum;ends=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);fori=1:size(f,1)y=f(i,:);ends=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);fori=1:size(f,1)y=f(i,:);ends=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);fori=1:size(f,1)x=f(i,:);ends=x'.*hamming(256);sk=fft(s,256);ssk=log(abs(sk));cn=ifft(ssk,256);%求出语音信号倒谱c(n)l(1)=1;fori=2:63l(i)=2;end%产生低时窗l()fori=64:256l(i)=0;endy=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);fori=1:size(f,1)x=f(i,:);ends=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子函数源程序：functionlpcc=lpc2lpcc(lpc)n_lpc=30;n_lpcc=40;lpcc=zeros(n_lpcc,1);lpcc(1)=lpc(1);forn=2:n_lpclpcc(n)=lpc(n);forl=1:n-1lpcc(n)=lpcc(n)+lpc(1)*lpcc(n-1)*(n-1)/n;endendforn=n_lpc+1:n_lpcclpcc(n)=0;forl=1:n_lpclpcc(n)=lpcc(n)+