语音信号处理

1、数字信号处理语音信号处理一、 任务要求选择一个语音信号作为分析的对象，或录制一段各人自己的语音信号，对其进行频谱分析；利用MATLAB中的随机函数产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析；设计FIR和IIR数字滤波器， 并对被噪声污染的语音信号进行滤波， 分析滤波后信号的时域和频域特征，回放语音信号。二、 语音信号处理总流程信号的采集语音信号分析含噪信号合成数字滤波器设计滤波结果显示分析图1 语音信号处理总流程三、 原始信号采集及频谱分析1 语音信号的读入与打开在MATLAB中，y,fs,bits=wavread(' E:dwje.wav ');用于读取

2、语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。可以使用sound(y,fs,bits)，用于对声音的回放。向量y则就代表了一个信号，也即一个复杂的“函数表达式”，也可以说像处理一个信号的表达式一样处理这个声音信号。2 语音信号的频频分析利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。图2 原始信号时域和频域分析对语音信号进行FFT变换就是为了得到它的频域的图形，便于从图中观察出信号的幅度等特性，从图中可以看出语音信号的截止频率为5000Hz。四、 语音信号加噪与频谱分析在Matlab中人为设计一个固定频率6000Hz的余弦序列噪声干扰信号d=0.0

3、5*cos(2*pi*6000*t)'。然后将噪声与原始语音信号叠加Y=y+d，叠加两信号的维数必须相同，否则无法叠加。对噪音信号进行频谱变换得到其频谱图,从图中可以看出干扰信号，在6000Hz频点处有一高峰,其中5500Hz 正是本设计所要利用的。图3 叠加后信号时域和频域分析五、 设计IIR滤波器对加噪语音滤波1 利用模拟滤波器设计IIR数字低通滤波器的设计流程如图4所示：确定技术指标技术指标转换：数字到模拟双线性变换法：W=2/T *tan(w/2)实现模拟滤波器设计将模拟滤波器转换为数字滤波器图4 IIR滤波器设计流程2 滤波器参数选取由图3观察可得，通带截止频率Wp=5000

4、Hz，阻带截止频率Ws=5500Hz，一般情况下通带最大衰减Rp=0.1，阻带最小衰减Rs=40。根据设计流程设计IIR滤波器如图5所示：图5 IIR低通滤波器3 滤波前后波形分析图6 原始信号与滤波后信号波形比较图7 原始信号与滤波后信号频谱比较由图6和图7可以看出滤波前后语音信号的波形和频谱几乎一样，所以该滤波器设计合格。六、 设计FIR滤波器对加噪语音滤波1 利用窗函数设计FIR数字低通滤波器的设计流程如图8所示：图8 FIR滤波器设计流程2 窗函数类型选取常用的窗函数一共有矩形窗、汉宁窗、海明窗、凯泽窗。这几种窗函数的基本参数比较如表1所示，根据实际阻带最小衰减选择合适的窗函数设计FI

5、R滤波器。表1 几种窗函数的基本参数比较最小阻带衰减只由窗形决定，不受N的影响，而过渡带宽则随N的增加而减小。3 窗长度N的确定通过求得的过渡带宽以及matlab所提供的ceil()函数，可以很方便的求得窗口长度N，N=ceil(6.6*pi/Tr_width)。4 实际滤波器的设计通过Matlab中设计的函数ideallp()求得理想冲激响应，即hd=ideallp(Wc,N)。通过窗长度N，计算设计理想低通的hdn，求出汉宁窗序列。然后给序列加窗hn=hdn.*wn。计算频率响应：H,w=freqz(h,1,512)。设计完成后滤波器如下图9所示：图9 FIR滤波器5 滤波前后语音信号对比

6、分析滤波前后原始语音信号波形的波形和频谱比较如图10和图11所示，所以该波器设计符合题目要求。图10 FIR滤波后语音和原语音信号波形比较图11 FIR滤波后语音和原语音信号频谱比较附录：（源程序代码）clc;clear;fs=22050,bits=1024;y,fs,bits=wavread('E:dwje.wav'); %y:语音数据；fs:采样频率；bits:采样点数% sound(y,fs,bits); %话音回放figure(1); %建立subplot(2,1,1); %对图形窗口进行分割plot(y);xlabel('time n');ylabel

