基于MATLAB的语音信号滤波处理.doc

基于MATLAB的语音信号滤波处理

题 目： 基于MATLAB的语音信号滤波处理
课 程： 数字信号处理
学 院： 电气工程学院
班 级：
学 生：
指导教师：

二O一三年十二月

目录CONTENTS

摘要
一、引言
二、正文
1.设计要求
2.设计步骤
3.设计内容
4.简易GUI设计

三、结论
四、收获与心得
五、附录

一、引言
随着Matlab仿真技术的推广，我们可以在计算机上对声音信号进行处理，甚至是模拟。通过计算机作图，采样，我们可以更加直观的了解语音信号的性质，通过matlab编程，调用相关的函数，我们可以非常方便的对信号进行运算和处理。

二、正文
2.1 设计要求
在有噪音的环境中录制语音，并设计滤波器去除噪声。
2.2 设计步骤
1. 分析原始信号，画出原始信号频谱图及时频图，确定滤波器类型及相关指标 ；
2. 按照类型及指标要求设计出滤波器，画出滤波器幅度和相位响应，分析该滤波器是否符合要求；
3. 用所设计的滤波器对原始信号进行滤波处理，画出滤波后信号的频谱图及时频图；
4. 对滤波前的信号进行分析比对，评估所设计滤波器性能。
2.3 设计内容
1.原始信号分析

分析信号的谱图可知，噪音在1650HZ和3300HZ附近的能量较高，而人声的能量基本位于1000HZ以下。因此，可以设计低通滤波器对信号进行去噪处理。
2.IIR滤波器设计
用双线性变换法分别设计了巴特沃斯低通滤波器和椭圆低通滤波器和带阻滤波器：
①巴特沃斯滤波器
fp=800;fs=1300;rs=35;rp=0.5;
程序代码如下：
fp=800;fs=1300;rs=35;rp=0.5;Fs=44100;
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=buttord(wp,ws,rp,rs,'s');
[b,a]=butter(n,wn,'s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);

②椭圆低通滤波器
fp=1300;fs=1600;rs=60;rp=0.5;
程序代码如下：
fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=ellipord(wp,ws,rp,rs,'s');
[b,a]=ellip(n,rp,rs,wn,'s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);

③带阻滤波器
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6
fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6
程序代码如下：
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6;Fs=44100;
fp=[fp1,fp2];fs=[fs1,fs2];
wp=2*Fs*tan(2*pi*fp/(2*Fs));
ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=buttord(wp,ws,rp,rs,'s');
[b,a]=butter(n,wn,'stop','s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);

fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6;Fs=44100;
fp1=[fp3,fp4];fs1=[fs3,fs4];
wp1=2*Fs*tan(2*pi*fp1/(2*Fs));
ws1=2*Fs*tan(2*pi*fs1/(2*Fs));
[n1,wn1]=buttord(wp1,ws1,rp,rs,'s');
[b1,a1]=butter(n1,wn1,'stop','s');
[num1,den1]=bilinear(b1,a1,Fs);
[h1,w1]=freqz(num1,den1,512,Fs);

3. FIR滤波器
①加hamming窗
n=100;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),Hamming(n+1));
[h,w]=freqz(b,1,512,Fs);

②加hanning窗
n=;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),Hanning(n+1));
[h,w]=freqz(b,1,512,Fs);

③加blackman窗
n=100;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),blackman(n+1));
[h,w]=freqz(b,1,512,Fs);

4. 滤波前后比对
①巴特沃斯低通滤波器滤波后

②椭圆低通滤波器滤波后

③带阻滤波器

④加hamming窗

⑤加hanning窗

⑥加blackman窗

2.4简易GUI界面设计
为了便于操作和演示，设计了如下的简易GUI界面。

三、 结论
由以上谱图分析可知，经过滤波器滤波后，信号中的高频杂音明显被抑制，而人声成分大部分被保留，起到了预期的滤波作用。
对比所设计的两种滤波器，椭圆滤波器在过渡带相对较窄的情况下，能满足相对较高阻带衰减。
四、收获与心得
本次设计大概进行了一周的时间，语音信号处理的是目前比较流行且十分有趣的，在编程实现的过程中还是遇到了很多困难。在前期的准备工作中，查阅了大量资料，以完善我们的理论知识。我们为了完成本次设计，我们通过查阅相关书籍以及matlab中的帮助，选用不同的matlab函数，尝试不同的参数。经过接近一个礼拜的反复调试，最终基本的实现了设计任务。
虽然遇到了很多困难，但是我们在设计过程中都有收获很大。本次设计将信号与系统课上学习的知识用于实践，让我们对对语音信号处理更深入的了解，也让我们加深了对滤波器相关内容的理解，同时也使得我们的Matlab能力有了很大的提高。

