MATLAB的有噪声的语音信号分析与处理设计武汉理工大学

1、 录目 4）1课程设计目的（ 4）.（2课程设计基本要求 4）.（3课程设计内容. ）4课程设计实现（4 . ）（4.（1）语音信号的采集 ）（52）语音信号的时频分析.（ 5）语音信号加噪与频谱分析.（3 ）.（8IIR（4）设计FIR和数字滤波器 . 5）用滤波器对加噪语音信号进行滤波（ ）比较滤波前后语音信号的波形及频谱6（ . ）回放语音信号（7 11）.、心得体会5.（ ）11（.、参考文献6 的有噪声的语音信号分析与处理设计基于MATLAB 课程设计目的1综合运用数字信号处理的理论知识对语音信号进行时频分析和滤波器设计，从作为编程工具进行计算机实现，通过理论推导得出相应结论，再利用

2、 MATLAB 而加深对所学知识的理解，建立概念。 课程设计基本要求2 的程序设计方法。的使用，掌握 MATLAB 学会 1 ) MATLAB 环境下语音信号采集的方法。掌握在 Windows 2 ) 掌握数字信号处理的基本概念、基本理论和基本方法。 3 ) 数字滤波器的方法。和IIR 4 ) 掌握 MATLAB 设计 FIR 对信号进行分析和处理。学会用 MATLAB 5 ) 课程设计内容 31）选择一个语音信号作为分析对象，或录制一段语音信号； 2）对语音信号进行采样，画出采样后语音信号的时域波形和频谱图； 3）利用MATLAB中的随机函数产生噪声加入到语音信号中，使语音信号被污染，然后进

3、行频谱分析； 4）设计FIR和IIR数字滤波器，并对被噪声污染的语音信号进行滤波，画出滤波前后信号的时域波形和频谱，并对滤波前后的信号进行比较，分析信号的变化； 5）回放语音信号。 4课程设计实现 （1）语音信号的采集 利用PC 机上的声卡和WINDOWS 操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插口上,启动录音机。按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名“speech”保存入g : MATLAB work 中。可以看到,文件存储器的后缀默认为. wav ,这是WINDOWS 操作

4、系统规定的声音文件存的标准。 （2）语音信号的时频分析 Matlab软件平台下，利用wavread函数对语音信号进行采样，记住采样频率和采样点数 Wavread 函数调用格式 y=wavread（file）%读取file所规定的wav文件，返回采样值放在向量y中。 y,fs,nbits=wavread(file) %采样值放在向量y中，fs表示采样频率（hz），nbits表示采样位数。 中。y点的采样值放在向量N读取钱%）N，file（y=wavready=wavread（file，N1,N2）%读取从N1到N2点的采样值放在向量y中。 对语音信号speech.wav进行采样其程序如下： y,

