语音信号处理实验

实验一特征提取1.1语谱图程序语音信号随时间变化的频谱特性可以用语谱图直观的表示，语谱图的纵坐标对应频率，横坐标对应时间，而图像的黑百度对应于信号的能量。所以，声道的谐振频率在图上就表示成为黑带，浊音部分则以出现条纹图形为其特征，这是因为此时的时域波形有周期性，而在浊音的时间间隔内图形显得很致密。clearall;[x,sr]=wavread('***.wav');%sr为采样频率if(size(x,1)size(x,2))%size(x,1)为x的行数，size(x,2)为x的列数x=x';ends=length(x);w=round(44*sr/1000);%窗长,取离44*sr/100最近的整数n=w;%fft的点数ov=w/2;%50%的重叠h=w-ov;%win=hanning(n)';%哈宁窗win=hamming(n)';%哈宁窗c=1;ncols=1+fix((s-n)/h);%fix函数是将(s-n)/h的小数舎去d=zeros((1+n/2),ncols);forb=0:h:(s-n)u=win.*x((b+1):(b+n));t=fft(u);d(:,c)=t(1:(1+n/2))';c=c+1;endtt=[0:h:(s-n)]/sr;ff=[0:(n/2)]*sr/n;imagesc(tt/1000,ff/1000,20*log10(abs(d)));colormap(gray);axisxyxlabel('时间/s');ylabel('频率/kHz');1.2预加重（高频提升）对输入的语音信号进行预加重，其目的是提升语音信号的高频部分，去除口唇辐射的影响，增强语音的高频分辨率。-1H(Z)=1-z[x,sr]=wavread('***.wav');%读数据ee=x(200:455);%选取原始文件e的第200到455点的语音,也可选其他样点r=fft(ee,1024);%对信号ee进行1024点傅立叶变换r1=abs(r);%对r取绝对值r1表示频谱的幅度值pinlv=(0:1:255)*8000/512%点和频率的对应关系yuanlai=20*log10(r1)%对幅值取对数signal(1:256)=yuanlai(1:256);%取256个点，目的是画图的时候，维数一致[h1,f1]=freqz([1,-0.98],[1],256,4000);%高通滤波器pha=angle(h1);%高通滤波器的相位H1=abs(h1);%高通滤波器的幅值r2(1:256)=r(1:256)u=r2.*h1'%将信号频域与高通滤波器频域相乘相当于在时域的卷积u2=abs(u)%取幅度绝对值u3=20*log10(u2)%对幅值取对数un=filter([1,-0.98],[1],ee)%un为经过高频提升后的时域信号figure(1);subplot(211);plot(f1,H1);title('高通滤波器的幅频响应');xlabel('频率/Hz');ylabel('幅度');subplot(212);plot(pha);title('高通滤波器的相位响应');xlabel('频率/Hz');ylabel('角度/radians');figure(2);subplot(211);plot(ee);title('原始语音信号');xlabel('样点数');ylabel('幅度');axis([0256-3*10^42*10^4]);subplot(212);plot(real(un));title('经高通滤波后的语音信号');xlabel('样点数');ylabel('幅度');axis([0256-1*10^41*10^4]);figure(3);subplot(211);plot(pinlv,ee);title('原始语音信号频谱');xlabel('频率/Hz');ylabel('幅度/dB');subplot(212);plot(pinlv,u3);title('经高通滤波后的语音信号频谱');xlabel('频率/Hz');ylabel('幅度/dB');1.3短时能量语音信号的能量随时间而变化，清音和浊音之间的能量差别相当显著，因此对短时能量进行分析，可以描述语音的这种特征变化的情况。2(1)[()()]nnmnNExmwnmn时刻语音信号的短时平均能量，其中N为窗长。[x,sr]=wavread('*.wav');%读入语音文件%计算N=50，帧移=50时的语音能量s=fra(50,50,x);s2=s.^2;%一帧内各样点的能量energy=sum(s2,2);%求一帧能量subplot(2,2,1)%定义画图数量和布局plot(energy);%画N=50时的语音能量图xlabel('帧数')%横坐标ylabel('短时能量E')%纵坐标legend('N=50')%曲线标识axis([0,1500,0,2*10^10])%定义横纵坐标范围%计算N=100，帧移=100时的语音能量s=fra(100,100,x);s2=s.^2;energy=sum(s2,2);subplot(2,2,2)plot(energy)%画N=100时的语音能量图xlabel('帧数')ylabel('短时能量E')legend('N=100')axis([0,750,0,4*10^10])%定义横纵坐标范围%计算N=400，帧移=400时的语音能量s=fra(400,400,x);s2=s.^2;energy=sum(s2,2);subplot(2,2,3)plot(energy)%画N=400时的语音能量图xlabel('帧数')ylabel('短时能量E')legend('N=400')axis([0,190,0,1.5*10^11])%定义横纵坐标范围%计算N=800，帧移=800时的语音能量s=fra(800,800,x);s2=s.^2;energy=sum(s2,2);subplot(2,2,4)plot(energy)%画N=800时的语音能量图xlabel('帧数')ylabel('短时能量E')legend('N=800')axis([0,95,0,3*10^11])%定义横纵坐标范围分帧函数：fra().m的定义functionf=fra(len,inc,x)fh=fix(((size(x,1)-len)/inc)+1);f=zeros(fh,len);i=1;n=1;whilei=fhj=1;whilej=lenf(i,j)=x(n);j=j+1;n=n+1;endn=n-len+inc;i=i+1;end1.4短时平均过零率短时平均过零率是指每帧语音信号通过零值电平的次数。在离散域，如果相邻的采样值具有不同符号就称为发生了过零，因此可以计算过零的次数。(1)sgn[()]sgn[(1)]()nnmnNZxmxmwnmclearall[x1,sr]=wavread('*.wav');%读入语音文件x=awgn(x1,15,'measured');%加入15dB的噪声s=fra(220,110,x);%分帧，帧移110zcr=zcro(s);%求过零率figure(1);subplot(2,1,1)plot(x);title('原始信号');xlabel('样点数');ylabel('幅度');axis([0,39760,-2*10^4,2*10^4]);subplot(2,1,2)plot(zcr);xlabel('帧数');ylabel('过零次数');title('原始信号的过零率');axis([0,360,0,200]);过零率函数：zcro.mfunctionf=zcro(x)f=zeros(size(x,1),1);%生成全零矩阵fori=1:size(x,1)z=x(i,:);%提取一行数据forj=1:(length(z)-1);ifz(j)*z(j+1)0;f(i)=f(i)+1;endendend实验二基音周期估计基音周期是表征语音信号本质特征的参数，在语音编码、语音合成和语音识别中有重要的作用。语音信号基音周期的估计方法很多，本实验采用两种最基本的方法：基于短时自相关法和基于短时平均幅度差法。2.1基于短时自相关法的基音周期估计自相关函数的性质是：如果是一个周期信号，那这个信号的自相关函数也是同周期的信号，且在信号周期的整数倍处，自相关值有最大值。语音信号的浊音信号具有准周期特性，其自相关函数在基音周期的整数倍处取最大值，计算两相邻最大峰值之间的距离，就可以估计出基音周期。为了减小运算量，突出反映基因周期的信息，压缩无关信息，可以对语音信号进行中心削波处理，再进行自相关计算。中心削波[a,count]=wavread('*.wav');L=length(a);%测定语音的长度m=max(a);fori=1:La(i)=a(i)/m;%数据归一化end%找到归一化以后数据的最大值和最小值m=max(a);%找到最大的正值n=min(a);%找到最小的负值%为保证幅度值与横坐标轴对称，采用计算公式是n+(m-n)/2,合并为(m+n)/2ht=(m+n)/2;fori=1:L;%数据中心下移保持和横坐标轴对称a(i)=a(i)-ht;endfigure(1);%画第一幅图subplot(2,1,1);%第一个子图plot(a,'k');axis([0,1711,-1,1]);%确定横、纵坐标的范围title('中心削波前语音波形');%图标题xlabel('样点数');%横坐标ylabel('幅度');%纵坐标coeff=0.7;%中心削波函数系数取0.7th0=max(a)*coeff;%求中心削波函数门限(threshold)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(2,1,2);%第二个子图plot(a,'k');axis([0,1711,-1,1]);%确定横、纵坐标的范围title('中心削波后语音波形');%图标题xlabel('样点数');%横坐标ylabel('幅度');%纵坐标%没有经过中心削波的修正自相关计算[b,count]=wavread('*.wav')；N=320;%选择的窗长A=[];fork=1:320;%选择延迟长度sum=0;form=1:N;sum=sum+b(m)*b(m+k-1);%计算自相关endA(k)=sum;endfork=1:320B(k)=A(k)/A(1);%自相关归一化endfigure(2);%画第二幅图subplot(2,1,1);%第一个子图plot(B,'k');title('中心削波前修正自相关');%图标题xlabel('延时k');%横坐标ylabel('幅度');%纵坐标axis([0,320,-1,1]);2.2基于AMDF法的基音周期估计如果信号时周期信号，则相距为周期的整数倍上的抽样点的幅度值是相等的，二者差值为零。对于浊音信号，在基音周期整数倍上，这个值不是零，因此可以通过计算短时平均幅度差函数中相邻谷值间的距离进行基音周期的估计。10()()()Nnmrkxmxmk[b,count]=wavread('*.wav')b1=b(1:640);N=320;%选择的窗长A=[];fork=1:320;sum=0;form=1:N;sum=sum+abs(b1(m)-b1(m+k-1));endA(k)=sum;ends=b(1:320)figure(1)subplot(2,1,1)plot(s);xlabel('样点数')ylabel('幅度')title('浊音信号')axis(