参考文献

《应用matlab实现信号分析和处理》 科学出版社
附录
1 巴特沃斯低通滤波器
fp=800;fs=1300;rs=35;rp=0.5;Fs=44100;
wp=2*Fs*tan(2*pi*fp/(2*Fs));
ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=buttord(wp,ws,rp,rs,'s');
[b,a]=butter(n,wn,'s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);
figure(1)
%subplot(3,1,1)
plot(w,abs(h));
xlabel('频率/Hz');ylabel('幅值');
title('巴特沃斯低通滤波器幅度特性');
axis([0,5000,0,1.2])
grid on;
figure(2)
%subplot(3,1,2)
plot(w,20*log10(abs(h)));
xlabel('频率/Hz');ylabel('幅值db');
title('巴特沃斯低通滤波器幅度特性db');
axis([0,5000,-90,10]);
grid on;
figure(3)
plot(w,180/pi*unwrap(angle(h)));
xlabel('频率/Hz');ylabel('相位');
title('巴特沃斯低通滤波器相位特性');
axis([0,5000,-1000,10])
grid on;

[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1);
f1=Fs*(0:N1-1)/N1;
t1=(0:N1-1)/Fs;
figure(4)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])

y=filter(num,den,x1);

sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2);
f2=Fs*(0:N2-1)/N2;
t2=(0:N2-1)/Fs;
figure(5)
plot(f2,abs(Y2))
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

2椭圆低通滤波器

fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
wp=2*Fs*tan(2*pi*fp/(2*Fs));
ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=ellipord(wp,ws,rp,rs,'s');
[b,a]=ellip(n,rp,rs,wn,'s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);
figure(1)
plot(w,abs(h));
xlabel('频率/Hz');ylabel('幅值');
title('椭圆低通滤波器幅度特性');
axis([0,5000,0,1.2])
grid on;
figure(2)
plot(w,20*log10(abs(h)));
xlabel('频率/Hz');ylabel('幅值db');
title('椭圆低通滤波器幅度特性db');
axis([0,5000,-90,10]);
grid on;
figure(3)
plot(w,180/pi*unwrap(angle(h)));
xlabel('频率/Hz');ylabel('相位');
title('椭圆低通滤波器相位特性');
axis([0,5000,-1000,10])
grid on;

[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1); %对信号做N点FFT变换
f1=Fs*(0:N1-1)/N1;
t1=(0:N1-1)/Fs;
figure(4)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])

y=filter(num,den,x1);

sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2); %对信号做N点FFT变换
f2=Fs*(0:N2-1)/N2;
t2=(0:N2-1)/Fs;
figure(5)
plot(f2,abs(Y2))
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

3.带阻滤波器
fp1=800;fp2=2300;fs1=1300;fs2=1800;rs=30;rp=0.6;Fs=44100;
fp=[fp1,fp2];fs=[fs1,fs2];
wp=2*Fs*tan(2*pi*fp/(2*Fs));
ws=2*Fs*tan(2*pi*fs/(2*Fs));%wap=2*tan(wp/2)/Ts
[n,wn]=buttord(wp,ws,rp,rs,'s');
[b,a]=butter(n,wn,'stop','s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);

fp3=2800;fp4=4000;fs3=3200;fs4=3700;rs=30;rp=0.6;Fs=44100;
fp1=[fp3,fp4];fs1=[fs3,fs4];
wp1=2*Fs*tan(2*pi*fp1/(2*Fs));
ws1=2*Fs*tan(2*pi*fs1/(2*Fs));%wap=2*tan(wp/2)/Ts
[n1,wn1]=buttord(wp1,ws1,rp,rs,'s');
[b1,a1]=butter(n1,wn1,'stop','s');
[num1,den1]=bilinear(b1,a1,Fs);
[h1,w1]=freqz(num1,den1,512,Fs);
figure(1)
plot(w,abs(h));
xlabel('频率/Hz');ylabel('幅值');
title('巴特沃斯带阻滤波器幅度特性');
axis([0,5000,0,1.2])
grid on;
figure(2)
plot(w,20*log10(abs(h)));
xlabel('频率/Hz');ylabel('幅值db');
title('巴特沃斯带阻滤波器幅度特性db');
axis([0,5000,-90,10]);
grid on;
figure(3)
plot(w,180/pi*unwrap(angle(h)));
xlabel('频率/Hz');ylabel('相位');
title('巴特沃斯带阻滤波器相位特性');
axis([0,5000,-1000,10])
grid on;

figure(4)
plot(w1,abs(h1));
xlabel('频率/Hz');ylabel('幅值');
title('巴特沃斯带阻滤波器幅度特性');
axis([0,5000,0,1.2])
grid on;
figure(5)
plot(w1,20*log10(abs(h1)));
xlabel('频率/Hz');ylabel('幅值db');
title('巴特沃斯带阻滤波器幅度特性db');
axis([0,5000,-90,10]);
grid on;
figure(6)
plot(w1,180/pi*unwrap(angle(h1)));
xlabel('频率/Hz');ylabel('相位');
title('巴特沃斯带阻滤波器相位特性');
axis([0,5000,-1000,10])
grid on;

[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1); %对信号做N点FFT变换
f1=Fs*(0:N1-1)/N1;
t1=(0:N1-1)/Fs;
figure(7)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])