7、('amplitude');title('原始信号波形');y1=fft(y,1024); %fft就是把时域上的曲线变成频域上的曲线，所以fft变换后横坐标是频率fx=fs*(0:511)/1024; %将横轴变为频率轴fy=abs(y1(1:512); %将纵轴变为频率幅度轴subplot(2,1,2);plot(fx,fy);xlabel('频率Hz');ylabel('幅度');title('原始信号频谱');t=(0:length(y)-1)/fs;d=0.05*cos(2*pi*6000*t)'

8、Y=y+d; % sound(Y,fs,bits); figure(2);subplot(2,1,1);plot(Y); xlabel('time n');ylabel('amplitude');title('叠噪信号时域波形');y2=fft(Y,1024); %将横轴变为频率轴fy=abs(y2(1:512); %将纵轴变为频率幅度轴subplot(2,1,2);plot(fx,fy);xlabel('频率Hz');ylabel('幅度');title('叠噪信号频谱');Ts=1/fs;Wp=

9、2*pi*5000/fs;Ws=2*pi*5500/fs;Rp=0.1;Rs=40; %通带最大衰减(db)，阻带最大衰减(db) Wp1=2/Ts*tan(Wp/2);Ws1=2/Ts*tan(Ws/2); %将模拟指标转换成数字指标N,Wn=buttord(Wp1,Ws1,Rp,Rs,'s'); %选择滤波器的最小阶数(估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc)Z,P,K=buttap(N); %创建Butterworth低通滤波器原型Bap,Aap=zp2tf(Z,P,K); %将零极点增益转换为普遍分子，分母b,a=lp2lp(Bap,A

10、ap,Wn); %将普遍的分子和分母转换为以Wn为截止频率bz,az=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的H,W=freqz(bz,az); %求频率响应figure(3);plot(W*fs/(2*pi),abs(H); %绘制Butterworth低通滤波器频率响应曲线grid; %添加图格xlabel('频率/Hz');ylabel('频率响应幅度');title('IIR低通滤波器')f1=filter(bz,az,y);figure(4);subplot(2,1,1);plot(y);titl

11、e('IIR滤波前的时域波形');subplot(2,1,2);plot(f1);title('IIR滤波后的时域波形');% sound(f1,22050);F0=fft(f1,1024);fx=fs*(0:511)/1024; %将横轴变为频率轴 /此处出错fy=abs(y1(1:512); %将纵轴变为频率幅度轴figure(5)subplot(2,1,1); %创建两行一列绘图区间的第1绘图区间plot(fx,fy);xlabel('频率Hz');ylabel('幅度');title('原始信号频谱');s

12、ubplot(2,1,2); %创建两行一列绘图区间的第2绘图区间 plot(fx,abs(F0(1:512);xlabel('频率');ylabel('幅值');title('滤波后的频谱');%窗函数法设计滤波器fs=22050; Wp=2*pi*5000/fs;Ws=2*pi*5500/fs;rp=0.1;rs=40;%指标参数 Bt=Ws-Wp; %过渡带宽度 N=ceil(6.2*pi/Bt); %由阻带最小衰减指标确定选hanning窗，并确定其阶数 wc=(Wp+Ws)/2;r=(N-1)/2; %理想低通截止频率wc n=0:N-

13、1;hdn=sin(wc*(n-r)./(pi*(n-r); %计算理想低通的hdn window=hanning(N+1); hdn(16)=wc/pi; %在n=(N-1)/2=15点为0/0型，直接赋值 wn=0.5*(1-cos(2*pi*n/(N-1); %求窗序列wn hn=hdn.*wn; %加窗 H,w=freqz(window,1); figure(6); subplot(2,2,1); stem(window); xlabel('频率Hz');ylabel('幅度');title('Hanning窗函数');grid on;

14、subplot(2,2,2);plot(w/pi,20*log(abs(H)/abs(H(1);axis(0 1 -400 0);xlabel('w/pi');ylabel('幅度(db)');title('Hanning窗函数的频谱');grid;hn1=fir1(N,wc/pi,window);h2,w=freqz(hn1,1,512);subplot(2,2,3);stem(hn1);xlabel('n');ylabel('hn');title('Hanning窗函数的单位脉冲响应');grid;subplot(2,2,4);plot(w/pi,20*log(abs(h2)/abs(h2(1);xlabel('w/pi');ylabe