5、fs,nbits=wavered (speech); %把语音信号进行加载入Matlab 仿真软件平台中 fs = 8000 nbits = 16 首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在matlab中利用fft对信号进行快速傅里叶变换，得到信号的频谱特性。 其程序如下： y,fs,nbits=wavread (speech); sound(y,fs,nbits); %回放语音信号 n = length (y) ; %求出语音信号的长度 Y=fft(y,n); %傅里叶变换 subplot(2,1,1);plot(y);title(原始信号波形); subplot(2,1,2

6、);plot(abs(Y);title(原始信号频谱) 程序结果如下图： （3）语音信号加噪与频谱分析 利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。 其程序如下： y,fs,nbits=wavread (speech); n = length (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 sound(s); subplot(2,1,1); plot(s);title(加噪语音信号的时域波形); S=fft(s); %傅里叶变换 s

7、ubplot(2,1,2); plot(abs(S);title(加噪语音信号的频域波形) 程序结果如下图： (4) 设计FIR和IIR数字滤波器 根据语音信号的特点给出有关滤波器的新能指标： 低通滤波器的性能指标 fp=1000Hz，fc=1200Hz，As=100db ,Ap=1dB 高通滤波器的性能指标 fp=4800Hz，fc=5000Hz，As=100dB，Ap=1dB； 带通滤波器的性能指标 fp1=1200Hz，fp2=3000hZ，fc1=1000Hz，fc2=3200Hz， As=100dB，Ap=1dB 在Matlab中，可以利用函数fir1设计FIR滤波器，利用函数but

8、ter,cheby1和ellip设计IIR滤波器，利用Matlab中的函数freqz画出各步步器的频率响应。 分析如下：函数fir1默认的设计滤波器的方法为窗函数法，其中可选的窗函数有Rectangular Barlrtt Hamming Hann Blackman窗，其相应的都有实现函数。 函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。 IIR低通： Ft=8000; Fp=1000; Fs=1200; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fp=2*Ft*tan(wp/2); fs=

9、2*Fs*tan(wp/2); n11,wn11=buttord(wp,ws,1,50,s); %求低通滤波器的阶数和截止频率 b11,a11=butter(n11,wn11,s); %求S域的频率响应的参数 num11,den11=bilinear(b11,a11,0.5); %利用双线性变换实现频率响应S域到Z域的变换 h,w=freqz(num11,den11); %根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h); legend(用butter设计); grid IIR带通： Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; Ft=

10、8000; wp1=tan(pi*Fp1/Ft); %带通到低通滤波器的转换 wp2=tan(pi*Fp2/Ft); ws1=tan(pi*Fs1/Ft); ws2=tan(pi*Fs2/Ft); w=wp1*wp2/ws2; bw=wp2-wp1; wp=1; ws=(wp1*wp2-w.2)/(bw*w); n12,wn12=buttord(wp,ws,1,50,s); %求低通滤波器阶数和截止频率 b12,a12=butter(n12,wn12,s); %求S域的频率响应参数 num2,den2=lp2bp(b12,a12,sqrt(wp1*wp2),bw);%将S域低通参数转为带通的

11、num12,den12=bilinear(num2,den2,0.5);%利用双线性变换实现频率响应S域到Z域的转换 h,w=freqz(num12,den12);%根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h); axis(0 4000 0 1.5); legend(用butter设计); grid; IIR高通： Ft=8000; Fp=4000; Fs=3500; 高通到低通滤波器参数转换wp1=tan(pi*Fp/Ft);%ws1=tan(pi*Fs/Ft); wp=1; ws=wp1*wp/ws1; 求模拟的低通滤波器阶数和截止频率n13,wn13=cheb

12、1ord(wp,ws,1,50,s); % S域的频率响应的参数 b13,a13=cheby1(n13,1,wn13,s); %求 域低通参数转为高通的Snum,den=lp2hp(b13,a13,wn13);%将域ZS域到利用双线性变换实现频率响应num13,den13=bilinear(num,den,0.5); % 转换h,w=freqz(num13,den13); plot(w*21000*0.5/pi,abs(h); ); 高通滤波器title(IIR); 设计cheby1legend(用axis(0 12000 0 1.5); grid; FIR低通 :用窗函数设计低通滤波器的程序

13、如下Ft=8000; Fp=1000; Fs=1200; wp=2*Fp/Ft; ws=2*Fs/Ft; rp=1; rs=50; p=1-10.(-rp/20); %通带阻带波纹 s=10.(-rs/20); fpts=wp ws; mag=1 0; dev=p s; n21,wn21,beta,ftype=kaiserord(fpts,mag,dev);%由kaiserord求滤波器的阶数和截止频率 b21=fir1(n21,wn21,Kaiser(n21+1,beta); %由fir1设计滤波器 h,w=freqz(b21,1); %得到频率响应 plot(w/pi,abs(h); ti

14、tle(FIR低通滤波器); grid; 带通滤波器：FIRFp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; Ft=8000; %带通到低通滤波器参数转换wp1=tan(pi*Fp1/Ft); wp2=tan(pi*Fp2/Ft); ws1=tan(pi*Fs1/Ft); ws2=tan(pi*Fs2/Ft); w=wp1*wp2/ws2; bw=wp2-wp1; wp=1; ws=(wp*wp2-w.2)/(bw*w); 求低通滤波器阶数和截止频率 %n22,wn22=buttord(wp,ws,1,50,s); S求域的频率响应的参数 %b22,a22=but

15、ter(n22,wn22,s); S域低通参数转为带通的num2,den2=lp2bp(b22,a22,sqrt(wp1*wp2),bw); %将 域的转换S域到Znum22,den22=bilinear(num2,den2,0.5);%利用双线性变换实现 根据参数求出频率响应 % h,w=freqz(num22,den22); plot(w*8000*0.5/pi,abs(h); axis(0 4000 0 1.5); ); butterlegend(用设计grid FIR高通： Ft=8001; Fp=4000; Fs=3500; wp=2*Fp/Ft; ws=2*Fs/Ft; rp=1;

16、 rs=50; p=1-10.(-rp/20); %通带阻带波纹 s=10.(-rs/20); fpts=ws wp; mag=0 1; dev=p s; n23,wn23,beta,ftype=kaiserord(fpts,mag,dev); b23=fir1(n23,wn23,high,Kaiser(n23+1,beta); %由fir1设计滤波器 h,w=freqz(b23,1); %得到频率响应 plot(w*12000*0.5/pi,abs(h); title(FIR高通滤波器); axis(3000 6000 0 1.2); grid （5）用滤波器对加噪语音信号进行滤波 用自己设

17、计的各滤波器分别对加噪的语音信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。 双线性变换法： i.低通滤波器 z11=filter(num11,den11,s); ii.带通滤波器 z12=filter(num12,den12,s); iii.高通滤波器 z13=filter(num13,den13,s); 窗函数法 i.低通滤波器 z21=fftfilt(b21,s); ii.带通滤波器 z22=fftfilt(b22,s); iii.高通滤波器 z23=fftfilt(b23,s); （ 6）比较滤波前后语音

18、信号的波形及频谱 双线性变换法 i.低通滤波器 y,fs,nbits=wavread (speech); n = length (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z11=filter(num11,den11,s); sound(z11); m11=fft(z11); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subplot(2,2,2); plot(abs(m11),r

19、); title(滤波后信号的频谱); grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z11); title(滤波后的信号波形); grid; ii.带通滤波器 y,fs,nbits=wavread (speech); n = length (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z12=filter(num12,den12,s); sound(z12);

20、m12=fft(z12); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subplot(2,2,2); plot(abs(m12),r); title(滤波后信号的频谱); grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z12); title(滤波后的信号波形); grid; 高通滤波器iii.y,fs,nbits=wavread (speech); 求出语音信号的长度 % n = length (y) ;

21、随机函数产生噪声 %noise=0.01*randn(n,2); %语音信号加入噪声 s=y+noise; 傅里叶变换 % S=fft(s); z13=filter(num13,den13,s); sound(z13); 求滤波后的信号 % m13=fft(z13); subplot(2,2,1); plot(abs(S),g); ); title(滤波前信号的频谱grid; subplot(2,2,2); plot(abs(m13),r); ); 滤波后信号的频谱title(grid; subplot(2,2,3); plot(s); ); 滤波前信号的波形title(grid; subpl

22、ot(2,2,4); plot(z13); title(滤波后的信号波形); grid; 窗函数法 i.低通滤波器 y,fs,nbits=wavread (speech); n = length (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z21=fftfilt(b21,s); sound(z21); m21=fft(z21); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subp

23、lot(2,2,2); plot(abs(m21),r); title(滤波后信号的频谱); grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z21); title(滤波后的信号波形); grid; ii.带通滤波器 y,fs,nbits=wavread (speech); n = length (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z22=fftfilt(

24、b22,s); sound(z22); m22=fft(z22); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subplot(2,2,2); plot(abs(m22),r); title(滤波后信号的频谱); grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z22); title(滤波后的信号波形); grid; iii.高通滤波器 y,fs,nbits=wavread (speech); n = leng

25、th (y) ; %求出语音信号的长度 noise=0.01*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z23=fftfilt(b23,s); sound(z23); m23=fft(z23); %求滤波后的信号 subplot(2,2,1); plot(abs(S),g); title(滤波前信号的频谱); grid; subplot(2,2,2); plot(abs(m23),r); title(滤波后信号的频谱); grid; subplot(2,2,3); plot(s); title(滤波前信号的波形); grid; subplot(2,2,4); plot(z23); title(滤波后的信号波形); grid; 分析：加入噪声后回放的声音与原始的语音信号有明显的不同,其伴随较尖锐的干扰啸叫声。从含噪语音