y1=filter(num,den,x1);
y=filter(num1,den1,y1);

sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2);
f2=Fs*(0:N2-1)/N2;
t2=(0:N2-1)/Fs;
figure(8)
plot(f2,abs(Y2))
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

4. 加hamming窗
n=100;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),Hamming(n+1));
[h,w]=freqz(b,1,512,Fs);
a=num2str(a);
[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);
sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1);
f1=Fs*(0:N1-1)/N1;
y=fftfilt(b,x1);
sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2);
f2=Fs*(0:N2-1)/N2;
figure(4)
subplot(2,1,1)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])
subplot(2,1,2)
plot(f2,abs(Y2));
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

5. 加hanning窗
n=100;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),Hanning(n+1));
[h,w]=freqz(b,1,512,Fs);
a=num2str(a);
[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);
sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1);
f1=Fs*(0:N1-1)/N1;
y=fftfilt(b,x1);
sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2);
f2=Fs*(0:N2-1)/N2;
figure(4)
subplot(2,1,1)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])
subplot(2,1,2)
plot(f2,abs(Y2));
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

6. 加blackman窗
n=100;fp=1000;Fs=44100;
b=fir1(n,fp/(Fs/2),blackman(n+1));
[h,w]=freqz(b,1,512,Fs);
a=num2str(a);
[s1,Fs,bits]=wavread('D:\222.wav');
x1=s1(:,1);
sound(x1,Fs,bits);
N1=length(x1);
Y1=fft(x1,N1);
f1=Fs*(0:N1-1)/N1;
y=fftfilt(b,x1);
sound(y,Fs,bits);
N2=length(y);
Y2=fft(y,N2);
f2=Fs*(0:N2-1)/N2;
figure(4)
subplot(2,1,1)
plot(f1,abs(Y1))
xlabel('频率/Hz');ylabel('幅度');
title('原始信号频谱');
grid on;axis([0 6000 0 400])
subplot(2,1,2)
plot(f2,abs(Y2));
xlabel('频率/Hz');ylabel('幅度');
title('过滤后信号的频谱');
grid on;axis([0 6000 0 100])

7.巴特沃斯低通滤波时频分析

clear
Fs=44100;
[s1,Fs,bits]=wavread('D:\222.wav');
x=s1(:,1);x1=x(40001:64576)
l=length(x1);
N=6;
m=downsample(x1,N);%降抽样后 Fs=7350hz
z=length(m);%4096
[tfr, t, f] = tfrstft(m);
figure(13)
contour(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
xlabel('时间t');
ylabel('频率f');
title('等高线图（滤波前）');
grid on;
figure(14)
mesh(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
xlabel('时间t');
ylabel('频率f');
zlabel('幅值A');
title('三维图（滤波前）');
grid on;

fp=1300;fs=1600;rs=60;rp=0.5;Fs=44100;
wp=2*Fs*tan(2*pi*fp/(2*Fs));ws=2*Fs*tan(2*pi*fs/(2*Fs));
[n,wn]=ellipord(wp,ws,rp,rs,'s');
[b,a]=ellip(n,rp,rs,wn,'s');
[num,den]=bilinear(b,a,Fs);
[h,w]=freqz(num,den,512,Fs);

y=filter(num,den,x);
x1=y(40001:64576)
l=length(x1);
N=6;
m=downsample(x1,N);%降抽样后 Fs=7350hz
z=length(m);%4096
[tfr,t,f] = tfrstft(m);
figure(15)
contour(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
xlabel('时间t');
ylabel('频率f');
title('等高线图（滤波后）');
grid on;
figure(16)
mesh(t,(Fs/N)*(0:z/2-1)/z,abs(tfr(1:z/2,:)).^2);
xlabel('时间t');
ylabel('频率f');
zlabel('幅值A');
title('三维图（滤波后）');
grid on;

33